Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MultiVector/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_MultiVector& C_GEMM ) {
52 
53  // For given values of TransA, TransB, alpha and beta, a (possibly
54  // zero) filled Epetra_MultiVector C, and allocated
55  // Epetra_MultiVectors A, B and C_GEMM this routine will generate values for
56  // Epetra_MultiVectors 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 MultiVectors
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, &A_Stride);
89  B.ExtractView(&Bp, &B_Stride);
90  C.ExtractView(&Cp, &C_Stride);
91  C_GEMM.ExtractView(&C_GEMMp, &C_GEMM_Stride);
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  }
274 int BuildMultiVectorTests (Epetra_MultiVector & C, const double alpha,
276  Epetra_MultiVector& sqrtA,
278  Epetra_MultiVector& C_alphaA,
279  Epetra_MultiVector& C_alphaAplusB,
280  Epetra_MultiVector& 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_MultiVector& 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_MultiVector (the this object), allocated double * arguments dotvec_AB,
293  // norm1_A, and norm2_A, and allocated Epetra_MultiVectors A, sqrtA,
294  // B, C_alpha, C_alphaAplusB and C_plusB, this method will generate values for
295  // Epetra_MultiVectors 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, fj; // 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  {
396  for (j = 0; j <A_ncols ; j++)
397  {
398  for (i = 0; i<A_nrows; i++)
399  {
400  fi = i+1; // Get float version of i and j, offset by 1.
401  fj = j+1;
402  Ap[j][i] = (fi*sfinv)*fj;
403  sqrtAp[j][i] = std::sqrt(Ap[j][i]);
404  }
405  }
406  }
407  else
408  {
409  for (j = 0; j <A_ncols ; j++)
410  {
411  for (i = 0; i<A_nrows; i++)
412  {
413  fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
414  fj = j+1;
415  Ap[j][i] = (fi*sfinv)*fj;
416  sqrtAp[j][i] = std::sqrt(Ap[j][i]);
417  }
418  }
419  }
420 
421  // Define B depending on TransB and B_is_local
422 
423  if (B_is_local)
424  {
425  for (j = 0; j <B_ncols ; j++)
426  {
427  for (i = 0; i<B_nrows; i++)
428  {
429  fi = i+1; // Get float version of i and j, offset by 1.
430  fj = j+1;
431  Bp[j][i] = 1.0/((fi*sfinv)*fj);
432  }
433  }
434  }
435  else
436  {
437  for (j = 0; j <B_ncols ; j++)
438  {
439  for (i = 0; i<B_nrows; i++)
440  {
441  fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
442  fj = j+1;
443  Bp[j][i] = 1.0/((fi*sfinv)*fj);
444  }
445  }
446  }
447 
448  // Generate C_alphaA = alpha * A
449 
450  for (j = 0; j <A_ncols ; j++)
451  for (i = 0; i<A_nrows; i++)
452  C_alphaAp[j][i] = alpha * Ap[j][i];
453 
454  // Generate C_alphaA = alpha * A + B
455 
456  for (j = 0; j <A_ncols ; j++)
457  for (i = 0; i<A_nrows; i++)
458  C_alphaAplusBp[j][i] = alpha * Ap[j][i] + Bp[j][i];
459 
460  // Generate C_plusB = this + B
461 
462  for (j = 0; j <A_ncols ; j++)
463  for (i = 0; i<A_nrows; i++)
464  C_plusBp[j][i] = Cp[j][i] + Bp[j][i];
465 
466  // Generate dotvec_AB. Because B(i,j) = 1/A(i,j), dotvec[i] = C.GlobalLength()
467 
468  for (i=0; i< A.NumVectors(); i++) dotvec_AB[i] = C.GlobalLength();
469 
470  // For the next two results we want to be careful how we do arithmetic
471  // to avoid very large numbers.
472  // We are computing sfinv*(C.GlobalLength()*(C.GlobalLength()+1)/2)
473 
474  double result = C.GlobalLength();
475  result *= sfinv;
476  result /= 2.0;
477  result *= (double)(C.GlobalLength()+1);
478 
479  // Generate norm1_A. Can use formula for sum of first n integers.
480 
481  for (i=0; i< A.NumVectors(); i++)
482  // m1_A[i] = (i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2;
483  norm1_A[i] = result * ((double) (i+1));
484 
485  // Generate norm2_sqrtA. Can use formula for sum of first n integers.
486 
487  for (i=0; i< A.NumVectors(); i++)
488  // norm2_sqrtA[i] = std::sqrt((double) ((i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2));
489  norm2_sqrtA[i] = std::sqrt(result * ((double) (i+1)));
490 
491  // Generate norminf_A, minval_A, maxval_A, meanval_A.
492 
493  for (i=0; i< A.NumVectors(); i++)
494  {
495  norminf_A[i] = (double) (i+1);
496  minval_A[i] = (double) (i+1)/ (double) A.GlobalLength();
497  maxval_A[i] = (double) (i+1);
498  meanval_A[i] = norm1_A[i]/((double) (A.GlobalLength()));
499  }
500 
501  // Define weights and expected weighted norm
502  for (i=0; i< A.NumVectors(); i++)
503  {
504  double ip1 = (double) i+1;
505  normw_A[i] = ip1;
506  for (j=0; j<A_nrows; j++) Weightsp[i][j] = Ap[i][j]/ip1;
507  }
508 
509 
510  delete A_Map;
511  delete B_Map;
512  delete [] A_MyGlobalElements;
513  delete [] B_MyGlobalElements;
514 
515  return(0);
516 }
517 
518 //========================================================================
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int GlobalLength() const
Returns the global vector length of vectors in the 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 MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int BuildMatrixTests(Epetra_MultiVector &C, const char TransA, const char TransB, const double alpha, Epetra_MultiVector &A, Epetra_MultiVector &B, const double beta, Epetra_MultiVector &C_GEMM)
int IndexBase() const
Index base for this map.
int ExtractView(double **A, int *MyLDA) const
Set user-provided addresses of A and MyLDA.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int BuildMultiVectorTests(Epetra_MultiVector &C, const double alpha, Epetra_MultiVector &A, Epetra_MultiVector &sqrtA, Epetra_MultiVector &B, Epetra_MultiVector &C_alphaA, Epetra_MultiVector &C_alphaAplusB, Epetra_MultiVector &C_plusB, double *const dotvec_AB, double *const norm1_A, double *const norm2_sqrtA, double *const norminf_A, double *const normw_A, Epetra_MultiVector &Weights, double *const minval_A, double *const maxval_A, double *const meanval_A)
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.