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 2001 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       int ierr;
00102     // Loop through all trans cases using a variety of values for alpha and beta
00103     for (i=0; i<4; i++)  {
00104   ierr = 0;
00105   char transa = 'N'; if (i>1) transa = 'T';
00106   char transb = 'N'; if (i%2!=0) transb = 'T';
00107   double alpha = (double) i+1;
00108   double beta  = (double) (i/2);
00109   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00110   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
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   ierr += Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
00123   ierr += Cp->Update(-1.0, C_GEMM, 1.0);
00124   ierr += Cp->Norm2(residual);
00125 
00126   if (verbose && ierr==0)
00127     {
00128       cout << "XXXXX Replicated Local MultiVector GEMM tests";
00129       if (strided)
00130       cout << " (Strided Multivectors)" << endl;
00131       else
00132       cout << " (Non-Strided Multivectors)" << endl;
00133       cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00134      <<", transb = " << transb;
00135     }
00136   if (ierr==0 && BadResidual(verbose,residual, NumVectors)) return(-1);
00137       }
00138 
00139       }
00140     for (i=0; i<NumVectors; i++)
00141       {
00142   delete [] App[i];
00143   delete [] Bpp[i];
00144   delete [] Cpp[i];
00145       }
00146     delete [] App;
00147     delete [] Bpp;
00148     delete [] Cpp;
00149     }
00150       
00151     // ====================================
00152     // Case 5  (A, B distributed C  local)
00153     // ====================================
00154 
00155     // Construct MultiVectors
00156   {
00157     Epetra_MultiVector A(Map, NumVectors);
00158     Epetra_MultiVector B(Map, NumVectors);
00159     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00160     Epetra_MultiVector C(Map2d, NumVectors);
00161     Epetra_MultiVector C_GEMM(Map2d, NumVectors);
00162 
00163     char transa = 'T';
00164     char transb = 'N';
00165     double alpha = 2.0;
00166     double beta  = 1.0;
00167     EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00168     ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
00169     ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00170     ierr += C.Update(-1.0, C_GEMM, 1.0);
00171     ierr += C.Norm2(residual);
00172 
00173     if (verbose && ierr==0)
00174       {
00175   cout << "XXXXX Generalized 2D dot product via GEMM call     " << endl;
00176   cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00177        <<", transb = " << transb;
00178       }
00179     if (BadResidual(verbose,residual, NumVectors)) return(-1);
00180     
00181     
00182   }      
00183     // ====================================
00184     // Case 6-7  (A, C distributed, B local)
00185     // ====================================
00186 
00187     // Construct MultiVectors
00188   {
00189     Epetra_MultiVector A(Map, NumVectors);
00190     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00191     Epetra_MultiVector B(Map2d, NumVectors);
00192     Epetra_MultiVector C(Map, NumVectors);
00193     Epetra_MultiVector C_GEMM(Map, NumVectors);
00194 
00195     for (i=0; i<2; i++)
00196       {
00197   char transa = 'N';
00198   char transb = 'N'; if (i>0) transb = 'T';
00199   double alpha = 2.0;
00200   double beta  = 1.1;
00201   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00202   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00203   ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00204   ierr += C.Update(-1.0, C_GEMM, 1.0);
00205   ierr += C.Norm2(residual);
00206   
00207   if (verbose)
00208     {
00209       cout << "XXXXX Generalized 2D vector update via GEMM call     " << endl;
00210       cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00211      <<", transb = " << transb;
00212     }
00213   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00214       }
00215 
00216     
00217   }
00218     // ====================================
00219     // LocalMap Tests
00220     // ====================================
00221 
00222     // Construct MultiVectors
00223   {
00224     
00225         int localLength = 10;
00226         double *localMinValue = new double[localLength];
00227         double *localMaxValue = new double[localLength];
00228         double *localNorm1 = new double[localLength];
00229         double *localDot = new double[localLength];
00230         double *localNorm2 = new double[localLength];
00231         double *localMeanValue = new double[localLength];
00232         Epetra_LocalMap MapSmall(localLength, IndexBase, Comm);
00233         Epetra_MultiVector A(MapSmall, NumVectors);
00234 
00235         double doubleLocalLength = (double) localLength;
00236         for (int j=0; j< NumVectors; j++) {
00237           for (i=0; i< localLength-1; i++) A[j][i] = (double) (i+1);
00238           A[j][localLength-1] = (double) (localLength+j); // Only the last value differs across multivectors
00239           localMinValue[j] = A[j][0]; // Increasing values
00240           localMaxValue[j] = A[j][localLength-1];
00241           localNorm1[j] = (doubleLocalLength-1.0)*(doubleLocalLength)/2.0+A[j][localLength-1];
00242           localDot[j] = (doubleLocalLength-1.0)*(doubleLocalLength)*(2.0*(doubleLocalLength-1.0)+1.0)/6.0+A[j][localLength-1]*A[j][localLength-1];
00243           localNorm2[j] = std::sqrt(localDot[j]);
00244           localMeanValue[j] = localNorm1[j]/doubleLocalLength;
00245         }
00246   ierr += A.MinValue(residual);
00247         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMinValue[j]);
00248   if (verbose) cout << "XXXXX MinValue" << endl;
00249   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00250 
00251   ierr += A.MaxValue(residual);
00252         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMaxValue[j]);
00253   if (verbose) cout << "XXXXX MaxValue" << endl;
00254   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00255 
00256   ierr += A.Norm1(residual);
00257         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm1[j]);
00258   if (verbose) cout << "XXXXX Norm1" << endl;
00259   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00260 
00261   ierr += A.Dot(A,residual);
00262         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localDot[j]);
00263   if (verbose) cout << "XXXXX Dot" << endl;
00264   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00265 
00266   ierr += A.Norm2(residual);
00267         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localNorm2[j]);
00268   if (verbose) cout << "XXXXX Norm2" << endl;
00269   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00270 
00271   ierr += A.MeanValue(residual);
00272         for (int j=0; j<NumVectors; j++) residual[j] = std::abs(residual[j] - localMeanValue[j]);
00273   if (verbose) cout << "XXXXX MeanValue" << endl;
00274   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00275 
00276         delete [] localMinValue;
00277         delete [] localMaxValue;
00278         delete [] localNorm1;
00279         delete [] localDot;
00280         delete [] localNorm2;
00281         delete [] localMeanValue;
00282 
00283   }
00284 
00285     delete [] residual;
00286     
00287     return(ierr);
00288   }
00289 
00290 int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
00291 {
00292   const Epetra_Comm & Comm = Map.Comm();
00293   int ierr = 0, i;
00294   double *residual = new double[NumVectors];
00295   
00296   Epetra_BLAS BLAS;
00297   /* get number of processors and the name of this processor */
00298   
00299   // int NumProc = Comm.getNumProc();
00300   int MyPID   = Comm.MyPID();
00301   
00302   // Construct MultiVectors
00303   
00304   Epetra_MultiVector A(Map, NumVectors);
00305   Epetra_MultiVector sqrtA(Map, NumVectors);
00306   Epetra_MultiVector B(Map, NumVectors);
00307   Epetra_MultiVector C(Map, NumVectors);
00308   Epetra_MultiVector C_alphaA(Map, NumVectors);
00309   Epetra_MultiVector C_alphaAplusB(Map, NumVectors);
00310   Epetra_MultiVector C_plusB(Map, NumVectors);
00311   Epetra_MultiVector Weights(Map, NumVectors);
00312   
00313   // Construct double vectors
00314   double *dotvec_AB   = new double[NumVectors];
00315   double *norm1_A     = new double[NumVectors];
00316   double *norm2_sqrtA = new double[NumVectors];
00317   double *norminf_A = new double[NumVectors];
00318   double *normw_A = new double[NumVectors];
00319   double *minval_A = new double[NumVectors];
00320   double *maxval_A = new double[NumVectors];
00321   double *meanval_A = new double[NumVectors];
00322   
00323   // Generate data 
00324 
00325   
00326   EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
00327   double alpha = 2.0;
00328   BuildMultiVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
00329            C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A, 
00330            normw_A, Weights, minval_A, maxval_A, meanval_A);
00331 
00332   int err = 0;
00333   if (verbose) cout << "XXXXX Testing alpha * A     ";
00334   // Test alpha*A
00335   Epetra_MultiVector alphaA(A); // Copy of A
00336   EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
00337   EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
00338   EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
00339 
00340   if (err) ierr += err;
00341   else {
00342     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00343   }
00344 
00345   err = 0;
00346   if (verbose) cout << "XXXXX Testing C = alpha * A + B      ";
00347   // Test alpha*A + B
00348   Epetra_MultiVector alphaAplusB(A); // Copy of A
00349   EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
00350   EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
00351   EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
00352   
00353   if (err) ierr += err;
00354   else {
00355     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00356   }
00357 
00358   err = 0;
00359   if (verbose) cout << "XXXXX Testing C += B      ";
00360   // Test + B
00361   Epetra_MultiVector plusB(C); // Copy of C
00362   EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
00363   EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
00364   EPETRA_TEST_ERR(plusB.Norm2(residual),err);
00365   
00366   if (err) ierr += err;
00367   else {
00368     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00369   }
00370 
00371   err = 0;
00372   if (verbose) cout << "XXXXX Testing A.dotProd(B)     ";
00373   // Test A.dotvec(B)
00374   double *dotvec = residual;
00375   EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
00376   BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
00377   
00378   if (err) ierr += err;
00379   else {
00380     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00381   }
00382 
00383   err = 0;
00384   if (verbose) cout << "XXXXX Testing norm1_A      ";
00385   // Test A.norm1()
00386   double *norm1 = residual;
00387   EPETRA_TEST_ERR(A.Norm1(norm1),err);
00388   BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
00389   
00390   if (err) ierr += err;
00391   else {
00392     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00393   }  
00394 
00395   err = 0;
00396   if (verbose) cout << "XXXXX Testing norm2_sqrtA     ";
00397   // Test sqrtA.norm2()
00398   double *norm2 = residual;
00399   EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
00400   BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
00401   
00402   if (err) ierr += err;
00403   else {
00404     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00405   }
00406 
00407   err = 0;
00408   if (verbose) cout << "XXXXX Testing norminf_A     ";
00409   // Test A.norminf()
00410   double *norminf = residual;
00411   EPETRA_TEST_ERR(A.NormInf(norminf),err);
00412   BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
00413   
00414   if (err) ierr += err;
00415   else {
00416     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00417   }
00418   
00419   err = 0;
00420   if (verbose) cout << "XXXXX Testing normw_A     ";
00421   // Test A.NormWeighted()
00422   double *normw = residual;
00423   EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
00424   BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
00425   
00426   if (err) ierr += err;
00427   else {
00428     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00429   }
00430   
00431   err = 0;
00432   if (verbose) cout << "XXXXX Testing minval_A     ";
00433   // Test A.MinValue()
00434   double *minval = residual;
00435   EPETRA_TEST_ERR(A.MinValue(minval),err);
00436   BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
00437   
00438   if (err) ierr += err;
00439   else {
00440     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00441   }
00442 
00443   err = 0;
00444   if (verbose) cout << "XXXXX Testing maxval_A     ";
00445   // Test A.MaxValue()
00446   double *maxval = residual;
00447   EPETRA_TEST_ERR(A.MaxValue(maxval),err);
00448   BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
00449 
00450   if (err) ierr += err;
00451   else {
00452     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00453   }
00454 
00455   err = 0;
00456   if (verbose) cout << "XXXXX Testing meanval_A     ";
00457   // Test A.MeanValue()
00458   double *meanval = residual;
00459   EPETRA_TEST_ERR(A.MeanValue(meanval),err);
00460   BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
00461   
00462   if (err) ierr += err;
00463   else {
00464     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00465   }
00466 
00467   err = 0;
00468   if (verbose) cout << "XXXXX Testing abs_A     ";
00469   // Test A.Abs()
00470   Epetra_MultiVector Abs_A = A;
00471   EPETRA_TEST_ERR(Abs_A.Abs(A),err);
00472   EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
00473   EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
00474 
00475   if (err) ierr += err;
00476   else {
00477     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00478   }
00479   
00480   err = 0;
00481   if (verbose) cout << "XXXXX Testing random_A (Test1) ";
00482   // Test A.Random()
00483   Epetra_MultiVector Rand1_A(A);
00484   Epetra_MultiVector Rand2_A(A);
00485   EPETRA_TEST_ERR(Rand1_A.Random(),err);
00486   EPETRA_TEST_ERR(Rand2_A.Random(),err);
00487   // Rand2_A = Rand1_A - Rand2_A (should be nonzero since Random() should give different vectors > 0)
00488   EPETRA_TEST_ERR(Rand2_A.Update(1.0, Rand1_A, -1.0),err); 
00489   EPETRA_TEST_ERR(Rand2_A.Norm2(residual),err);
00490 
00491   if (err) ierr += err;
00492   else {
00493     EPETRA_TEST_ERR(BadResidual1(verbose,residual, NumVectors),ierr);
00494   }
00495 
00496   err = 0;
00497   if (verbose) cout << "XXXXX Testing random_A (Test2) ";
00498 
00499   // Next test that each column of the multivector is different from all other columns by testing the first value
00500   // of each vector against the first value of every other vector.
00501   int randvalsdiffer = 1; // Assume they all differ
00502   for (i=0; i< NumVectors; i++) 
00503     for (int j=i+1; j<NumVectors; j++) 
00504       if (Rand1_A[i][0]==Rand1_A[j][0]) randvalsdiffer = 0; // make false if equal
00505   int allrandvals = 0;
00506   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00507 
00508   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
00509   int locerr = err;
00510   Comm.MinAll(&locerr, &err, 1);
00511 
00512   if (verbose)
00513     if (err==0) {
00514       cout << "\t Checked OK" << endl;
00515     } else {
00516       cout << "\t Checked Failed" << endl;
00517     }
00518 
00519   err = 0;
00520   if (verbose) cout << "XXXXX Testing random_A (Test3) ";
00521 
00522   // Next test that the first element on each processor of the first column of Rand1_A is different from all others
00523   // First we will gather them all to PE 0
00524   
00525 
00526   Epetra_Map RandstartsMap(-1, 1, 0, Comm); // This Map has a single element on each PE
00527   int itmp = 0;
00528   int nproc = Comm.NumProc();
00529   if (MyPID==0) itmp = nproc;
00530   Epetra_Map AllrandstartsMap(nproc, itmp, 0, Comm); // Map has NumProc elements on PE 0, none elsewhere
00531   Epetra_MultiVector Randstarts(RandstartsMap, NumVectors);
00532   Epetra_MultiVector Allrandstarts(AllrandstartsMap, NumVectors);
00533   for (i=0; i< NumVectors; i++) Randstarts[i][0] = Rand1_A[i][0]; // Load first value of local multivector
00534 
00535   Epetra_Import Randimporter(AllrandstartsMap,RandstartsMap);
00536   EPETRA_TEST_ERR(Allrandstarts.Import(Randstarts,Randimporter,Insert),err);
00537   // cout << "Randstarts = " << Randstarts << endl << "Allrandstarts = " << Allrandstarts << endl;
00538   // Allrandstarts now contains the first values for each local section of Rand1_A.
00539   // Next test that this is true.
00540   randvalsdiffer = 1; // Assume they all differ
00541   if (MyPID==0) {
00542     for (i=0; i< NumVectors; i++) 
00543       for (int irand=0; irand<nproc; irand++)
00544   for (int jrand=irand+1; jrand<nproc; jrand++) 
00545     if (Allrandstarts[i][irand]==Allrandstarts[i][jrand]) randvalsdiffer = 0; // make false if equal
00546   }
00547   allrandvals = 0;
00548   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00549 
00550   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error 
00551   locerr = err;
00552   Comm.MinAll(&locerr, &err, 1);
00553   if (verbose)
00554     if (err==0) {
00555       cout << "\t Checked OK" << endl;
00556     } else {
00557       cout << "\t Checked Failed" << endl;
00558     }
00559 
00560   // Delete everything
00561   
00562   delete [] dotvec_AB;
00563   delete [] norm1_A;
00564   delete [] norm2_sqrtA;
00565   delete [] norminf_A;
00566   delete [] normw_A;
00567   delete [] minval_A;
00568   delete [] maxval_A;
00569   delete [] meanval_A;
00570   delete [] residual;
00571 
00572   //*******************************************************************
00573   // Post-construction modification tests
00574   //*******************************************************************
00575   
00576   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a multivector"
00577         <<endl<<endl;
00578 
00579   err = 0;
00580 
00581   Epetra_MultiVector X(Map, NumVectors);
00582   X.Random();
00583 
00584   // Pick middle range values for GID, LID and Vector Index
00585   int testGID = Map.NumGlobalElements()/2;
00586   int testVecIndex = NumVectors/2;
00587 
00588   int GIDSize = 1;
00589   int LIDOfGID = 0;
00590   int FirstEntryOfGID = 0;
00591 
00592   if (Map.MyGID(testGID)) {
00593     LIDOfGID = Map.LID(testGID);
00594     GIDSize = Map.ElementSize(LIDOfGID);
00595     FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00596   }
00597 
00598   // ========================================================================
00599   // Test int ReplaceGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00600   // ========================================================================
00601 
00602   double newGIDValue = 4.0;
00603   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00604 
00605   if (Map.MyGID(testGID)) {
00606     if (X[testVecIndex][FirstEntryOfGID]!=newGIDValue) err++;
00607     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00608           <<  X[testVecIndex][FirstEntryOfGID]
00609           << " should = " << newGIDValue << endl;
00610   }
00611   else
00612     if (locerr!=1) err++; // Test for GID out of range error (=1)
00613 
00614   // ========================================================================
00615   // Test int ReplaceGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00616   // ========================================================================
00617   newGIDValue = 8.0;
00618   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00619 
00620   if (Map.MyGID(testGID)) {
00621     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=newGIDValue) err++;
00622     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00623           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00624           << " should = " << newGIDValue << endl;
00625   }
00626   else
00627     if (locerr!=1) err++; // Test for GID out of range error (=1)
00628   
00629   // ========================================================================
00630   // Test int SumIntoGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00631   // ========================================================================
00632 
00633   newGIDValue = 1.0;
00634   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00635   locerr = X.SumIntoGlobalValue(testGID, testVecIndex, newGIDValue);
00636   if (Map.MyGID(testGID)) {
00637     if (X[testVecIndex][FirstEntryOfGID]!=(newGIDValue+newGIDValue)) err++;
00638     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00639           <<  X[testVecIndex][FirstEntryOfGID]
00640           << " should = " << newGIDValue << endl;
00641   }
00642   else 
00643     if (locerr!=1) err++; // Test for GID out of range error (=1)
00644     
00645   // ========================================================================
00646   // Test int SumIntoGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00647   // ========================================================================
00648     
00649   newGIDValue = 1.0;
00650   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00651   locerr = X.SumIntoGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00652 
00653   if (Map.MyGID(testGID)) {
00654     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=(newGIDValue+newGIDValue)) err++;
00655     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00656           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00657           << " should = " << newGIDValue << endl;
00658   }
00659   else
00660     if (locerr!=1) err++; // Test for GID out of range error (=1)
00661 
00662   // ========================================================================
00663   // Test Local "My" versions of same routine (less complicated)
00664   // ========================================================================
00665   
00666   // Pick middle range values for LID
00667   int testLID = Map.NumMyElements()/2;
00668 
00669   int LIDSize = Map.ElementSize(testLID);
00670   int FirstEntryOfLID = Map.FirstPointInElement(testLID);
00671 
00672 
00673   double newLIDValue = 4.0;
00674   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00675 
00676   if (X[testVecIndex][FirstEntryOfLID]!=newLIDValue) err++;
00677   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00678         <<  X[testVecIndex][FirstEntryOfLID]
00679         << " should = " << newLIDValue << endl;
00680   
00681   newLIDValue = 8.0;
00682   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00683   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=newLIDValue) err++;
00684   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00685         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00686         << " should = " << newLIDValue << endl;
00687   newLIDValue = 1.0;
00688   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00689   locerr = X.SumIntoMyValue(testLID, testVecIndex, newLIDValue);
00690   if (X[testVecIndex][FirstEntryOfLID]!=(newLIDValue+newLIDValue)) err++;
00691   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00692         <<  X[testVecIndex][FirstEntryOfLID]
00693         << " should = " << newLIDValue << endl;
00694   newLIDValue = 2.0;
00695   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00696   locerr = X.SumIntoMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00697   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00698         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00699         << " should = " << newLIDValue << endl;
00700   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=(newLIDValue+newLIDValue)) err++;
00701 
00702   ierr += err;
00703 
00704   // ========================================================================
00705   // Test Post-construction modification of an Epetra_Vector using a vector
00706   // our multivector X
00707   // ========================================================================
00708 
00709   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a vector"
00710         << endl << endl;
00711 
00712   Epetra_Vector * x = X(testVecIndex);
00713 
00714   int NumEntries = 2;
00715   double * VecValues = new double[NumEntries];
00716   int * VecGIDs = new int[NumEntries];
00717   VecGIDs[0] = testGID;
00718   VecGIDs[1] = testGID+1; // Some pathological chance that these GIDs are not valid
00719 
00720   // ========================================================================
00721   // Test int ReplaceGlobalValues (int NumEntries, double *Values, int *Indices)
00722   // ========================================================================
00723 
00724   VecValues[0] = 2.0; VecValues[1] = 4.0;
00725   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00726 
00727   for (i=0; i<NumEntries; i++) {
00728     testGID = VecGIDs[i];
00729     if (Map.MyGID(testGID)) {
00730       LIDOfGID = Map.LID(testGID);
00731       GIDSize = EPETRA_MIN(GIDSize,Map.ElementSize(LIDOfGID)); // Need this value below
00732       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00733       if ((*x)[FirstEntryOfGID]!=VecValues[i]) err++;
00734       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00735       << (*x)[FirstEntryOfGID] 
00736       << " should = " << VecValues[i] << endl;
00737     }
00738     else
00739       if (locerr!=1) err++; // Test for GID out of range error (=1)
00740   }
00741 
00742 
00743   // ========================================================================
00744   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00745   // ========================================================================
00746 
00747   VecValues[0] = 4.0; VecValues[1] = 8.0;
00748   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00749 
00750   for (i=0; i<NumEntries; i++) {
00751     testGID = VecGIDs[i];
00752     if (Map.MyGID(testGID)) {
00753       LIDOfGID = Map.LID(testGID);
00754       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00755       if ((*x)[FirstEntryOfGID+GIDSize-1]!=VecValues[i]) err++;
00756       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00757       << (*x)[FirstEntryOfGID+GIDSize-1] 
00758       << " should = " << VecValues[i] << endl;
00759     }
00760     else
00761       if (locerr!=1) err++; // Test for GID out of range error (=1)
00762   }
00763 
00764   // ========================================================================
00765   // Test int SumIntoGlobalValues (int NumEntries, double *Values, int *Indices)
00766   // ========================================================================
00767 
00768   VecValues[0] = 1.0; VecValues[1] = 2.0;
00769   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00770   locerr = x->SumIntoGlobalValues(NumEntries, VecValues, VecGIDs);
00771 
00772   for (i=0; i<NumEntries; i++) {
00773     testGID = VecGIDs[i];
00774     if (Map.MyGID(testGID)) {
00775       LIDOfGID = Map.LID(testGID);
00776       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00777       if ((*x)[FirstEntryOfGID]!=(VecValues[i]+VecValues[i])) err++;
00778       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00779       << (*x)[FirstEntryOfGID] 
00780       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00781     }
00782     else
00783       if (locerr!=1) err++; // Test for GID out of range error (=1)
00784   }
00785   // ========================================================================
00786   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00787   // ========================================================================
00788 
00789   VecValues[0] = 1.0; VecValues[1] = 2.0;
00790   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00791   locerr = x->SumIntoGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00792 
00793   for (i=0; i<NumEntries; i++) {
00794     testGID = VecGIDs[i];
00795     if (Map.MyGID(testGID)) {
00796       LIDOfGID = Map.LID(testGID);
00797       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00798       if ((*x)[FirstEntryOfGID+GIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00799       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00800       << (*x)[FirstEntryOfGID+GIDSize-1] 
00801       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00802     }
00803     else
00804       if (locerr!=1) err++; // Test for GID out of range error (=1)
00805   }
00806 
00807   // ========================================================================
00808   // Test Local "My" versions of same routine (less complicated)
00809   // ========================================================================
00810   int * VecLIDs = new int[NumEntries];
00811   VecLIDs[0] = testLID;
00812   VecLIDs[1] = testLID+1; // Some pathological chance that these LIDs are not valid
00813 
00814   VecValues[0] = 2.0; VecValues[1] = 4.0;
00815   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00816 
00817   for (i=0; i<NumEntries; i++) {
00818     testLID = VecLIDs[i];
00819     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00820     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00821     if ((*x)[FirstEntryOfLID]!=VecValues[i]) err++;
00822     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00823           << (*x)[FirstEntryOfLID] 
00824           << " should = " << VecValues[i] << endl;
00825   }
00826 
00827   VecValues[0] = 4.0; VecValues[1] = 8.0;
00828   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00829 
00830   for (i=0; i<NumEntries; i++) {
00831     testLID = VecLIDs[i];
00832     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00833     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00834     if ((*x)[FirstEntryOfLID+LIDSize-1]!=VecValues[i]) err++;
00835     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00836           << (*x)[FirstEntryOfLID+LIDSize-1] 
00837           << " should = " << VecValues[i] << endl;
00838   }
00839 
00840   VecValues[0] = 1.0; VecValues[1] = 1.0;
00841   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00842   locerr = x->SumIntoMyValues(NumEntries, VecValues, VecLIDs);
00843 
00844   for (i=0; i<NumEntries; i++) {
00845     testLID = VecLIDs[i];
00846     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00847     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00848     if ((*x)[FirstEntryOfLID]!=(VecValues[i]+VecValues[i])) err++;
00849     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00850           << (*x)[FirstEntryOfLID] 
00851           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00852   }
00853 
00854   VecValues[0] = 2.0; VecValues[1] = 4.0;
00855   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00856   locerr = x->SumIntoMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00857 
00858   for (i=0; i<NumEntries; i++) {
00859     testLID = VecLIDs[i];
00860     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00861     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00862     if ((*x)[FirstEntryOfLID+LIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00863     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00864           << (*x)[FirstEntryOfLID+LIDSize-1] 
00865           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00866   }
00867 
00868     delete [] VecValues;
00869     delete [] VecGIDs;
00870     delete [] VecLIDs;
00871 
00872   return(ierr);
00873 }
00874 
00875 int BadResidual(bool verbose, double * Residual, int NumVectors)
00876 {
00877   double threshold = 5.0E-6;
00878   int ierr = 0;
00879   for (int i=0; i<NumVectors; i++) {
00880     if (Residual[i]>threshold) {
00881       ierr = 1;// Output will be more useful after returning from this method
00882       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i];
00883     }
00884   }
00885   if (verbose)
00886     if (ierr==0) cout << "\t Checked OK" << endl;
00887   
00888   return(ierr);
00889 }
00890 
00891 // This version tests to make sure residuals are large (when we want vectors to be different)
00892 int BadResidual1(bool verbose, double * Residual, int NumVectors)
00893 {
00894   double threshold = 5.0E-6;
00895   int ierr = 0;
00896   for (int i=0; i<NumVectors; i++) {
00897     if (Residual[i]<threshold) {
00898       ierr = 1;// Output will be more useful after returning from this method
00899       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i] << "  Should be larger";
00900     }
00901   }
00902   if (verbose)
00903     if (ierr==0) cout << "\t Checked OK" << endl;
00904   
00905   return(ierr);
00906 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines