Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
MultiVector/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#include "Epetra_Vector.h"
48#include "Epetra_IntVector.h"
49#include "Epetra_Import.h"
50
51 int MatrixTests(const Epetra_BlockMap & Map, const Epetra_LocalMap & LocalMap, int NumVectors,
52 bool verbose)
53 {
54 const Epetra_Comm & Comm = Map.Comm();
55 int ierr = 0, i;
56 int IndexBase = 0;
57 double *residual = new double[NumVectors];
58
59 /* get ID of this processor */
60
61
62 // Test GEMM first. 7 cases:
63
64 // Num
65 // OPERATIONS case Notes
66 // 1) C(local) = A^X(local) * B^X(local) 4 (X=Trans or Not, No Comm needed)
67 // 2) C(local) = A^T(distr) * B (distr) 1 (2D dot product, replicate C)
68 // 3) C(distr) = A (distr) * B^X(local) 2 (2D vector update, no Comm needed)
69
70 // ==================================================================
71 // Case 1 through 4 (A, B, C all local) Strided and non-strided cases
72 // ==================================================================
73
74 // Construct MultiVectors
75
76 {
77 Epetra_MultiVector A(LocalMap, NumVectors);
78 Epetra_MultiVector B(LocalMap, NumVectors);
79 Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
80 Epetra_MultiVector C(Map2d, NumVectors);
81 Epetra_MultiVector C_GEMM(Map2d, NumVectors);
82
83 double **App, **Bpp, **Cpp;
84
85 Epetra_MultiVector *Ap, *Bp, *Cp;
86
87 // For testing non-strided mode, create MultiVectors that are scattered throughout memory
88
89 App = new double *[NumVectors];
90 Bpp = new double *[NumVectors];
91 Cpp = new double *[NumVectors];
92 for (i=0; i<NumVectors; i++) App[i] = new double[A.MyLength()+i];
93 for (i=0; i<NumVectors; i++) Bpp[i] = new double[B.MyLength()+i];
94 for (i=0; i<NumVectors; i++) Cpp[i] = new double[C.MyLength()+i];
95
96 Epetra_MultiVector A1(View, LocalMap, App, NumVectors);
97 Epetra_MultiVector B1(View, LocalMap, Bpp, NumVectors);
98 Epetra_MultiVector C1(View, Map2d, Cpp, NumVectors);
99
100 for (int strided = 0; strided<2; strided++) {
101
102 // Loop through all trans cases using a variety of values for alpha and beta
103 for (i=0; i<4; i++) {
104 char transa = 'N'; if (i>1) transa = 'T';
105 char transb = 'N'; if (i%2!=0) transb = 'T';
106 double alpha = (double) i+1;
107 double beta = (double) (i/2);
108 EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
109 int localierr = BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
110 if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
111 if (strided)
112 {
113 Ap = &A; Bp = &B; Cp = &C;
114 }
115 else
116 {
117 A.ExtractCopy(App); Ap = &A1;
118 B.ExtractCopy(Bpp); Bp = &B1;
119 C.ExtractCopy(Cpp); Cp = &C1;
120 }
121
122 localierr = Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
123 if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
124 ierr += Cp->Update(-1.0, C_GEMM, 1.0);
125 ierr += Cp->Norm2(residual);
126
127 if (verbose)
128 {
129 cout << "XXXXX Replicated Local MultiVector GEMM tests";
130 if (strided)
131 cout << " (Strided Multivectors)" << endl;
132 else
133 cout << " (Non-Strided Multivectors)" << endl;
134 cout << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
135 <<", transb = " << transb;
136 }
137 if (BadResidual(verbose,residual, NumVectors)) return(-1);
138 }
139 }
140 }
141
142 }
143 for (i=0; i<NumVectors; i++)
144 {
145 delete [] App[i];
146 delete [] Bpp[i];
147 delete [] Cpp[i];
148 }
149 delete [] App;
150 delete [] Bpp;
151 delete [] Cpp;
152 }
153
154 // ====================================
155 // Case 5 (A, B distributed C local)
156 // ====================================
157
158 // Construct MultiVectors
159 {
160 Epetra_MultiVector A(Map, NumVectors);
161 Epetra_MultiVector B(Map, NumVectors);
162 Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
163 Epetra_MultiVector C(Map2d, NumVectors);
164 Epetra_MultiVector C_GEMM(Map2d, NumVectors);
165
166 char transa = 'T';
167 char transb = 'N';
168 double alpha = 2.0;
169 double beta = 1.0;
170 EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
171 ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
172 int localierr = C.Multiply(transa, transb, alpha, A, B, beta);
173 if (localierr!=-2) { // -2 means the shapes didn't match
174 ierr += C.Update(-1.0, C_GEMM, 1.0);
175 ierr += C.Norm2(residual);
176
177 if (verbose)
178 {
179 cout << "XXXXX Generalized 2D dot product via GEMM call " << endl;
180 cout << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
181 <<", transb = " << transb;
182 }
183 if (BadResidual(verbose,residual, NumVectors)) return(-1);
184 }
185
186 }
187 // ====================================
188 // Case 6-7 (A, C distributed, B local)
189 // ====================================
190
191 // Construct MultiVectors
192 {
193 Epetra_MultiVector A(Map, NumVectors);
194 Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
195 Epetra_MultiVector B(Map2d, NumVectors);
196 Epetra_MultiVector C(Map, NumVectors);
197 Epetra_MultiVector C_GEMM(Map, NumVectors);
198
199 for (i=0; i<2; i++)
200 {
201 char transa = 'N';
202 char transb = 'N'; if (i>0) transb = 'T';
203 double alpha = 2.0;
204 double beta = 1.1;
205 EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers
206 ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
207 ierr += C.Multiply(transa, transb, alpha, A, B, beta);
208 ierr += C.Update(-1.0, C_GEMM, 1.0);
209 ierr += C.Norm2(residual);
210
211 if (verbose)
212 {
213 cout << "XXXXX Generalized 2D vector update via GEMM call " << endl;
214 cout << " alpha = " << alpha << ", beta = " << beta <<", transa = "<<transa
215 <<", transb = " << transb;
216 }
217 if (BadResidual(verbose,residual, NumVectors)) return(-1);
218 }
219
220
221 }
222 // ====================================
223 // LocalMap Tests
224 // ====================================
225
226 // Construct MultiVectors
227 {
228
229 int localLength = 10;
230 double *localMinValue = new double[localLength];
231 double *localMaxValue = new double[localLength];
232 double *localNorm1 = new double[localLength];
233 double *localDot = new double[localLength];
234 double *localNorm2 = new double[localLength];
235 double *localMeanValue = new double[localLength];
236 Epetra_LocalMap MapSmall(localLength, IndexBase, Comm);
237 Epetra_MultiVector A(MapSmall, NumVectors);
238
239 double doubleLocalLength = (double) localLength;
240 for (int j=0; j< NumVectors; j++) {
241 for (i=0; i< localLength-1; i++) A[j][i] = (double) (i+1);
242 A[j][localLength-1] = (double) (localLength+j); // Only the last value differs across multivectors
243 localMinValue[j] = A[j][0]; // Increasing values
244 localMaxValue[j] = A[j][localLength-1];
245 localNorm1[j] = (doubleLocalLength-1.0)*(doubleLocalLength)/2.0+A[j][localLength-1];
246 localDot[j] = (doubleLocalLength-1.0)*(doubleLocalLength)*(2.0*(doubleLocalLength-1.0)+1.0)/6.0+A[j][localLength-1]*A[j][localLength-1];
247 localNorm2[j] = std::sqrt(localDot[j]);
248 localMeanValue[j] = localNorm1[j]/doubleLocalLength;
249 }
250 ierr += A.MinValue(residual);
251 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMinValue[j]);
252 if (verbose) cout << "XXXXX MinValue" << endl;
253 if (BadResidual(verbose,residual, NumVectors)) return(-1);
254
255 ierr += A.MaxValue(residual);
256 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMaxValue[j]);
257 if (verbose) cout << "XXXXX MaxValue" << endl;
258 if (BadResidual(verbose,residual, NumVectors)) return(-1);
259
260 ierr += A.Norm1(residual);
261 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm1[j]);
262 if (verbose) cout << "XXXXX Norm1" << endl;
263 if (BadResidual(verbose,residual, NumVectors)) return(-1);
264
265 ierr += A.Dot(A,residual);
266 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localDot[j]);
267 if (verbose) cout << "XXXXX Dot" << endl;
268 if (BadResidual(verbose,residual, NumVectors)) return(-1);
269
270 ierr += A.Norm2(residual);
271 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm2[j]);
272 if (verbose) cout << "XXXXX Norm2" << endl;
273 if (BadResidual(verbose,residual, NumVectors)) return(-1);
274
275 ierr += A.MeanValue(residual);
276 for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMeanValue[j]);
277 if (verbose) cout << "XXXXX MeanValue" << endl;
278 if (BadResidual(verbose,residual, NumVectors)) return(-1);
279
280 delete [] localMinValue;
281 delete [] localMaxValue;
282 delete [] localNorm1;
283 delete [] localDot;
284 delete [] localNorm2;
285 delete [] localMeanValue;
286
287 }
288
289 delete [] residual;
290
291 return(ierr);
292 }
293
294int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
295{
296 const Epetra_Comm & Comm = Map.Comm();
297 int ierr = 0, i;
298 double *residual = new double[NumVectors];
299
300 Epetra_BLAS BLAS;
301 /* get number of processors and the name of this processor */
302
303 // int NumProc = Comm.getNumProc();
304 int MyPID = Comm.MyPID();
305
306 // Construct MultiVectors
307
308 Epetra_MultiVector A(Map, NumVectors);
309 Epetra_MultiVector sqrtA(Map, NumVectors);
310 Epetra_MultiVector B(Map, NumVectors);
311 Epetra_MultiVector C(Map, NumVectors);
312 Epetra_MultiVector C_alphaA(Map, NumVectors);
313 Epetra_MultiVector C_alphaAplusB(Map, NumVectors);
314 Epetra_MultiVector C_plusB(Map, NumVectors);
315 Epetra_MultiVector Weights(Map, NumVectors);
316
317 // Construct double vectors
318 double *dotvec_AB = new double[NumVectors];
319 double *norm1_A = new double[NumVectors];
320 double *norm2_sqrtA = new double[NumVectors];
321 double *norminf_A = new double[NumVectors];
322 double *normw_A = new double[NumVectors];
323 double *minval_A = new double[NumVectors];
324 double *maxval_A = new double[NumVectors];
325 double *meanval_A = new double[NumVectors];
326
327 // Generate data
328
329
330 EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
331 double alpha = 2.0;
332 BuildMultiVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
333 C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A,
334 normw_A, Weights, minval_A, maxval_A, meanval_A);
335
336 int err = 0;
337 if (verbose) cout << "XXXXX Testing alpha * A ";
338 // Test alpha*A
339 Epetra_MultiVector alphaA(A); // Copy of A
340 EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
341 EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
342 EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
343
344 if (err) ierr += err;
345 else {
346 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
347 }
348
349 err = 0;
350 if (verbose) cout << "XXXXX Testing C = alpha * A + B ";
351 // Test alpha*A + B
352 Epetra_MultiVector alphaAplusB(A); // Copy of A
353 EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
354 EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
355 EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
356
357 if (err) ierr += err;
358 else {
359 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
360 }
361
362 err = 0;
363 if (verbose) cout << "XXXXX Testing C += B ";
364 // Test + B
365 Epetra_MultiVector plusB(C); // Copy of C
366 EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
367 EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
368 EPETRA_TEST_ERR(plusB.Norm2(residual),err);
369
370 if (err) ierr += err;
371 else {
372 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
373 }
374
375 err = 0;
376 if (verbose) cout << "XXXXX Testing A.dotProd(B) ";
377 // Test A.dotvec(B)
378 double *dotvec = residual;
379 EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
380 BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
381
382 if (err) ierr += err;
383 else {
384 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
385 }
386
387 err = 0;
388 if (verbose) cout << "XXXXX Testing norm1_A ";
389 // Test A.norm1()
390 double *norm1 = residual;
391 EPETRA_TEST_ERR(A.Norm1(norm1),err);
392 BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
393
394 if (err) ierr += err;
395 else {
396 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
397 }
398
399 err = 0;
400 if (verbose) cout << "XXXXX Testing norm2_sqrtA ";
401 // Test sqrtA.norm2()
402 double *norm2 = residual;
403 EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
404 BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
405
406 if (err) ierr += err;
407 else {
408 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
409 }
410
411 err = 0;
412 if (verbose) cout << "XXXXX Testing norminf_A ";
413 // Test A.norminf()
414 double *norminf = residual;
415 EPETRA_TEST_ERR(A.NormInf(norminf),err);
416 BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
417
418 if (err) ierr += err;
419 else {
420 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
421 }
422
423 err = 0;
424 if (verbose) cout << "XXXXX Testing normw_A ";
425 // Test A.NormWeighted()
426 double *normw = residual;
427 EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
428 BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
429
430 if (err) ierr += err;
431 else {
432 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
433 }
434
435 err = 0;
436 if (verbose) cout << "XXXXX Testing minval_A ";
437 // Test A.MinValue()
438 double *minval = residual;
439 EPETRA_TEST_ERR(A.MinValue(minval),err);
440 BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
441
442 if (err) ierr += err;
443 else {
444 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
445 }
446
447 err = 0;
448 if (verbose) cout << "XXXXX Testing maxval_A ";
449 // Test A.MaxValue()
450 double *maxval = residual;
451 EPETRA_TEST_ERR(A.MaxValue(maxval),err);
452 BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
453
454 if (err) ierr += err;
455 else {
456 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
457 }
458
459 err = 0;
460 if (verbose) cout << "XXXXX Testing meanval_A ";
461 // Test A.MeanValue()
462 double *meanval = residual;
463 EPETRA_TEST_ERR(A.MeanValue(meanval),err);
464 BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
465
466 if (err) ierr += err;
467 else {
468 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
469 }
470
471 err = 0;
472 if (verbose) cout << "XXXXX Testing abs_A ";
473 // Test A.Abs()
474 Epetra_MultiVector Abs_A = A;
475 EPETRA_TEST_ERR(Abs_A.Abs(A),err);
476 EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
477 EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
478
479 if (err) ierr += err;
480 else {
481 EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
482 }
483
484 err = 0;
485 if (verbose) cout << "XXXXX Testing random_A (Test1) ";
486 // Test A.Random()
487 Epetra_MultiVector Rand1_A(A);
488 Epetra_MultiVector Rand2_A(A);
489 EPETRA_TEST_ERR(Rand1_A.Random(),err);
490 EPETRA_TEST_ERR(Rand2_A.Random(),err);
491 // Rand2_A = Rand1_A - Rand2_A (should be nonzero since Random() should give different vectors > 0)
492 EPETRA_TEST_ERR(Rand2_A.Update(1.0, Rand1_A, -1.0),err);
493 EPETRA_TEST_ERR(Rand2_A.Norm2(residual),err);
494
495 if (err) ierr += err;
496 else {
497 EPETRA_TEST_ERR(BadResidual1(verbose,residual, NumVectors),ierr);
498 }
499
500 err = 0;
501 if (verbose) cout << "XXXXX Testing random_A (Test2) ";
502
503 // Next test that each column of the multivector is different from all other columns by testing the first value
504 // of each vector against the first value of every other vector.
505 int randvalsdiffer = 1; // Assume they all differ
506 for (i=0; i< NumVectors; i++)
507 for (int j=i+1; j<NumVectors; j++)
508 if (Rand1_A[i][0]==Rand1_A[j][0]) randvalsdiffer = 0; // make false if equal
509 int allrandvals = 0;
510 Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
511
512 EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
513 int locerr = err;
514 Comm.MinAll(&locerr, &err, 1);
515
516 if (verbose) {
517 if (err==0) {
518 cout << "\t Checked OK" << endl;
519 } else {
520 cout << "\t Checked Failed" << endl;
521 }
522 }
523 err = 0;
524 if (verbose) cout << "XXXXX Testing random_A (Test3) ";
525
526 // Next test that the first element on each processor of the first column of Rand1_A is different from all others
527 // First we will gather them all to PE 0
528
529
530 Epetra_Map RandstartsMap(-1, 1, 0, Comm); // This Map has a single element on each PE
531 int itmp = 0;
532 int nproc = Comm.NumProc();
533 if (MyPID==0) itmp = nproc;
534 Epetra_Map AllrandstartsMap(nproc, itmp, 0, Comm); // Map has NumProc elements on PE 0, none elsewhere
535 Epetra_MultiVector Randstarts(RandstartsMap, NumVectors);
536 Epetra_MultiVector Allrandstarts(AllrandstartsMap, NumVectors);
537 for (i=0; i< NumVectors; i++) Randstarts[i][0] = Rand1_A[i][0]; // Load first value of local multivector
538
539 Epetra_Import Randimporter(AllrandstartsMap,RandstartsMap);
540 EPETRA_TEST_ERR(Allrandstarts.Import(Randstarts,Randimporter,Insert),err);
541 // cout << "Randstarts = " << Randstarts << endl << "Allrandstarts = " << Allrandstarts << endl;
542 // Allrandstarts now contains the first values for each local section of Rand1_A.
543 // Next test that this is true.
544 randvalsdiffer = 1; // Assume they all differ
545 if (MyPID==0) {
546 for (i=0; i< NumVectors; i++)
547 for (int irand=0; irand<nproc; irand++)
548 for (int jrand=irand+1; jrand<nproc; jrand++)
549 if (Allrandstarts[i][irand]==Allrandstarts[i][jrand]) randvalsdiffer = 0; // make false if equal
550 }
551 allrandvals = 0;
552 Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
553
554 EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
555 locerr = err;
556 Comm.MinAll(&locerr, &err, 1);
557 if (verbose) {
558 if (err==0) {
559 cout << "\t Checked OK" << endl;
560 } else {
561 cout << "\t Checked Failed" << endl;
562 }
563 }
564
565 // Delete everything
566
567 delete [] dotvec_AB;
568 delete [] norm1_A;
569 delete [] norm2_sqrtA;
570 delete [] norminf_A;
571 delete [] normw_A;
572 delete [] minval_A;
573 delete [] maxval_A;
574 delete [] meanval_A;
575 delete [] residual;
576
577 //*******************************************************************
578 // Post-construction modification tests
579 //*******************************************************************
580
581 if (verbose) cout << "\n\nXXXXX Testing Post-construction modification of a multivector"
582 <<endl<<endl;
583
584 err = 0;
585
586 Epetra_MultiVector X(Map, NumVectors);
587 X.Random();
588
589 // Pick middle range values for GID, LID and Vector Index
590 int testGID = Map.NumGlobalElements()/2;
591 int testVecIndex = NumVectors/2;
592
593 int GIDSize = 1;
594 int LIDOfGID = 0;
595 int FirstEntryOfGID = 0;
596
597 if (Map.MyGID(testGID)) {
598 LIDOfGID = Map.LID(testGID);
599 GIDSize = Map.ElementSize(LIDOfGID);
600 FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
601 }
602
603 // ========================================================================
604 // Test int ReplaceGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
605 // ========================================================================
606
607 double newGIDValue = 4.0;
608 locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
609
610 if (Map.MyGID(testGID)) {
611 if (X[testVecIndex][FirstEntryOfGID]!=newGIDValue) err++;
612 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
613 << X[testVecIndex][FirstEntryOfGID]
614 << " should = " << newGIDValue << endl;
615 }
616 else
617 if (locerr!=1) err++; // Test for GID out of range error (=1)
618
619 // ========================================================================
620 // Test int ReplaceGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
621 // ========================================================================
622 newGIDValue = 8.0;
623 locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
624
625 if (Map.MyGID(testGID)) {
626 if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=newGIDValue) err++;
627 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
628 << X[testVecIndex][FirstEntryOfGID+GIDSize-1]
629 << " should = " << newGIDValue << endl;
630 }
631 else
632 if (locerr!=1) err++; // Test for GID out of range error (=1)
633
634 // ========================================================================
635 // Test int SumIntoGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
636 // ========================================================================
637
638 newGIDValue = 1.0;
639 locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
640 locerr = X.SumIntoGlobalValue(testGID, testVecIndex, newGIDValue);
641 if (Map.MyGID(testGID)) {
642 if (X[testVecIndex][FirstEntryOfGID]!=(newGIDValue+newGIDValue)) err++;
643 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
644 << X[testVecIndex][FirstEntryOfGID]
645 << " should = " << newGIDValue << endl;
646 }
647 else
648 if (locerr!=1) err++; // Test for GID out of range error (=1)
649
650 // ========================================================================
651 // Test int SumIntoGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
652 // ========================================================================
653
654 newGIDValue = 1.0;
655 locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
656 locerr = X.SumIntoGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
657
658 if (Map.MyGID(testGID)) {
659 if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=(newGIDValue+newGIDValue)) err++;
660 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
661 << X[testVecIndex][FirstEntryOfGID+GIDSize-1]
662 << " should = " << newGIDValue << endl;
663 }
664 else
665 if (locerr!=1) err++; // Test for GID out of range error (=1)
666
667 // ========================================================================
668 // Test Local "My" versions of same routine (less complicated)
669 // ========================================================================
670
671 // Pick middle range values for LID
672 int testLID = Map.NumMyElements()/2;
673
674 int LIDSize = Map.ElementSize(testLID);
675 int FirstEntryOfLID = Map.FirstPointInElement(testLID);
676
677
678 double newLIDValue = 4.0;
679 locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
680
681 if (X[testVecIndex][FirstEntryOfLID]!=newLIDValue) err++;
682 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
683 << X[testVecIndex][FirstEntryOfLID]
684 << " should = " << newLIDValue << endl;
685
686 newLIDValue = 8.0;
687 locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
688 if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=newLIDValue) err++;
689 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
690 << X[testVecIndex][FirstEntryOfLID+LIDSize-1]
691 << " should = " << newLIDValue << endl;
692 newLIDValue = 1.0;
693 locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
694 locerr = X.SumIntoMyValue(testLID, testVecIndex, newLIDValue);
695 if (X[testVecIndex][FirstEntryOfLID]!=(newLIDValue+newLIDValue)) err++;
696 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
697 << X[testVecIndex][FirstEntryOfLID]
698 << " should = " << newLIDValue << endl;
699 newLIDValue = 2.0;
700 locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
701 locerr = X.SumIntoMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
702 if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
703 << X[testVecIndex][FirstEntryOfLID+LIDSize-1]
704 << " should = " << newLIDValue << endl;
705 if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=(newLIDValue+newLIDValue)) err++;
706
707 ierr += err;
708
709 // ========================================================================
710 // Test Post-construction modification of an Epetra_Vector using a vector
711 // our multivector X
712 // ========================================================================
713
714 if (verbose) cout << "\n\nXXXXX Testing Post-construction modification of a vector"
715 << endl << endl;
716
717 Epetra_Vector * x = X(testVecIndex);
718
719 int NumEntries = 2;
720 double * VecValues = new double[NumEntries];
721 int * VecGIDs = new int[NumEntries];
722 VecGIDs[0] = testGID;
723 VecGIDs[1] = testGID+1; // Some pathological chance that these GIDs are not valid
724
725 // ========================================================================
726 // Test int ReplaceGlobalValues (int NumEntries, double *Values, int *Indices)
727 // ========================================================================
728
729 VecValues[0] = 2.0; VecValues[1] = 4.0;
730 locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
731
732 for (i=0; i<NumEntries; i++) {
733 testGID = VecGIDs[i];
734 if (Map.MyGID(testGID)) {
735 LIDOfGID = Map.LID(testGID);
736 GIDSize = EPETRA_MIN(GIDSize,Map.ElementSize(LIDOfGID)); // Need this value below
737 FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
738 if ((*x)[FirstEntryOfGID]!=VecValues[i]) err++;
739 if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
740 << (*x)[FirstEntryOfGID]
741 << " should = " << VecValues[i] << endl;
742 }
743 else
744 if (locerr!=1) err++; // Test for GID out of range error (=1)
745 }
746
747
748 // ========================================================================
749 // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
750 // ========================================================================
751
752 VecValues[0] = 4.0; VecValues[1] = 8.0;
753 locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
754
755 for (i=0; i<NumEntries; i++) {
756 testGID = VecGIDs[i];
757 if (Map.MyGID(testGID)) {
758 LIDOfGID = Map.LID(testGID);
759 FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
760 if ((*x)[FirstEntryOfGID+GIDSize-1]!=VecValues[i]) err++;
761 if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
762 << (*x)[FirstEntryOfGID+GIDSize-1]
763 << " should = " << VecValues[i] << endl;
764 }
765 else
766 if (locerr!=1) err++; // Test for GID out of range error (=1)
767 }
768
769 // ========================================================================
770 // Test int SumIntoGlobalValues (int NumEntries, double *Values, int *Indices)
771 // ========================================================================
772
773 VecValues[0] = 1.0; VecValues[1] = 2.0;
774 locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
775 locerr = x->SumIntoGlobalValues(NumEntries, VecValues, VecGIDs);
776
777 for (i=0; i<NumEntries; i++) {
778 testGID = VecGIDs[i];
779 if (Map.MyGID(testGID)) {
780 LIDOfGID = Map.LID(testGID);
781 FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
782 if ((*x)[FirstEntryOfGID]!=(VecValues[i]+VecValues[i])) err++;
783 if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
784 << (*x)[FirstEntryOfGID]
785 << " should = " << (VecValues[i]+VecValues[i]) << endl;
786 }
787 else
788 if (locerr!=1) err++; // Test for GID out of range error (=1)
789 }
790 // ========================================================================
791 // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
792 // ========================================================================
793
794 VecValues[0] = 1.0; VecValues[1] = 2.0;
795 locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
796 locerr = x->SumIntoGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
797
798 for (i=0; i<NumEntries; i++) {
799 testGID = VecGIDs[i];
800 if (Map.MyGID(testGID)) {
801 LIDOfGID = Map.LID(testGID);
802 FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
803 if ((*x)[FirstEntryOfGID+GIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
804 if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
805 << (*x)[FirstEntryOfGID+GIDSize-1]
806 << " should = " << (VecValues[i]+VecValues[i]) << endl;
807 }
808 else
809 if (locerr!=1) err++; // Test for GID out of range error (=1)
810 }
811
812 // ========================================================================
813 // Test Local "My" versions of same routine (less complicated)
814 // ========================================================================
815 int * VecLIDs = new int[NumEntries];
816 VecLIDs[0] = testLID;
817 VecLIDs[1] = testLID+1; // Some pathological chance that these LIDs are not valid
818
819 VecValues[0] = 2.0; VecValues[1] = 4.0;
820 locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
821
822 for (i=0; i<NumEntries; i++) {
823 testLID = VecLIDs[i];
824 LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
825 FirstEntryOfLID = Map.FirstPointInElement(testLID);
826 if ((*x)[FirstEntryOfLID]!=VecValues[i]) err++;
827 if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
828 << (*x)[FirstEntryOfLID]
829 << " should = " << VecValues[i] << endl;
830 }
831
832 VecValues[0] = 4.0; VecValues[1] = 8.0;
833 locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
834
835 for (i=0; i<NumEntries; i++) {
836 testLID = VecLIDs[i];
837 LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
838 FirstEntryOfLID = Map.FirstPointInElement(testLID);
839 if ((*x)[FirstEntryOfLID+LIDSize-1]!=VecValues[i]) err++;
840 if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
841 << (*x)[FirstEntryOfLID+LIDSize-1]
842 << " should = " << VecValues[i] << endl;
843 }
844
845 VecValues[0] = 1.0; VecValues[1] = 1.0;
846 locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
847 locerr = x->SumIntoMyValues(NumEntries, VecValues, VecLIDs);
848
849 for (i=0; i<NumEntries; i++) {
850 testLID = VecLIDs[i];
851 LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
852 FirstEntryOfLID = Map.FirstPointInElement(testLID);
853 if ((*x)[FirstEntryOfLID]!=(VecValues[i]+VecValues[i])) err++;
854 if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
855 << (*x)[FirstEntryOfLID]
856 << " should = " << (VecValues[i]+VecValues[i]) << endl;
857 }
858
859 VecValues[0] = 2.0; VecValues[1] = 4.0;
860 locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
861 locerr = x->SumIntoMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
862
863 for (i=0; i<NumEntries; i++) {
864 testLID = VecLIDs[i];
865 LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
866 FirstEntryOfLID = Map.FirstPointInElement(testLID);
867 if ((*x)[FirstEntryOfLID+LIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
868 if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
869 << (*x)[FirstEntryOfLID+LIDSize-1]
870 << " should = " << (VecValues[i]+VecValues[i]) << endl;
871 }
872
873 delete [] VecValues;
874 delete [] VecGIDs;
875 delete [] VecLIDs;
876
877 return(ierr);
878}
879
880int BadResidual(bool verbose, double * Residual, int NumVectors)
881{
882 double threshold = 5.0E-6;
883 int ierr = 0;
884 for (int i=0; i<NumVectors; i++) {
885 if (Residual[i]>threshold) {
886 ierr = 1;// Output will be more useful after returning from this method
887 if (verbose) cout << endl << " Residual[" << i <<"] = " << Residual[i];
888 }
889 }
890 if (verbose)
891 if (ierr==0) cout << "\t Checked OK" << endl;
892
893 return(ierr);
894}
895
896// This version tests to make sure residuals are large (when we want vectors to be different)
897int BadResidual1(bool verbose, double * Residual, int NumVectors)
898{
899 double threshold = 5.0E-6;
900 int ierr = 0;
901 for (int i=0; i<NumVectors; i++) {
902 if (Residual[i]<threshold) {
903 ierr = 1;// Output will be more useful after returning from this method
904 if (verbose) cout << endl << " Residual[" << i <<"] = " << Residual[i] << " Should be larger";
905 }
906 }
907 if (verbose)
908 if (ierr==0) cout << "\t Checked OK" << endl;
909
910 return(ierr);
911}
@ Insert
#define EPETRA_MIN(x, y)
@ View
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)
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 BadResidual(bool verbose, double *Residual, int NumVectors)
int BadResidual1(bool verbose, double *Residual, int NumVectors)
int MatrixTests(const Epetra_BlockMap &Map, const Epetra_LocalMap &LocalMap, int NumVectors, bool verbose)
int MultiVectorTests(const Epetra_BlockMap &Map, int NumVectors, 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.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
int NumGlobalElements() const
Number of elements across all processors.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
int FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
virtual int NumProc() const =0
Returns total number of processes.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
virtual int MyPID() const =0
Return my process ID.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
Definition: Epetra_Import.h:63
Epetra_LocalMap: A class for replicating vectors and matrices across multiple processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
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 SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (GlobalRow, VectorIndex) location.
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 ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (GlobalRow, VectorIndex) location with ScalarValue.
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 ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (MyRow, VectorIndex) location with ScalarValue.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
int SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
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.
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int ReplaceGlobalValues(int NumEntries, const double *Values, const int *Indices)
Replace values in a vector with a given indexed list of values, indices are in global index space.
int ReplaceMyValues(int NumEntries, const double *Values, const int *Indices)
Replace values in a vector with a given indexed list of values, indices are in local index space.
int SumIntoMyValues(int NumEntries, const double *Values, const int *Indices)
Sum values into a vector with a given indexed list of values, indices are in local index space.
int SumIntoGlobalValues(int NumEntries, const double *Values, const int *Indices)
Sum values into a vector with a given indexed list of values, indices are in global index space.
#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)