Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
234int 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
471int 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}
@ View
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 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 MatrixTests(const Epetra_BlockMap &Map, const Epetra_LocalMap &LocalMap, bool verbose)
int BadResidual(bool verbose, double *Residual)
int VectorTests(const Epetra_BlockMap &Map, bool verbose)
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
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)
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
int NormInf(double *Result) const
Compute Inf-norm of each vector in multi-vector.
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
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.
int Norm1(double *Result) const
Compute 1-norm of each vector in multi-vector.
int Update(double ScalarA, const Epetra_MultiVector &A, double ScalarThis)
Update multi-vector values with scaled values of A, this = ScalarThis*this + ScalarA*A.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
int Random()
Set multi-vector values to random numbers.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
static void SetTracebackMode(int TracebackModeValue)
Set the value of the Epetra_Object error traceback report mode.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ExtractCopy(double *V) const
Put vector values into user-provided array.
#define EPETRA_TEST_ERR(a, b)
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)