Epetra Package Browser (Single Doxygen Collection)  Development
Vector/ExecuteTestProblems.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 "Epetra_BLAS.h"
44 #include "ExecuteTestProblems.h"
45 #include "BuildTestProblems.h"
46 #include "Epetra_Comm.h"
47 
48 #ifdef HAVE_EPETRA_TEUCHOS
49 # include "Teuchos_VerboseObject.hpp"
50 #endif
51 
52  int MatrixTests(const Epetra_BlockMap & Map, const Epetra_LocalMap & LocalMap,
53  bool verbose)
54  {
55 
56 #ifdef HAVE_EPETRA_TEUCHOS
57  Teuchos::RCP<Teuchos::FancyOStream>
58  fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
59  std::ostream &out = *fancyOut;
60 #else
61  std::ostream &out = std::cout;
62 #endif
63 
64  int NumVectors = 1;
65  const Epetra_Comm & Comm = Map.Comm();
66  int ierr = 0, i;
67  int IndexBase = 0;
68  double *residual = new double[NumVectors];
69 
70  /* get ID of this processor */
71 
72  // Test GEMM first. 7 cases:
73 
74  // Num
75  // OPERATIONS case Notes
76  // 1) C(local) = A^X(local) * B^X(local) 4 (X=Trans or Not, No Comm needed)
77  // 2) C(local) = A^T(distr) * B (distr) 1 (2D dot product, replicate C)
78  // 3) C(distr) = A (distr) * B^X(local) 2 (2D vector update, no Comm needed)
79 
80  // ==================================================================
81  // Case 1 through 4 (A, B, C all local) Strided and non-strided cases
82  // ==================================================================
83 
84  // Construct Vectors
85 
86  {
87  Epetra_Vector A(LocalMap);
88  Epetra_Vector B(LocalMap);
89  Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
90  Epetra_Vector C(Map2d);
91  Epetra_Vector C_GEMM(Map2d);
92 
93  double **App, **Bpp, **Cpp;
94 
95  Epetra_Vector *Ap, *Bp, *Cp;
96 
97  // For testing non-strided mode, create Vectors that are scattered throughout memory
98 
99  App = new double *[NumVectors];
100  Bpp = new double *[NumVectors];
101  Cpp = new double *[NumVectors];
102  for (i=0; i<NumVectors; i++) App[i] = new double[A.MyLength()+i];
103  for (i=0; i<NumVectors; i++) Bpp[i] = new double[B.MyLength()+i];
104  for (i=0; i<NumVectors; i++) Cpp[i] = new double[C.MyLength()+i];
105 
106  Epetra_Vector A1(View, LocalMap, App[0]);
107  Epetra_Vector B1(View, LocalMap, Bpp[0]);
108  Epetra_Vector C1(View, Map2d, Cpp[0]);
109 
110  for (int strided = 0; strided<2; strided++){
111  // Loop through all trans cases using a variety of values for alpha and beta
112  for (i=0; i<4; i++){
113  int localierr = 0;
114  char transa = 'N'; if (i>1) transa = 'T';
115  char transb = 'N'; if (i%2!=0) transb = 'T';
116  double alpha = (double) i+1;
117  double beta = (double) (i/2);
118  EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
119  localierr = BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
120  if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
121  if (strided)
122  {
123  Ap = &A; Bp = &B; Cp = &C;
124  }
125  else
126  {
127  A.ExtractCopy(App[0]); Ap = &A1;
128  B.ExtractCopy(Bpp[0]); Bp = &B1;
129  C.ExtractCopy(Cpp[0]); Cp = &C1;
130  }
131 
132  localierr = Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
133  if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
134  ierr += Cp->Update(-1.0, C_GEMM, 1.0);
135  ierr += Cp->Norm2(residual);
136 
137  if (verbose && ierr==0)
138  {
139  out << "XXXXX Replicated Local Vector GEMM tests";
140  if (strided)
141  out << " (Strided Multivectors)" << endl;
142  else
143  out << " (Non-Strided Multivectors)" << endl;
144  out << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
145  <<", transb = " << transb;
146  }
147  if (ierr==0 && BadResidual(verbose,residual)) return(-1);
148  }
149  }
150  }
151  }
152  for (i=0; i<NumVectors; i++)
153  {
154  delete [] App[i];
155  delete [] Bpp[i];
156  delete [] Cpp[i];
157  }
158  delete [] App;
159  delete [] Bpp;
160  delete [] Cpp;
161  }
162 
163  // ====================================
164  // Case 5 (A, B distributed C local)
165  // ====================================
166 
167  // Construct Vectors
168  {
169  Epetra_Vector A(Map);
170  Epetra_Vector B(Map);
171  Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
172  Epetra_Vector C(Map2d);
173  Epetra_Vector C_GEMM(Map2d);
174 
175  char transa = 'T';
176  char transb = 'N';
177  double alpha = 2.0;
178  double beta = 1.0;
179  EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
180  ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
181  ierr += C.Multiply(transa, transb, alpha, A, B, beta);
182  ierr += C.Update(-1.0, C_GEMM, 1.0);
183  ierr += C.Norm2(residual);
184 
185  if (verbose && ierr==0)
186  {
187  out << "XXXXX Generalized 2D dot product via GEMM call " << endl;
188  out << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
189  <<", transb = " << transb;
190  }
191  if (BadResidual(verbose,residual)) return(-1);
192 
193 
194  }
195  // ====================================
196  // Case 6-7 (A, C distributed, B local)
197  // ====================================
198 
199  // Construct Vectors
200  {
201  Epetra_Vector A(Map);
202  Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
203  Epetra_Vector B(Map2d);
204  Epetra_Vector C(Map);
205  Epetra_Vector C_GEMM(Map);
206 
207  for (i=0; i<2; i++)
208  {
209  char transa = 'N';
210  char transb = 'N'; if (i>0) transb = 'T';
211  double alpha = 2.0;
212  double beta = 1.1;
213  EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
214  ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
215  ierr += C.Multiply(transa, transb, alpha, A, B, beta);
216  ierr += C.Update(-1.0, C_GEMM, 1.0);
217  ierr += C.Norm2(residual);
218 
219  if (verbose)
220  {
221  out << "XXXXX Generalized 2D vector update via GEMM call " << endl;
222  out << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
223  <<", transb = " << transb;
224  }
225  if (BadResidual(verbose,residual)) return(-1);
226  }
227 
228  delete [] residual;
229 
230  return(ierr);
231  }
232  }
233 
234 int VectorTests(const Epetra_BlockMap & Map, bool verbose)
235 {
236  int NumVectors = 1;
237  int ierr = 0;
238  double *residual = new double[NumVectors];
239 
240 #ifdef HAVE_EPETRA_TEUCHOS
241  Teuchos::RCP<Teuchos::FancyOStream>
242  fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
243  std::ostream &out = *fancyOut;
244 #else
245  std::ostream &out = std::cout;
246 #endif
247 
248  Epetra_BLAS BLAS;
249  /* get number of processors and the name of this processor */
250 
251  // Construct Vectors
252 
253  Epetra_Vector A(Map);
254  Epetra_Vector sqrtA(Map);
255  Epetra_Vector B(Map);
256  Epetra_Vector C(Map);
257  Epetra_Vector C_alphaA(Map);
258  Epetra_Vector C_alphaAplusB(Map);
259  Epetra_Vector C_plusB(Map);
260  Epetra_Vector Weights(Map);
261 
262  // Construct double vectors
263  double *dotvec_AB = new double[NumVectors];
264  double *norm1_A = new double[NumVectors];
265  double *norm2_sqrtA = new double[NumVectors];
266  double *norminf_A = new double[NumVectors];
267  double *normw_A = new double[NumVectors];
268  double *minval_A = new double[NumVectors];
269  double *maxval_A = new double[NumVectors];
270  double *meanval_A = new double[NumVectors];
271 
272  A.SetTracebackMode(1);
273 
274  // Generate data
275 
276 
277  EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
278  double alpha = 2.0;
279  BuildVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
280  C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A,
281  normw_A, Weights, minval_A, maxval_A, meanval_A);
282 
283  int err = 0;
284 
285  // Test range checking for operator[](int)
286 #ifdef HAVE_EPETRA_ARRAY_BOUNDS_CHECK
287  try {
288  if (verbose) out << "XXXXX Testing operator[A.MyLength()] bounds check ";
289  A[A.MyLength()]; // Should throw!
290  // If we get here the test failed!
291  EPETRA_TEST_ERR( 1, ierr );
292  }
293  catch ( const int& err_code ) {
294  if (verbose) out << "\t Checked OK" << endl;
295  EPETRA_TEST_ERR( ( err_code == -99 ? 0 : 1 ), ierr );
296  }
297  try {
298  if (verbose) out << "XXXXX Testing operator[-1] bounds check ";
299  A[-1]; // Should throw!
300  // If we get here the test failed!
301  EPETRA_TEST_ERR( 1, ierr );
302  }
303  catch ( const int& err_code ) {
304  if (verbose) out << "\t Checked OK" << endl;
305  EPETRA_TEST_ERR( ( err_code == -99 ? 0 : 1 ), ierr );
306  }
307 #endif
308 
309  if (verbose) out << "XXXXX Testing alpha * A ";
310  // Test alpha*A
311  Epetra_Vector alphaA(A); // Copy of A
312  EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
313  EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
314  EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
315 
316  if (err) ierr += err;
317  else {
318  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
319  }
320 
321  err = 0;
322  if (verbose) out << "XXXXX Testing C = alpha * A + B ";
323  // Test alpha*A + B
324  Epetra_Vector alphaAplusB(A); // Copy of A
325  EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
326  EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
327  EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
328 
329  if (err) ierr += err;
330  else {
331  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
332  }
333 
334  err = 0;
335  if (verbose) out << "XXXXX Testing C += B ";
336  // Test + B
337  Epetra_Vector plusB(C); // Copy of C
338  EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
339  EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
340  EPETRA_TEST_ERR(plusB.Norm2(residual),err);
341 
342  if (err) ierr += err;
343  else {
344  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
345  }
346 
347  err = 0;
348  if (verbose) out << "XXXXX Testing A.dotProd(B) ";
349  // Test A.dotvec(B)
350  double *dotvec = residual;
351  EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
352  BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
353 
354  if (err) ierr += err;
355  else {
356  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
357  }
358 
359  err = 0;
360  if (verbose) out << "XXXXX Testing norm1_A ";
361  // Test A.norm1()
362  double *norm1 = residual;
363  EPETRA_TEST_ERR(A.Norm1(norm1),err);
364  BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
365 
366  if (err) ierr += err;
367  else {
368  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
369  }
370 
371  err = 0;
372  if (verbose) out << "XXXXX Testing norm2_sqrtA ";
373  // Test sqrtA.norm2()
374  double *norm2 = residual;
375  EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
376  BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
377 
378  if (err) ierr += err;
379  else {
380  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
381  }
382 
383  err = 0;
384  if (verbose) out << "XXXXX Testing norminf_A ";
385  // Test A.norminf()
386  double *norminf = residual;
387  EPETRA_TEST_ERR(A.NormInf(norminf),ierr);
388  BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
389 
390  if (err) ierr += err;
391  else {
392  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
393  }
394 
395  err = 0;
396  if (verbose) out << "XXXXX Testing normw_A ";
397  // Test A.NormWeighted()
398  double *normw = residual;
399  EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
400  BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
401 
402  if (err) ierr += err;
403  else {
404  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
405  }
406 
407  err = 0;
408  if (verbose) out << "XXXXX Testing minval_A ";
409  // Test A.MinValue()
410  double *minval = residual;
411  EPETRA_TEST_ERR(A.MinValue(minval),err);
412  BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
413 
414  if (err) ierr += err;
415  else {
416  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
417  }
418 
419  err = 0;
420  if (verbose) out << "XXXXX Testing maxval_A ";
421  // Test A.MaxValue()
422  double *maxval = residual;
423  EPETRA_TEST_ERR(A.MaxValue(maxval),err);
424  BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
425 
426  if (err) ierr += err;
427  else {
428  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
429  }
430 
431  err = 0;
432  if (verbose) out << "XXXXX Testing meanval_A ";
433  // Test A.MeanValue()
434  double *meanval = residual;
435  EPETRA_TEST_ERR(A.MeanValue(meanval),err);
436  BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
437 
438  if (err) ierr += err;
439  else {
440  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
441  }
442 
443  err = 0;
444  if (verbose) out << "XXXXX Testing abs_A ";
445  // Test A.Abs()
446  Epetra_Vector Abs_A = A;
447  EPETRA_TEST_ERR(Abs_A.Abs(A),err);
448  EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
449  EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
450 
451  if (err) ierr += err;
452  else {
453  EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
454  }
455 
456  // Delete everything
457 
458  delete [] dotvec_AB;
459  delete [] norm1_A;
460  delete [] norm2_sqrtA;
461  delete [] norminf_A;
462  delete [] normw_A;
463  delete [] minval_A;
464  delete [] maxval_A;
465  delete [] meanval_A;
466  delete [] residual;
467 
468  return(ierr);
469 }
470 
471 int BadResidual(bool verbose, double * Residual)
472 {
473 
474 #ifdef HAVE_EPETRA_TEUCHOS
475  Teuchos::RCP<Teuchos::FancyOStream>
476  fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
477  std::ostream &out = *fancyOut;
478 #else
479  std::ostream &out = std::cout;
480 #endif
481 
482  double threshold = 5.0E-6;
483  int ierr = 0;
484  if (Residual[0]>threshold) {
485  ierr = 1;// Output will be more useful after returning from method
486  if (verbose) out << endl << " Residual = " << Residual[0];
487  }
488  if (verbose)
489  if (ierr==0) out << "\t Checked OK" << endl;
490 
491  return(ierr);
492 }
MatrixTests
int MatrixTests(const Epetra_BlockMap &Map, const Epetra_LocalMap &LocalMap, bool verbose)
Definition: Vector/ExecuteTestProblems.cpp:52
B1
Epetra_BlockMap::Comm
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Definition: Epetra_BlockMap.h:770
Teuchos::VerboseObjectBase::getDefaultOStream
static RCP< FancyOStream > getDefaultOStream()
VectorTests
int VectorTests(const Epetra_BlockMap &Map, bool verbose)
Definition: Vector/ExecuteTestProblems.cpp:234
View
Definition: Epetra_DataAccess.h:57
Epetra_BLAS::AXPY
void AXPY(const int N, const float ALPHA, const float *X, float *Y, const int INCX=1, const int INCY=1) const
Epetra_BLAS vector update function (SAXPY)
Definition: Epetra_BLAS.cpp:117
Epetra_BLAS.h
Epetra_MultiVector::Abs
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
Definition: Epetra_MultiVector.cpp:1158
EPETRA_TEST_ERR
#define EPETRA_TEST_ERR(a, b)
Definition: epetra_test_err.h:55
Epetra_Comm
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
Teuchos::RCP
BuildMatrixTests
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)
Definition: MultiVector/BuildTestProblems.cpp:45
Epetra_MultiVector::Scale
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
Definition: Epetra_MultiVector.cpp:1229
ExecuteTestProblems.h
BadResidual
int BadResidual(bool verbose, double *Residual)
Definition: Vector/ExecuteTestProblems.cpp:471
Epetra_MultiVector::Norm2
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
Definition: Epetra_MultiVector.cpp:1547
Epetra_BLAS
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
Epetra_MultiVector::Multiply
int Multiply(char TransA, char TransB, double ScalarAB, const Epetra_MultiVector &A, const Epetra_MultiVector &B, double ScalarThis)
Matrix-Matrix multiplication, this = ScalarThis*this + ScalarAB*A*B.
Definition: Epetra_MultiVector.cpp:1936
Epetra_BlockMap
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
Definition: Epetra_BlockMap.h:194
Epetra_Vector
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
Definition: Epetra_Vector.h:142
Epetra_Comm.h
C
BuildVectorTests
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)
Definition: Vector/BuildTestProblems.cpp:274
A
Residual
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
Definition: test/SerialDense/cxx_main.cpp:649
Epetra_LocalMap
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
Definition: Epetra_LocalMap.h:89
BuildTestProblems.h
Teuchos_VerboseObject.hpp
B
Epetra_MultiVector::Update
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
Definition: Epetra_MultiVector.cpp:1276
Generated by   doxygen 1.8.16