Epetra Package Browser (Single Doxygen Collection) Development
MultiVector/ExecuteTestProblems.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: Linear Algebra Services Package
00005 //                 Copyright 2011 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00043 #include "Epetra_BLAS.h"
00044 #include "ExecuteTestProblems.h"
00045 #include "BuildTestProblems.h"
00046 #include "Epetra_Comm.h"
00047 #include "Epetra_Vector.h"
00048 #include "Epetra_IntVector.h"
00049 #include "Epetra_Import.h"
00050 
00051   int MatrixTests(const Epetra_BlockMap & Map, const Epetra_LocalMap & LocalMap, int NumVectors,
00052           bool verbose)
00053   {
00054     const Epetra_Comm & Comm = Map.Comm();
00055     int ierr = 0, i;
00056     int IndexBase = 0;
00057     double *residual = new double[NumVectors];
00058 
00059     /* get ID of this processor */
00060 
00061 
00062     // Test GEMM first.  7 cases:
00063 
00064     //                                       Num
00065     //     OPERATIONS                        case  Notes
00066     // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No Comm needed)
00067     // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
00068     // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no Comm needed)
00069 
00070     // ==================================================================
00071     // Case 1 through 4 (A, B, C all local) Strided and non-strided cases
00072     // ==================================================================
00073 
00074     // Construct MultiVectors
00075 
00076     {
00077     Epetra_MultiVector A(LocalMap, NumVectors);
00078     Epetra_MultiVector B(LocalMap, NumVectors);
00079     Epetra_LocalMap  Map2d(NumVectors, IndexBase, Comm);
00080     Epetra_MultiVector C(Map2d, NumVectors);
00081     Epetra_MultiVector C_GEMM(Map2d, NumVectors);
00082 
00083     double **App, **Bpp, **Cpp;
00084 
00085     Epetra_MultiVector *Ap, *Bp, *Cp;
00086 
00087     // For testing non-strided mode, create MultiVectors that are scattered throughout memory
00088 
00089     App = new double *[NumVectors];
00090     Bpp = new double *[NumVectors];
00091     Cpp = new double *[NumVectors];
00092     for (i=0; i<NumVectors; i++) App[i] = new double[A.MyLength()+i];
00093     for (i=0; i<NumVectors; i++) Bpp[i] = new double[B.MyLength()+i];
00094     for (i=0; i<NumVectors; i++) Cpp[i] = new double[C.MyLength()+i];
00095 
00096     Epetra_MultiVector A1(View, LocalMap, App, NumVectors);
00097     Epetra_MultiVector B1(View, LocalMap, Bpp, NumVectors);
00098     Epetra_MultiVector C1(View, Map2d, Cpp, NumVectors);
00099 
00100     for (int strided = 0; strided<2; strided++) {
00101 
00102     // Loop through all trans cases using a variety of values for alpha and beta
00103     for (i=0; i<4; i++)  {
00104   char transa = 'N'; if (i>1) transa = 'T';
00105   char transb = 'N'; if (i%2!=0) transb = 'T';
00106   double alpha = (double) i+1;
00107   double beta  = (double) (i/2);
00108   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00109   int localierr = BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00110   if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
00111     if (strided)
00112       {
00113         Ap = &A; Bp = &B; Cp = &C;
00114       }
00115     else
00116       {
00117         A.ExtractCopy(App); Ap = &A1;
00118         B.ExtractCopy(Bpp); Bp = &B1;
00119         C.ExtractCopy(Cpp); Cp = &C1;
00120       }
00121   
00122     localierr = Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
00123     if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
00124       ierr += Cp->Update(-1.0, C_GEMM, 1.0);
00125       ierr += Cp->Norm2(residual);
00126   
00127       if (verbose)
00128         {
00129     cout << "XXXXX Replicated Local MultiVector GEMM tests";
00130     if (strided)
00131       cout << " (Strided Multivectors)" << endl;
00132     else
00133       cout << " (Non-Strided Multivectors)" << endl;
00134     cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00135          <<", transb = " << transb;
00136         }
00137       if (BadResidual(verbose,residual, NumVectors)) return(-1);
00138     }
00139   }
00140       }
00141 
00142       }
00143     for (i=0; i<NumVectors; i++)
00144       {
00145   delete [] App[i];
00146   delete [] Bpp[i];
00147   delete [] Cpp[i];
00148       }
00149     delete [] App;
00150     delete [] Bpp;
00151     delete [] Cpp;
00152     }
00153 
00154     // ====================================
00155     // Case 5  (A, B distributed C  local)
00156     // ====================================
00157 
00158     // Construct MultiVectors
00159   {
00160     Epetra_MultiVector A(Map, NumVectors);
00161     Epetra_MultiVector B(Map, NumVectors);
00162     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00163     Epetra_MultiVector C(Map2d, NumVectors);
00164     Epetra_MultiVector C_GEMM(Map2d, NumVectors);
00165 
00166     char transa = 'T';
00167     char transb = 'N';
00168     double alpha = 2.0;
00169     double beta  = 1.0;
00170     EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00171     ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
00172     int localierr = C.Multiply(transa, transb, alpha, A, B, beta);
00173     if (localierr!=-2) { // -2 means the shapes didn't match
00174       ierr += C.Update(-1.0, C_GEMM, 1.0);
00175       ierr += C.Norm2(residual);
00176 
00177       if (verbose)
00178   {
00179     cout << "XXXXX Generalized 2D dot product via GEMM call     " << endl;
00180     cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00181          <<", transb = " << transb;
00182   }
00183       if (BadResidual(verbose,residual, NumVectors)) return(-1);
00184     }
00185 
00186   }
00187     // ====================================
00188     // Case 6-7  (A, C distributed, B local)
00189     // ====================================
00190 
00191     // Construct MultiVectors
00192   {
00193     Epetra_MultiVector A(Map, NumVectors);
00194     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00195     Epetra_MultiVector B(Map2d, NumVectors);
00196     Epetra_MultiVector C(Map, NumVectors);
00197     Epetra_MultiVector C_GEMM(Map, NumVectors);
00198 
00199     for (i=0; i<2; i++)
00200       {
00201   char transa = 'N';
00202   char transb = 'N'; if (i>0) transb = 'T';
00203   double alpha = 2.0;
00204   double beta  = 1.1;
00205   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00206   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00207   ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00208   ierr += C.Update(-1.0, C_GEMM, 1.0);
00209   ierr += C.Norm2(residual);
00210   
00211   if (verbose)
00212     {
00213       cout << "XXXXX Generalized 2D vector update via GEMM call     " << endl;
00214       cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00215      <<", transb = " << transb;
00216     }
00217   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00218       }
00219 
00220 
00221   }
00222     // ====================================
00223     // LocalMap Tests
00224     // ====================================
00225 
00226     // Construct MultiVectors
00227   {
00228 
00229         int localLength = 10;
00230         double *localMinValue = new double[localLength];
00231         double *localMaxValue = new double[localLength];
00232         double *localNorm1 = new double[localLength];
00233         double *localDot = new double[localLength];
00234         double *localNorm2 = new double[localLength];
00235         double *localMeanValue = new double[localLength];
00236         Epetra_LocalMap MapSmall(localLength, IndexBase, Comm);
00237         Epetra_MultiVector A(MapSmall, NumVectors);
00238 
00239         double doubleLocalLength = (double) localLength;
00240         for (int j=0; j< NumVectors; j++) {
00241           for (i=0; i< localLength-1; i++) A[j][i] = (double) (i+1);
00242           A[j][localLength-1] = (double) (localLength+j); // Only the last value differs across multivectors
00243           localMinValue[j] = A[j][0]; // Increasing values
00244           localMaxValue[j] = A[j][localLength-1];
00245           localNorm1[j] = (doubleLocalLength-1.0)*(doubleLocalLength)/2.0+A[j][localLength-1];
00246           localDot[j] = (doubleLocalLength-1.0)*(doubleLocalLength)*(2.0*(doubleLocalLength-1.0)+1.0)/6.0+A[j][localLength-1]*A[j][localLength-1];
00247           localNorm2[j] = std::sqrt(localDot[j]);
00248           localMeanValue[j] = localNorm1[j]/doubleLocalLength;
00249         }
00250   ierr += A.MinValue(residual);
00251         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMinValue[j]);
00252   if (verbose) cout << "XXXXX MinValue" << endl;
00253   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00254 
00255   ierr += A.MaxValue(residual);
00256         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMaxValue[j]);
00257   if (verbose) cout << "XXXXX MaxValue" << endl;
00258   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00259 
00260   ierr += A.Norm1(residual);
00261         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm1[j]);
00262   if (verbose) cout << "XXXXX Norm1" << endl;
00263   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00264 
00265   ierr += A.Dot(A,residual);
00266         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localDot[j]);
00267   if (verbose) cout << "XXXXX Dot" << endl;
00268   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00269 
00270   ierr += A.Norm2(residual);
00271         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm2[j]);
00272   if (verbose) cout << "XXXXX Norm2" << endl;
00273   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00274 
00275   ierr += A.MeanValue(residual);
00276         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMeanValue[j]);
00277   if (verbose) cout << "XXXXX MeanValue" << endl;
00278   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00279 
00280         delete [] localMinValue;
00281         delete [] localMaxValue;
00282         delete [] localNorm1;
00283         delete [] localDot;
00284         delete [] localNorm2;
00285         delete [] localMeanValue;
00286 
00287   }
00288 
00289     delete [] residual;
00290 
00291     return(ierr);
00292   }
00293 
00294 int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
00295 {
00296   const Epetra_Comm & Comm = Map.Comm();
00297   int ierr = 0, i;
00298   double *residual = new double[NumVectors];
00299 
00300   Epetra_BLAS BLAS;
00301   /* get number of processors and the name of this processor */
00302 
00303   // int NumProc = Comm.getNumProc();
00304   int MyPID   = Comm.MyPID();
00305 
00306   // Construct MultiVectors
00307 
00308   Epetra_MultiVector A(Map, NumVectors);
00309   Epetra_MultiVector sqrtA(Map, NumVectors);
00310   Epetra_MultiVector B(Map, NumVectors);
00311   Epetra_MultiVector C(Map, NumVectors);
00312   Epetra_MultiVector C_alphaA(Map, NumVectors);
00313   Epetra_MultiVector C_alphaAplusB(Map, NumVectors);
00314   Epetra_MultiVector C_plusB(Map, NumVectors);
00315   Epetra_MultiVector Weights(Map, NumVectors);
00316 
00317   // Construct double vectors
00318   double *dotvec_AB   = new double[NumVectors];
00319   double *norm1_A     = new double[NumVectors];
00320   double *norm2_sqrtA = new double[NumVectors];
00321   double *norminf_A = new double[NumVectors];
00322   double *normw_A = new double[NumVectors];
00323   double *minval_A = new double[NumVectors];
00324   double *maxval_A = new double[NumVectors];
00325   double *meanval_A = new double[NumVectors];
00326 
00327   // Generate data
00328 
00329 
00330   EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
00331   double alpha = 2.0;
00332   BuildMultiVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
00333            C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A,
00334            normw_A, Weights, minval_A, maxval_A, meanval_A);
00335 
00336   int err = 0;
00337   if (verbose) cout << "XXXXX Testing alpha * A     ";
00338   // Test alpha*A
00339   Epetra_MultiVector alphaA(A); // Copy of A
00340   EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
00341   EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
00342   EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
00343 
00344   if (err) ierr += err;
00345   else {
00346     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00347   }
00348 
00349   err = 0;
00350   if (verbose) cout << "XXXXX Testing C = alpha * A + B      ";
00351   // Test alpha*A + B
00352   Epetra_MultiVector alphaAplusB(A); // Copy of A
00353   EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
00354   EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
00355   EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
00356 
00357   if (err) ierr += err;
00358   else {
00359     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00360   }
00361 
00362   err = 0;
00363   if (verbose) cout << "XXXXX Testing C += B      ";
00364   // Test + B
00365   Epetra_MultiVector plusB(C); // Copy of C
00366   EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
00367   EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
00368   EPETRA_TEST_ERR(plusB.Norm2(residual),err);
00369 
00370   if (err) ierr += err;
00371   else {
00372     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00373   }
00374 
00375   err = 0;
00376   if (verbose) cout << "XXXXX Testing A.dotProd(B)     ";
00377   // Test A.dotvec(B)
00378   double *dotvec = residual;
00379   EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
00380   BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
00381 
00382   if (err) ierr += err;
00383   else {
00384     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00385   }
00386 
00387   err = 0;
00388   if (verbose) cout << "XXXXX Testing norm1_A      ";
00389   // Test A.norm1()
00390   double *norm1 = residual;
00391   EPETRA_TEST_ERR(A.Norm1(norm1),err);
00392   BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
00393 
00394   if (err) ierr += err;
00395   else {
00396     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00397   }
00398 
00399   err = 0;
00400   if (verbose) cout << "XXXXX Testing norm2_sqrtA     ";
00401   // Test sqrtA.norm2()
00402   double *norm2 = residual;
00403   EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
00404   BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
00405 
00406   if (err) ierr += err;
00407   else {
00408     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00409   }
00410 
00411   err = 0;
00412   if (verbose) cout << "XXXXX Testing norminf_A     ";
00413   // Test A.norminf()
00414   double *norminf = residual;
00415   EPETRA_TEST_ERR(A.NormInf(norminf),err);
00416   BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
00417 
00418   if (err) ierr += err;
00419   else {
00420     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00421   }
00422 
00423   err = 0;
00424   if (verbose) cout << "XXXXX Testing normw_A     ";
00425   // Test A.NormWeighted()
00426   double *normw = residual;
00427   EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
00428   BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
00429 
00430   if (err) ierr += err;
00431   else {
00432     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00433   }
00434 
00435   err = 0;
00436   if (verbose) cout << "XXXXX Testing minval_A     ";
00437   // Test A.MinValue()
00438   double *minval = residual;
00439   EPETRA_TEST_ERR(A.MinValue(minval),err);
00440   BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
00441 
00442   if (err) ierr += err;
00443   else {
00444     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00445   }
00446 
00447   err = 0;
00448   if (verbose) cout << "XXXXX Testing maxval_A     ";
00449   // Test A.MaxValue()
00450   double *maxval = residual;
00451   EPETRA_TEST_ERR(A.MaxValue(maxval),err);
00452   BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
00453 
00454   if (err) ierr += err;
00455   else {
00456     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00457   }
00458 
00459   err = 0;
00460   if (verbose) cout << "XXXXX Testing meanval_A     ";
00461   // Test A.MeanValue()
00462   double *meanval = residual;
00463   EPETRA_TEST_ERR(A.MeanValue(meanval),err);
00464   BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
00465 
00466   if (err) ierr += err;
00467   else {
00468     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00469   }
00470 
00471   err = 0;
00472   if (verbose) cout << "XXXXX Testing abs_A     ";
00473   // Test A.Abs()
00474   Epetra_MultiVector Abs_A = A;
00475   EPETRA_TEST_ERR(Abs_A.Abs(A),err);
00476   EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
00477   EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
00478 
00479   if (err) ierr += err;
00480   else {
00481     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00482   }
00483 
00484   err = 0;
00485   if (verbose) cout << "XXXXX Testing random_A (Test1) ";
00486   // Test A.Random()
00487   Epetra_MultiVector Rand1_A(A);
00488   Epetra_MultiVector Rand2_A(A);
00489   EPETRA_TEST_ERR(Rand1_A.Random(),err);
00490   EPETRA_TEST_ERR(Rand2_A.Random(),err);
00491   // Rand2_A = Rand1_A - Rand2_A (should be nonzero since Random() should give different vectors > 0)
00492   EPETRA_TEST_ERR(Rand2_A.Update(1.0, Rand1_A, -1.0),err);
00493   EPETRA_TEST_ERR(Rand2_A.Norm2(residual),err);
00494 
00495   if (err) ierr += err;
00496   else {
00497     EPETRA_TEST_ERR(BadResidual1(verbose,residual, NumVectors),ierr);
00498   }
00499 
00500   err = 0;
00501   if (verbose) cout << "XXXXX Testing random_A (Test2) ";
00502 
00503   // Next test that each column of the multivector is different from all other columns by testing the first value
00504   // of each vector against the first value of every other vector.
00505   int randvalsdiffer = 1; // Assume they all differ
00506   for (i=0; i< NumVectors; i++)
00507     for (int j=i+1; j<NumVectors; j++)
00508       if (Rand1_A[i][0]==Rand1_A[j][0]) randvalsdiffer = 0; // make false if equal
00509   int allrandvals = 0;
00510   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00511 
00512   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
00513   int locerr = err;
00514   Comm.MinAll(&locerr, &err, 1);
00515 
00516   if (verbose) {
00517     if (err==0) {
00518       cout << "\t Checked OK" << endl;
00519     } else {
00520       cout << "\t Checked Failed" << endl;
00521     }
00522   }
00523   err = 0;
00524   if (verbose) cout << "XXXXX Testing random_A (Test3) ";
00525 
00526   // Next test that the first element on each processor of the first column of Rand1_A is different from all others
00527   // First we will gather them all to PE 0
00528 
00529 
00530   Epetra_Map RandstartsMap(-1, 1, 0, Comm); // This Map has a single element on each PE
00531   int itmp = 0;
00532   int nproc = Comm.NumProc();
00533   if (MyPID==0) itmp = nproc;
00534   Epetra_Map AllrandstartsMap(nproc, itmp, 0, Comm); // Map has NumProc elements on PE 0, none elsewhere
00535   Epetra_MultiVector Randstarts(RandstartsMap, NumVectors);
00536   Epetra_MultiVector Allrandstarts(AllrandstartsMap, NumVectors);
00537   for (i=0; i< NumVectors; i++) Randstarts[i][0] = Rand1_A[i][0]; // Load first value of local multivector
00538 
00539   Epetra_Import Randimporter(AllrandstartsMap,RandstartsMap);
00540   EPETRA_TEST_ERR(Allrandstarts.Import(Randstarts,Randimporter,Insert),err);
00541   // cout << "Randstarts = " << Randstarts << endl << "Allrandstarts = " << Allrandstarts << endl;
00542   // Allrandstarts now contains the first values for each local section of Rand1_A.
00543   // Next test that this is true.
00544   randvalsdiffer = 1; // Assume they all differ
00545   if (MyPID==0) {
00546     for (i=0; i< NumVectors; i++)
00547       for (int irand=0; irand<nproc; irand++)
00548   for (int jrand=irand+1; jrand<nproc; jrand++)
00549     if (Allrandstarts[i][irand]==Allrandstarts[i][jrand]) randvalsdiffer = 0; // make false if equal
00550   }
00551   allrandvals = 0;
00552   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00553 
00554   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
00555   locerr = err;
00556   Comm.MinAll(&locerr, &err, 1);
00557   if (verbose) {
00558     if (err==0) {
00559       cout << "\t Checked OK" << endl;
00560     } else {
00561       cout << "\t Checked Failed" << endl;
00562     }
00563   }
00564 
00565   // Delete everything
00566 
00567   delete [] dotvec_AB;
00568   delete [] norm1_A;
00569   delete [] norm2_sqrtA;
00570   delete [] norminf_A;
00571   delete [] normw_A;
00572   delete [] minval_A;
00573   delete [] maxval_A;
00574   delete [] meanval_A;
00575   delete [] residual;
00576 
00577   //*******************************************************************
00578   // Post-construction modification tests
00579   //*******************************************************************
00580 
00581   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a multivector"
00582         <<endl<<endl;
00583 
00584   err = 0;
00585 
00586   Epetra_MultiVector X(Map, NumVectors);
00587   X.Random();
00588 
00589   // Pick middle range values for GID, LID and Vector Index
00590   int testGID = Map.NumGlobalElements()/2;
00591   int testVecIndex = NumVectors/2;
00592 
00593   int GIDSize = 1;
00594   int LIDOfGID = 0;
00595   int FirstEntryOfGID = 0;
00596 
00597   if (Map.MyGID(testGID)) {
00598     LIDOfGID = Map.LID(testGID);
00599     GIDSize = Map.ElementSize(LIDOfGID);
00600     FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00601   }
00602 
00603   // ========================================================================
00604   // Test int ReplaceGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00605   // ========================================================================
00606 
00607   double newGIDValue = 4.0;
00608   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00609 
00610   if (Map.MyGID(testGID)) {
00611     if (X[testVecIndex][FirstEntryOfGID]!=newGIDValue) err++;
00612     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00613           <<  X[testVecIndex][FirstEntryOfGID]
00614           << " should = " << newGIDValue << endl;
00615   }
00616   else
00617     if (locerr!=1) err++; // Test for GID out of range error (=1)
00618 
00619   // ========================================================================
00620   // Test int ReplaceGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00621   // ========================================================================
00622   newGIDValue = 8.0;
00623   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00624 
00625   if (Map.MyGID(testGID)) {
00626     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=newGIDValue) err++;
00627     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00628           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00629           << " should = " << newGIDValue << endl;
00630   }
00631   else
00632     if (locerr!=1) err++; // Test for GID out of range error (=1)
00633 
00634   // ========================================================================
00635   // Test int SumIntoGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00636   // ========================================================================
00637 
00638   newGIDValue = 1.0;
00639   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00640   locerr = X.SumIntoGlobalValue(testGID, testVecIndex, newGIDValue);
00641   if (Map.MyGID(testGID)) {
00642     if (X[testVecIndex][FirstEntryOfGID]!=(newGIDValue+newGIDValue)) err++;
00643     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00644           <<  X[testVecIndex][FirstEntryOfGID]
00645           << " should = " << newGIDValue << endl;
00646   }
00647   else
00648     if (locerr!=1) err++; // Test for GID out of range error (=1)
00649 
00650   // ========================================================================
00651   // Test int SumIntoGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00652   // ========================================================================
00653 
00654   newGIDValue = 1.0;
00655   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00656   locerr = X.SumIntoGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00657 
00658   if (Map.MyGID(testGID)) {
00659     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=(newGIDValue+newGIDValue)) err++;
00660     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00661           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00662           << " should = " << newGIDValue << endl;
00663   }
00664   else
00665     if (locerr!=1) err++; // Test for GID out of range error (=1)
00666 
00667   // ========================================================================
00668   // Test Local "My" versions of same routine (less complicated)
00669   // ========================================================================
00670 
00671   // Pick middle range values for LID
00672   int testLID = Map.NumMyElements()/2;
00673 
00674   int LIDSize = Map.ElementSize(testLID);
00675   int FirstEntryOfLID = Map.FirstPointInElement(testLID);
00676 
00677 
00678   double newLIDValue = 4.0;
00679   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00680 
00681   if (X[testVecIndex][FirstEntryOfLID]!=newLIDValue) err++;
00682   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00683         <<  X[testVecIndex][FirstEntryOfLID]
00684         << " should = " << newLIDValue << endl;
00685 
00686   newLIDValue = 8.0;
00687   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00688   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=newLIDValue) err++;
00689   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00690         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00691         << " should = " << newLIDValue << endl;
00692   newLIDValue = 1.0;
00693   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00694   locerr = X.SumIntoMyValue(testLID, testVecIndex, newLIDValue);
00695   if (X[testVecIndex][FirstEntryOfLID]!=(newLIDValue+newLIDValue)) err++;
00696   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00697         <<  X[testVecIndex][FirstEntryOfLID]
00698         << " should = " << newLIDValue << endl;
00699   newLIDValue = 2.0;
00700   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00701   locerr = X.SumIntoMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00702   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00703         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00704         << " should = " << newLIDValue << endl;
00705   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=(newLIDValue+newLIDValue)) err++;
00706 
00707   ierr += err;
00708 
00709   // ========================================================================
00710   // Test Post-construction modification of an Epetra_Vector using a vector
00711   // our multivector X
00712   // ========================================================================
00713 
00714   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a vector"
00715         << endl << endl;
00716 
00717   Epetra_Vector * x = X(testVecIndex);
00718 
00719   int NumEntries = 2;
00720   double * VecValues = new double[NumEntries];
00721   int * VecGIDs = new int[NumEntries];
00722   VecGIDs[0] = testGID;
00723   VecGIDs[1] = testGID+1; // Some pathological chance that these GIDs are not valid
00724 
00725   // ========================================================================
00726   // Test int ReplaceGlobalValues (int NumEntries, double *Values, int *Indices)
00727   // ========================================================================
00728 
00729   VecValues[0] = 2.0; VecValues[1] = 4.0;
00730   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00731 
00732   for (i=0; i<NumEntries; i++) {
00733     testGID = VecGIDs[i];
00734     if (Map.MyGID(testGID)) {
00735       LIDOfGID = Map.LID(testGID);
00736       GIDSize = EPETRA_MIN(GIDSize,Map.ElementSize(LIDOfGID)); // Need this value below
00737       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00738       if ((*x)[FirstEntryOfGID]!=VecValues[i]) err++;
00739       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00740       << (*x)[FirstEntryOfGID]
00741       << " should = " << VecValues[i] << endl;
00742     }
00743     else
00744       if (locerr!=1) err++; // Test for GID out of range error (=1)
00745   }
00746 
00747 
00748   // ========================================================================
00749   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00750   // ========================================================================
00751 
00752   VecValues[0] = 4.0; VecValues[1] = 8.0;
00753   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00754 
00755   for (i=0; i<NumEntries; i++) {
00756     testGID = VecGIDs[i];
00757     if (Map.MyGID(testGID)) {
00758       LIDOfGID = Map.LID(testGID);
00759       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00760       if ((*x)[FirstEntryOfGID+GIDSize-1]!=VecValues[i]) err++;
00761       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00762       << (*x)[FirstEntryOfGID+GIDSize-1]
00763       << " should = " << VecValues[i] << endl;
00764     }
00765     else
00766       if (locerr!=1) err++; // Test for GID out of range error (=1)
00767   }
00768 
00769   // ========================================================================
00770   // Test int SumIntoGlobalValues (int NumEntries, double *Values, int *Indices)
00771   // ========================================================================
00772 
00773   VecValues[0] = 1.0; VecValues[1] = 2.0;
00774   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00775   locerr = x->SumIntoGlobalValues(NumEntries, VecValues, VecGIDs);
00776 
00777   for (i=0; i<NumEntries; i++) {
00778     testGID = VecGIDs[i];
00779     if (Map.MyGID(testGID)) {
00780       LIDOfGID = Map.LID(testGID);
00781       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00782       if ((*x)[FirstEntryOfGID]!=(VecValues[i]+VecValues[i])) err++;
00783       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00784       << (*x)[FirstEntryOfGID]
00785       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00786     }
00787     else
00788       if (locerr!=1) err++; // Test for GID out of range error (=1)
00789   }
00790   // ========================================================================
00791   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00792   // ========================================================================
00793 
00794   VecValues[0] = 1.0; VecValues[1] = 2.0;
00795   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00796   locerr = x->SumIntoGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00797 
00798   for (i=0; i<NumEntries; i++) {
00799     testGID = VecGIDs[i];
00800     if (Map.MyGID(testGID)) {
00801       LIDOfGID = Map.LID(testGID);
00802       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00803       if ((*x)[FirstEntryOfGID+GIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00804       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00805       << (*x)[FirstEntryOfGID+GIDSize-1]
00806       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00807     }
00808     else
00809       if (locerr!=1) err++; // Test for GID out of range error (=1)
00810   }
00811 
00812   // ========================================================================
00813   // Test Local "My" versions of same routine (less complicated)
00814   // ========================================================================
00815   int * VecLIDs = new int[NumEntries];
00816   VecLIDs[0] = testLID;
00817   VecLIDs[1] = testLID+1; // Some pathological chance that these LIDs are not valid
00818 
00819   VecValues[0] = 2.0; VecValues[1] = 4.0;
00820   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00821 
00822   for (i=0; i<NumEntries; i++) {
00823     testLID = VecLIDs[i];
00824     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00825     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00826     if ((*x)[FirstEntryOfLID]!=VecValues[i]) err++;
00827     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00828           << (*x)[FirstEntryOfLID]
00829           << " should = " << VecValues[i] << endl;
00830   }
00831 
00832   VecValues[0] = 4.0; VecValues[1] = 8.0;
00833   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00834 
00835   for (i=0; i<NumEntries; i++) {
00836     testLID = VecLIDs[i];
00837     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00838     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00839     if ((*x)[FirstEntryOfLID+LIDSize-1]!=VecValues[i]) err++;
00840     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00841           << (*x)[FirstEntryOfLID+LIDSize-1]
00842           << " should = " << VecValues[i] << endl;
00843   }
00844 
00845   VecValues[0] = 1.0; VecValues[1] = 1.0;
00846   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00847   locerr = x->SumIntoMyValues(NumEntries, VecValues, VecLIDs);
00848 
00849   for (i=0; i<NumEntries; i++) {
00850     testLID = VecLIDs[i];
00851     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00852     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00853     if ((*x)[FirstEntryOfLID]!=(VecValues[i]+VecValues[i])) err++;
00854     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00855           << (*x)[FirstEntryOfLID]
00856           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00857   }
00858 
00859   VecValues[0] = 2.0; VecValues[1] = 4.0;
00860   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00861   locerr = x->SumIntoMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00862 
00863   for (i=0; i<NumEntries; i++) {
00864     testLID = VecLIDs[i];
00865     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00866     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00867     if ((*x)[FirstEntryOfLID+LIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00868     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00869           << (*x)[FirstEntryOfLID+LIDSize-1]
00870           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00871   }
00872 
00873     delete [] VecValues;
00874     delete [] VecGIDs;
00875     delete [] VecLIDs;
00876 
00877   return(ierr);
00878 }
00879 
00880 int BadResidual(bool verbose, double * Residual, int NumVectors)
00881 {
00882   double threshold = 5.0E-6;
00883   int ierr = 0;
00884   for (int i=0; i<NumVectors; i++) {
00885     if (Residual[i]>threshold) {
00886       ierr = 1;// Output will be more useful after returning from this method
00887       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i];
00888     }
00889   }
00890   if (verbose)
00891     if (ierr==0) cout << "\t Checked OK" << endl;
00892 
00893   return(ierr);
00894 }
00895 
00896 // This version tests to make sure residuals are large (when we want vectors to be different)
00897 int BadResidual1(bool verbose, double * Residual, int NumVectors)
00898 {
00899   double threshold = 5.0E-6;
00900   int ierr = 0;
00901   for (int i=0; i<NumVectors; i++) {
00902     if (Residual[i]<threshold) {
00903       ierr = 1;// Output will be more useful after returning from this method
00904       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i] << "  Should be larger";
00905     }
00906   }
00907   if (verbose)
00908     if (ierr==0) cout << "\t Checked OK" << endl;
00909 
00910   return(ierr);
00911 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines