Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Vector/BuildTestProblems.cpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Epetra: Linear Algebra Services Package
5// Copyright 2011 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42
43#include "BuildTestProblems.h"
44
46 const char TransA, const char TransB,
47 const double alpha,
50 const double beta,
51 Epetra_Vector& C_GEMM ) {
52
53 // For given values of TransA, TransB, alpha and beta, a (possibly
54 // zero) filled Epetra_Vector C, and allocated
55 // Epetra_Vectors A, B and C_GEMM this routine will generate values for
56 // Epetra_Vectors A, B and C_GEMM such that, if A, B and (this) are
57 // used with GEMM in this class, the results should match the results
58 // generated by this routine.
59
60 // Test for Strided multivectors (required for GEMM ops)
61
62 if (!A.ConstantStride() ||
63 !B.ConstantStride() ||
64 !C_GEMM.ConstantStride() ||
65 !C.ConstantStride()) return(-1); // Error
66
67 int i, j;
68 double fi, fj; // Used for casting loop variables to floats
69
70 // Get a view of the Vectors
71
72 double *Ap = 0;
73 double *Bp = 0;
74 double *Cp = 0;
75 double *C_GEMMp = 0;
76
77 int A_nrows = A.MyLength();
78 int A_ncols = A.NumVectors();
79 int B_nrows = B.MyLength();
80 int B_ncols = B.NumVectors();
81 int C_nrows = C.MyLength();
82 int C_ncols = C.NumVectors();
83 int A_Stride = 0;
84 int B_Stride = 0;
85 int C_Stride = 0;
86 int C_GEMM_Stride = 0;
87
88 A.ExtractView(&Ap);
89 B.ExtractView(&Bp);
90 C.ExtractView(&Cp);
91 C_GEMM.ExtractView(&C_GEMMp);
92
93 // Define some useful constants
94
95
96 int opA_ncols = (TransA=='N') ? A.NumVectors() : A.MyLength();
97 int opB_nrows = (TransB=='N') ? B.MyLength() : B.NumVectors();
98
99 int C_global_inner_dim = (TransA=='N') ? A.NumVectors() : A.GlobalLength();
100
101
102 bool A_is_local = (!A.DistributedGlobal());
103 bool B_is_local = (!B.DistributedGlobal());
104 bool C_is_local = (!C.DistributedGlobal());
105
106 int A_IndexBase = A.Map().IndexBase();
107 int B_IndexBase = B.Map().IndexBase();
108
109 // Build two new maps that we can use for defining global equation indices below
110 Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
111 Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());
112
113 int* A_MyGlobalElements = new int[A_nrows];
114 A_Map->MyGlobalElements(A_MyGlobalElements);
115 int* B_MyGlobalElements = new int[B_nrows];
116 B_Map->MyGlobalElements(B_MyGlobalElements);
117
118 // Check for compatible dimensions
119
120 if (C.MyLength() != C_nrows ||
121 opA_ncols != opB_nrows ||
122 C.NumVectors() != C_ncols ||
123 C.MyLength() != C_GEMM.MyLength() ||
124 C.NumVectors() != C_GEMM.NumVectors() ) {
125 delete A_Map;
126 delete B_Map;
127 delete [] A_MyGlobalElements;
128 delete [] B_MyGlobalElements;
129 return(-2); // Return error
130 }
131
132 bool Case1 = ( A_is_local && B_is_local && C_is_local); // Case 1 above
133 bool Case2 = (!A_is_local && !B_is_local && C_is_local && TransA=='T' );// Case 2
134 bool Case3 = (!A_is_local && B_is_local && !C_is_local && TransA=='N');// Case 3
135
136 // Test for meaningful cases
137
138 if (!(Case1 || Case2 || Case3)) {
139 delete A_Map;
140 delete B_Map;
141 delete [] A_MyGlobalElements;
142 delete [] B_MyGlobalElements;
143 return(-3); // Meaningless case
144 }
145
146 /* Fill A, B and C with values as follows:
147
148 If A_is_local is false:
149 A(i,j) = A_MyGlobalElements[i]*j, i=1,...,numLocalEquations, j=1,...,NumVectors
150 else
151 A(i,j) = i*j, i=1,...,numLocalEquations, j=1,...,NumVectors
152
153 If B_is_local is false:
154 B(i,j) = 1/(A_MyGlobalElements[i]*j), i=1,...,numLocalEquations, j=1,...,NumVectors
155
156 else
157 B(i,j) = 1/(i*j), i=1,...,numLocalEquations, j=1,...,NumVectors
158
159 In addition, scale each entry by GlobalLength for A and
160 1/GlobalLength for B--keeps the magnitude of entries in check
161
162
163 C_GEMM will depend on A_is_local and B_is_local. Three cases:
164
165 1) A_is_local true and B_is_local true:
166 C_GEMM will be local replicated and equal to A*B = i*NumVectors/j
167
168 2) A_is_local false and B_is_local false
169 C_GEMM will be local replicated = A(trans)*B(i,j) = i*numGlobalEquations/j
170
171 3) A_is_local false B_is_local true
172 C_GEMM will distributed global and equals A*B = A_MyGlobalElements[i]*NumVectors/j
173
174 */
175
176 // Define a scalar to keep magnitude of entries reasonable
177
178 double sf = C_global_inner_dim;
179 double sfinv = 1.0/sf;
180
181 // Define A depending on A_is_local
182
183 if (A_is_local)
184 {
185 for (j = 0; j <A_ncols ; j++)
186 for (i = 0; i<A_nrows; i++)
187 {
188 fi = i+1; // Get float version of i and j, offset by 1.
189 fj = j+1;
190 Ap[i + A_Stride*j] = (fi*sfinv)*fj;
191 }
192 }
193 else
194 {
195 for (j = 0; j <A_ncols ; j++)
196 for (i = 0; i<A_nrows; i++)
197 {
198 fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
199 fj = j+1;
200 Ap[i + A_Stride*j] = (fi*sfinv)*fj;
201 }
202 }
203
204 // Define B depending on TransB and B_is_local
205
206 if (B_is_local)
207 {
208 for (j = 0; j <B_ncols ; j++)
209 for (i = 0; i<B_nrows; i++)
210 {
211 fi = i+1; // Get float version of i and j, offset by 1.
212 fj = j+1;
213 Bp[i + B_Stride*j] = 1.0/((fi*sfinv)*fj);
214 }
215 }
216 else
217 {
218 for (j = 0; j <B_ncols ; j++)
219 for (i = 0; i<B_nrows; i++)
220 {
221 fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
222 fj = j+1;
223 Bp[i + B_Stride*j] = 1.0/((fi*sfinv)*fj);
224 }
225 }
226 // Define C_GEMM depending on A_is_local and B_is_local. C_GEMM is also a
227 // function of alpha, beta, TransA, TransB:
228
229 // C_GEMM = alpha*A(TransA)*B(TransB) + beta*C_GEMM
230
231 if (Case1)
232 {
233 for (j = 0; j <C_ncols ; j++)
234 for (i = 0; i<C_nrows; i++)
235 {
236 // Get float version of i and j, offset by 1.
237 fi = (i+1)*C_global_inner_dim;
238 fj = j+1;
239 C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
240 + beta * Cp[i + C_Stride*j];
241 }
242 }
243 else if (Case2)
244 {
245 for (j = 0; j <C_ncols ; j++)
246 for (i = 0; i<C_nrows; i++)
247 {
248 // Get float version of i and j, offset by 1.
249 fi = (i+1)*C_global_inner_dim;
250 fj = j+1;
251 C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
252 + beta * Cp[i + C_Stride*j];
253 }
254 }
255 else
256 {
257 for (j = 0; j <C_ncols ; j++)
258 for (i = 0; i<C_nrows; i++)
259 {
260 // Get float version of i and j.
261 fi = (A_MyGlobalElements[i]+1)*C_global_inner_dim;
262 fj = j+1;
263 C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
264 + beta * Cp[i + C_Stride*j];
265 }
266 }
267 delete A_Map;
268 delete B_Map;
269 delete []A_MyGlobalElements;
270 delete []B_MyGlobalElements;
271
272 return(0);
273 }
274int BuildVectorTests (Epetra_Vector & C, const double alpha,
275 Epetra_Vector& A,
276 Epetra_Vector& sqrtA,
277 Epetra_Vector& B,
278 Epetra_Vector& C_alphaA,
279 Epetra_Vector& C_alphaAplusB,
280 Epetra_Vector& C_plusB,
281 double* const dotvec_AB,
282 double* const norm1_A,
283 double* const norm2_sqrtA,
284 double* const norminf_A,
285 double* const normw_A,
286 Epetra_Vector& Weights,
287 double* const minval_A,
288 double* const maxval_A,
289 double* const meanval_A ) {
290
291 // For given values alpha and a (possibly zero) filled
292 // Epetra_Vector (the this object), allocated double * arguments dotvec_AB,
293 // norm1_A, and norm2_A, and allocated Epetra_Vectors A, sqrtA,
294 // B, C_alpha, C_alphaAplusB and C_plusB, this method will generate values for
295 // Epetra_Vectors A, B and all of the additional arguments on
296 // the list above such that, if A, B and (this) are used with the methods in
297 // this class, the results should match the results generated by this routine.
298 // Specifically, the results in dotvec_AB should match those from a call to
299 // A.dotProd (B,dotvec). Similarly for other routines.
300
301 int i,j;
302 double fi; // Used for casting loop variables to floats
303 // Define some useful constants
304
305 int A_nrows = A.MyLength();
306 int A_ncols = A.NumVectors();
307 int sqrtA_nrows = sqrtA.MyLength();
308 int sqrtA_ncols = sqrtA.NumVectors();
309 int B_nrows = B.MyLength();
310 int B_ncols = B.NumVectors();
311
312 double *Ap = 0;
313 double *sqrtAp = 0;
314 double *Bp = 0;
315 double *Cp = 0;
316 double *C_alphaAp = 0;
317 double *C_alphaAplusBp = 0;
318 double *C_plusBp = 0;
319 double *Weightsp = 0;
320
321 A.ExtractView(&Ap);
322 sqrtA.ExtractView(&sqrtAp);
323 B.ExtractView(&Bp);
324 C.ExtractView(&Cp);
325 C_alphaA.ExtractView(&C_alphaAp);
326 C_alphaAplusB.ExtractView(&C_alphaAplusBp);
327 C_plusB.ExtractView(&C_plusBp);
328 Weights.ExtractView(&Weightsp);
329
330 bool A_is_local = (A.MyLength() == A.GlobalLength());
331 bool B_is_local = (B.MyLength() == B.GlobalLength());
332 bool C_is_local = (C.MyLength() == C.GlobalLength());
333
334 int A_IndexBase = A.Map().IndexBase();
335 int B_IndexBase = B.Map().IndexBase();
336
337 // Build two new maps that we can use for defining global equation indices below
338 Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
339 Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());
340
341 int* A_MyGlobalElements = new int[A_nrows];
342 A_Map->MyGlobalElements(A_MyGlobalElements);
343 int* B_MyGlobalElements = new int[B_nrows];
344 B_Map->MyGlobalElements(B_MyGlobalElements);
345
346 // Check for compatible dimensions
347
348 if (C.MyLength() != A_nrows ||
349 A_nrows != B_nrows ||
350 C.NumVectors() != A_ncols ||
351 A_ncols != B_ncols ||
352 sqrtA_nrows != A_nrows ||
353 sqrtA_ncols != A_ncols ||
354 C.MyLength() != C_alphaA.MyLength() ||
355 C.NumVectors() != C_alphaA.NumVectors() ||
356 C.MyLength() != C_alphaAplusB.MyLength() ||
357 C.NumVectors() != C_alphaAplusB.NumVectors() ||
358 C.MyLength() != C_plusB.MyLength() ||
359 C.NumVectors() != C_plusB.NumVectors() ) return(-2); // Return error
360
361
362 bool Case1 = ( A_is_local && B_is_local && C_is_local); // Case 1
363 bool Case2 = (!A_is_local && !B_is_local && !C_is_local);// Case 2
364
365 // Test for meaningful cases
366
367 if (!(Case1 || Case2)) return(-3); // Meaningless case
368
369 /* Fill A and B with values as follows:
370
371 If A_is_local is false:
372 A(i,j) = A_MyGlobalElements[i]*j, i=1,...,numLocalEquations, j=1,...,NumVectors
373
374 else
375 A(i,j) = i*j, i=1,...,numLocalEquations, j=1,...,NumVectors
376
377 If B_is_local is false:
378 B(i,j) = 1/(A_MyGlobalElements[i]*j), i=1,...,numLocalEquations, j=1,...,NumVectors
379
380 else
381 B(i,j) = 1/(i*j), i=1,...,numLocalEquations, j=1,...,NumVectors
382
383 In addition, scale each entry by GlobalLength for A and
384 1/GlobalLength for B--keeps the magnitude of entries in check
385 */
386
387 //Define scale factor
388
389 double sf = A.GlobalLength();
390 double sfinv = 1.0/sf;
391
392 // Define A
393
394 if (A_is_local) {
395 for (i = 0; i<A_nrows; i++) {
396 fi = i+1; // Get float version of i and j, offset by 1.
397 Ap[i] = (fi*sfinv);
398 sqrtAp[i] = std::sqrt(Ap[i]);
399 }
400 }
401 else {
402 for (i = 0; i<A_nrows; i++) {
403 fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
404 Ap[i] = (fi*sfinv);
405 sqrtAp[i] = std::sqrt(Ap[i]);
406 }
407 }
408
409 // Define B depending on TransB and B_is_local
410
411 if (B_is_local) {
412 for (i = 0; i<B_nrows; i++) {
413 fi = i+1; // Get float version of i and j, offset by 1.
414 Bp[i] = 1.0/((fi*sfinv));
415 }
416 }
417 else {
418 for (i = 0; i<B_nrows; i++) {
419 fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
420 Bp[i] = 1.0/((fi*sfinv));
421 }
422 }
423
424 // Generate C_alphaA = alpha * A
425
426 for (i = 0; i<A_nrows; i++) C_alphaAp[i] = alpha * Ap[i];
427
428 // Generate C_alphaA = alpha * A + B
429
430 for (i = 0; i<A_nrows; i++) C_alphaAplusBp[i] = alpha * Ap[i] + Bp[i];
431
432 // Generate C_plusB = this + B
433
434 for (i = 0; i<A_nrows; i++) C_plusBp[i] = Cp[i] + Bp[i];
435
436 // Generate dotvec_AB. Because B(i,j) = 1/A(i,j), dotvec[i] = C.GlobalLength()
437
438 for (i=0; i< A.NumVectors(); i++) dotvec_AB[i] = C.GlobalLength();
439
440 // For the next two results we want to be careful how we do arithmetic
441 // to avoid very large numbers.
442 // We are computing sfinv*(C.GlobalLength()*(C.GlobalLength()+1)/2)
443
444 double result = C.GlobalLength();
445 result *= sfinv;
446 result /= 2.0;
447 result *= (double)(C.GlobalLength()+1);
448
449 // Generate norm1_A. Can use formula for sum of first n integers.
450
451 for (i=0; i< A.NumVectors(); i++)
452 // m1_A[i] = (i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2;
453 norm1_A[i] = result * ((double) (i+1));
454
455 // Generate norm2_sqrtA. Can use formula for sum of first n integers.
456
457 for (i=0; i< A.NumVectors(); i++)
458 // norm2_sqrtA[i] = std::sqrt((double) ((i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2));
459 norm2_sqrtA[i] = std::sqrt(result * ((double) (i+1)));
460
461 // Generate norminf_A, minval_A, maxval_A, meanval_A.
462
463 for (i=0; i< A.NumVectors(); i++)
464 {
465 norminf_A[i] = (double) (i+1);
466 minval_A[i] = (double) (i+1)/ (double) A.GlobalLength();
467 maxval_A[i] = (double) (i+1);
468 meanval_A[i] = norm1_A[i]/((double) (A.GlobalLength()));
469 }
470
471 // Define weights and expected weighted norm
472 normw_A[0] = 1;
473 for (j=0; j<A_nrows; j++) Weightsp[j] = Ap[j];
474
475 delete A_Map;
476 delete B_Map;
477 delete [] A_MyGlobalElements;
478 delete [] B_MyGlobalElements;
479
480 return(0);
481}
482
483//========================================================================
int BuildMatrixTests(Epetra_Vector &C, const char TransA, const char TransB, const double alpha, Epetra_Vector &A, Epetra_Vector &B, const double beta, Epetra_Vector &C_GEMM)
int BuildVectorTests(Epetra_Vector &C, const double alpha, Epetra_Vector &A, Epetra_Vector &sqrtA, Epetra_Vector &B, Epetra_Vector &C_alphaA, Epetra_Vector &C_alphaAplusB, Epetra_Vector &C_plusB, double *const dotvec_AB, double *const norm1_A, double *const norm2_sqrtA, double *const norminf_A, double *const normw_A, Epetra_Vector &Weights, double *const minval_A, double *const maxval_A, double *const meanval_A)
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int IndexBase() const
Index base for this map.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
int GlobalLength() const
Returns the global vector length of vectors in the multi-vector.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ExtractView(double **V) const
Set user-provided address of V.