Epetra Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
294 int 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 
880 int 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)
897 int 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 }
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumGlobalElements() const
Number of elements across all processors.
Epetra_Map: A class for partitioning vectors and matrices.
Definition: Epetra_Map.h:119
int Abs(const Epetra_MultiVector &A)
Puts element-wise absolute values of input Multi-vector in target.
int Random()
Set multi-vector values to random numbers.
int ReplaceMyValue(int MyRow, int VectorIndex, double ScalarValue)
Replace current value at the specified (MyRow, VectorIndex) location with ScalarValue.
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
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)
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
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)
int NormWeighted(const Epetra_MultiVector &Weights, double *Result) const
Compute Weighted 2-norm (RMS Norm) of each vector in multi-vector.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int MatrixTests(const Epetra_BlockMap &Map, const Epetra_LocalMap &LocalMap, int NumVectors, bool verbose)
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
#define EPETRA_MIN(x, y)
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
Definition: Epetra_Import.h:63
int ExtractCopy(double *A, int MyLDA) const
Put multi-vector values into user-provided two-dimensional array.
virtual int MyPID() const =0
Return my process ID.
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 BadResidual1(bool verbose, double *Residual, int NumVectors)
int Dot(const Epetra_MultiVector &A, double *Result) const
Computes dot product of each corresponding pair of vectors.
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 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 SumIntoGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (GlobalRow, VectorIndex) location.
Epetra_BLAS: The Epetra BLAS Wrapper Class.
Definition: Epetra_BLAS.h:70
int NumMyElements() const
Number of elements on the calling processor.
int MultiVectorTests(const Epetra_Map &Map, int NumVectors, bool verbose)
int Scale(double ScalarValue)
Scale the current values of a multi-vector, this = ScalarValue*this.
bool Residual(int N, int NRHS, double *A, int LDA, bool Transpose, double *X, int LDX, double *B, int LDB, double *resid)
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition: Epetra_Comm.h:73
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 FirstPointInElement(int LID) const
Returns the requested entry in the FirstPointInElementList; see FirstPointInElementList() for details...
int NormInf(double *Result) const
Compute Inf-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 SumIntoMyValue(int MyRow, int VectorIndex, double ScalarValue)
Adds ScalarValue to existing value at the specified (MyRow, VectorIndex) location.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int MeanValue(double *Result) const
Compute mean (average) value of each vector in multi-vector.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
int MinValue(double *Result) const
Compute minimum value of each vector in multi-vector.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
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 Norm1(double *Result) const
Compute 1-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.
virtual int NumProc() const =0
Returns total number of processes.
int Norm2(double *Result) const
Compute 2-norm of each vector in multi-vector.
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 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_LocalMap: A class for replicating vectors and matrices across multiple processors.
int BadResidual(bool verbose, double *Residual, int NumVectors)