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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #include "Epetra_BLAS.h"
00030 #include "ExecuteTestProblems.h"
00031 #include "BuildTestProblems.h"
00032 #include "Epetra_Comm.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_IntVector.h"
00035 #include "Epetra_Import.h"
00036 
00037   int MatrixTests(const Epetra_BlockMap & Map, const Epetra_LocalMap & LocalMap, int NumVectors,
00038           bool verbose)
00039   {
00040     const Epetra_Comm & Comm = Map.Comm();
00041     int ierr = 0, i;
00042     int IndexBase = 0;
00043     double *residual = new double[NumVectors];
00044 
00045     /* get ID of this processor */
00046 
00047 
00048     // Test GEMM first.  7 cases:
00049     
00050     //                                       Num
00051     //     OPERATIONS                        case  Notes
00052     // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No Comm needed) 
00053     // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
00054     // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no Comm needed)
00055 
00056     // ==================================================================
00057     // Case 1 through 4 (A, B, C all local) Strided and non-strided cases
00058     // ==================================================================
00059 
00060     // Construct MultiVectors
00061 
00062     {
00063     Epetra_MultiVector A(LocalMap, NumVectors);
00064     Epetra_MultiVector B(LocalMap, NumVectors);
00065     Epetra_LocalMap  Map2d(NumVectors, IndexBase, Comm);
00066     Epetra_MultiVector C(Map2d, NumVectors);
00067     Epetra_MultiVector C_GEMM(Map2d, NumVectors);
00068 
00069     double **App, **Bpp, **Cpp;
00070     
00071     Epetra_MultiVector *Ap, *Bp, *Cp;
00072 
00073     // For testing non-strided mode, create MultiVectors that are scattered throughout memory
00074 
00075     App = new double *[NumVectors];
00076     Bpp = new double *[NumVectors];
00077     Cpp = new double *[NumVectors];
00078     for (i=0; i<NumVectors; i++) App[i] = new double[A.MyLength()+i];
00079     for (i=0; i<NumVectors; i++) Bpp[i] = new double[B.MyLength()+i];
00080     for (i=0; i<NumVectors; i++) Cpp[i] = new double[C.MyLength()+i];
00081     
00082     Epetra_MultiVector A1(View, LocalMap, App, NumVectors);
00083     Epetra_MultiVector B1(View, LocalMap, Bpp, NumVectors);
00084     Epetra_MultiVector C1(View, Map2d, Cpp, NumVectors);
00085 
00086     for (int strided = 0; strided<2; strided++) {
00087       int ierr;
00088     // Loop through all trans cases using a variety of values for alpha and beta
00089     for (i=0; i<4; i++)  {
00090   ierr = 0;
00091   char transa = 'N'; if (i>1) transa = 'T';
00092   char transb = 'N'; if (i%2!=0) transb = 'T';
00093   double alpha = (double) i+1;
00094   double beta  = (double) (i/2);
00095   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00096   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00097   if (strided)
00098     {
00099       Ap = &A; Bp = &B; Cp = &C;
00100     }
00101   else
00102     {
00103       A.ExtractCopy(App); Ap = &A1;
00104       B.ExtractCopy(Bpp); Bp = &B1;
00105       C.ExtractCopy(Cpp); Cp = &C1;
00106     }
00107     
00108   ierr += Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
00109   ierr += Cp->Update(-1.0, C_GEMM, 1.0);
00110   ierr += Cp->Norm2(residual);
00111 
00112   if (verbose && ierr==0)
00113     {
00114       cout << "XXXXX Replicated Local MultiVector GEMM tests";
00115       if (strided)
00116       cout << " (Strided Multivectors)" << endl;
00117       else
00118       cout << " (Non-Strided Multivectors)" << endl;
00119       cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00120      <<", transb = " << transb;
00121     }
00122   if (ierr==0 && BadResidual(verbose,residual, NumVectors)) return(-1);
00123       }
00124 
00125       }
00126     for (i=0; i<NumVectors; i++)
00127       {
00128   delete [] App[i];
00129   delete [] Bpp[i];
00130   delete [] Cpp[i];
00131       }
00132     delete [] App;
00133     delete [] Bpp;
00134     delete [] Cpp;
00135     }
00136       
00137     // ====================================
00138     // Case 5  (A, B distributed C  local)
00139     // ====================================
00140 
00141     // Construct MultiVectors
00142   {
00143     Epetra_MultiVector A(Map, NumVectors);
00144     Epetra_MultiVector B(Map, NumVectors);
00145     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00146     Epetra_MultiVector C(Map2d, NumVectors);
00147     Epetra_MultiVector C_GEMM(Map2d, NumVectors);
00148 
00149     char transa = 'T';
00150     char transb = 'N';
00151     double alpha = 2.0;
00152     double beta  = 1.0;
00153     EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00154     ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
00155     ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00156     ierr += C.Update(-1.0, C_GEMM, 1.0);
00157     ierr += C.Norm2(residual);
00158 
00159     if (verbose && ierr==0)
00160       {
00161   cout << "XXXXX Generalized 2D dot product via GEMM call     " << endl;
00162   cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00163        <<", transb = " << transb;
00164       }
00165     if (BadResidual(verbose,residual, NumVectors)) return(-1);
00166     
00167     
00168   }      
00169     // ====================================
00170     // Case 6-7  (A, C distributed, B local)
00171     // ====================================
00172 
00173     // Construct MultiVectors
00174   {
00175     Epetra_MultiVector A(Map, NumVectors);
00176     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00177     Epetra_MultiVector B(Map2d, NumVectors);
00178     Epetra_MultiVector C(Map, NumVectors);
00179     Epetra_MultiVector C_GEMM(Map, NumVectors);
00180 
00181     for (i=0; i<2; i++)
00182       {
00183   char transa = 'N';
00184   char transb = 'N'; if (i>0) transb = 'T';
00185   double alpha = 2.0;
00186   double beta  = 1.1;
00187   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00188   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00189   ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00190   ierr += C.Update(-1.0, C_GEMM, 1.0);
00191   ierr += C.Norm2(residual);
00192   
00193   if (verbose)
00194     {
00195       cout << "XXXXX Generalized 2D vector update via GEMM call     " << endl;
00196       cout << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00197      <<", transb = " << transb;
00198     }
00199   if (BadResidual(verbose,residual, NumVectors)) return(-1);
00200       }
00201 
00202     delete [] residual;
00203     
00204     return(ierr);
00205   }
00206   }
00207 
00208 int MultiVectorTests(const Epetra_BlockMap & Map, int NumVectors, bool verbose)
00209 {
00210   const Epetra_Comm & Comm = Map.Comm();
00211   int ierr = 0, i;
00212   double *residual = new double[NumVectors];
00213   
00214   Epetra_BLAS BLAS;
00215   /* get number of processors and the name of this processor */
00216   
00217   // int NumProc = Comm.getNumProc();
00218   int MyPID   = Comm.MyPID();
00219   
00220   // Construct MultiVectors
00221   
00222   Epetra_MultiVector A(Map, NumVectors);
00223   Epetra_MultiVector sqrtA(Map, NumVectors);
00224   Epetra_MultiVector B(Map, NumVectors);
00225   Epetra_MultiVector C(Map, NumVectors);
00226   Epetra_MultiVector C_alphaA(Map, NumVectors);
00227   Epetra_MultiVector C_alphaAplusB(Map, NumVectors);
00228   Epetra_MultiVector C_plusB(Map, NumVectors);
00229   Epetra_MultiVector Weights(Map, NumVectors);
00230   
00231   // Construct double vectors
00232   double *dotvec_AB   = new double[NumVectors];
00233   double *norm1_A     = new double[NumVectors];
00234   double *norm2_sqrtA = new double[NumVectors];
00235   double *norminf_A = new double[NumVectors];
00236   double *normw_A = new double[NumVectors];
00237   double *minval_A = new double[NumVectors];
00238   double *maxval_A = new double[NumVectors];
00239   double *meanval_A = new double[NumVectors];
00240   
00241   // Generate data 
00242 
00243   
00244   EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
00245   double alpha = 2.0;
00246   BuildMultiVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
00247            C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A, 
00248            normw_A, Weights, minval_A, maxval_A, meanval_A);
00249 
00250   int err = 0;
00251   if (verbose) cout << "XXXXX Testing alpha * A     ";
00252   // Test alpha*A
00253   Epetra_MultiVector alphaA(A); // Copy of A
00254   EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
00255   EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
00256   EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
00257 
00258   if (err) ierr += err;
00259   else {
00260     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00261   }
00262 
00263   err = 0;
00264   if (verbose) cout << "XXXXX Testing C = alpha * A + B      ";
00265   // Test alpha*A + B
00266   Epetra_MultiVector alphaAplusB(A); // Copy of A
00267   EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
00268   EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
00269   EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
00270   
00271   if (err) ierr += err;
00272   else {
00273     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00274   }
00275 
00276   err = 0;
00277   if (verbose) cout << "XXXXX Testing C += B      ";
00278   // Test + B
00279   Epetra_MultiVector plusB(C); // Copy of C
00280   EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
00281   EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
00282   EPETRA_TEST_ERR(plusB.Norm2(residual),err);
00283   
00284   if (err) ierr += err;
00285   else {
00286     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00287   }
00288 
00289   err = 0;
00290   if (verbose) cout << "XXXXX Testing A.dotProd(B)     ";
00291   // Test A.dotvec(B)
00292   double *dotvec = residual;
00293   EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
00294   BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
00295   
00296   if (err) ierr += err;
00297   else {
00298     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00299   }
00300 
00301   err = 0;
00302   if (verbose) cout << "XXXXX Testing norm1_A      ";
00303   // Test A.norm1()
00304   double *norm1 = residual;
00305   EPETRA_TEST_ERR(A.Norm1(norm1),err);
00306   BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
00307   
00308   if (err) ierr += err;
00309   else {
00310     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00311   }  
00312 
00313   err = 0;
00314   if (verbose) cout << "XXXXX Testing norm2_sqrtA     ";
00315   // Test sqrtA.norm2()
00316   double *norm2 = residual;
00317   EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
00318   BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
00319   
00320   if (err) ierr += err;
00321   else {
00322     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00323   }
00324 
00325   err = 0;
00326   if (verbose) cout << "XXXXX Testing norminf_A     ";
00327   // Test A.norminf()
00328   double *norminf = residual;
00329   EPETRA_TEST_ERR(A.NormInf(norminf),err);
00330   BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
00331   
00332   if (err) ierr += err;
00333   else {
00334     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00335   }
00336   
00337   err = 0;
00338   if (verbose) cout << "XXXXX Testing normw_A     ";
00339   // Test A.NormWeighted()
00340   double *normw = residual;
00341   EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
00342   BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
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 minval_A     ";
00351   // Test A.MinValue()
00352   double *minval = residual;
00353   EPETRA_TEST_ERR(A.MinValue(minval),err);
00354   BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
00355   
00356   if (err) ierr += err;
00357   else {
00358     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00359   }
00360 
00361   err = 0;
00362   if (verbose) cout << "XXXXX Testing maxval_A     ";
00363   // Test A.MaxValue()
00364   double *maxval = residual;
00365   EPETRA_TEST_ERR(A.MaxValue(maxval),err);
00366   BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
00367 
00368   if (err) ierr += err;
00369   else {
00370     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00371   }
00372 
00373   err = 0;
00374   if (verbose) cout << "XXXXX Testing meanval_A     ";
00375   // Test A.MeanValue()
00376   double *meanval = residual;
00377   EPETRA_TEST_ERR(A.MeanValue(meanval),err);
00378   BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
00379   
00380   if (err) ierr += err;
00381   else {
00382     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00383   }
00384 
00385   err = 0;
00386   if (verbose) cout << "XXXXX Testing abs_A     ";
00387   // Test A.Abs()
00388   Epetra_MultiVector Abs_A = A;
00389   EPETRA_TEST_ERR(Abs_A.Abs(A),err);
00390   EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
00391   EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
00392 
00393   if (err) ierr += err;
00394   else {
00395     EPETRA_TEST_ERR(BadResidual(verbose,residual, NumVectors),ierr);
00396   }
00397   
00398   err = 0;
00399   if (verbose) cout << "XXXXX Testing random_A (Test1) ";
00400   // Test A.Random()
00401   Epetra_MultiVector Rand1_A(A);
00402   Epetra_MultiVector Rand2_A(A);
00403   EPETRA_TEST_ERR(Rand1_A.Random(),err);
00404   EPETRA_TEST_ERR(Rand2_A.Random(),err);
00405   // Rand2_A = Rand1_A - Rand2_A (should be nonzero since Random() should give different vectors > 0)
00406   EPETRA_TEST_ERR(Rand2_A.Update(1.0, Rand1_A, -1.0),err); 
00407   EPETRA_TEST_ERR(Rand2_A.Norm2(residual),err);
00408 
00409   if (err) ierr += err;
00410   else {
00411     EPETRA_TEST_ERR(BadResidual1(verbose,residual, NumVectors),ierr);
00412   }
00413 
00414   err = 0;
00415   if (verbose) cout << "XXXXX Testing random_A (Test2) ";
00416 
00417   // Next test that each column of the multivector is different from all other columns by testing the first value
00418   // of each vector against the first value of every other vector.
00419   int randvalsdiffer = 1; // Assume they all differ
00420   for (i=0; i< NumVectors; i++) 
00421     for (int j=i+1; j<NumVectors; j++) 
00422       if (Rand1_A[i][0]==Rand1_A[j][0]) randvalsdiffer = 0; // make false if equal
00423   int allrandvals = 0;
00424   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00425 
00426   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error
00427   int locerr = err;
00428   Comm.MinAll(&locerr, &err, 1);
00429 
00430   if (verbose)
00431     if (err==0) {
00432       cout << "\t Checked OK" << endl;
00433     } else {
00434       cout << "\t Checked Failed" << endl;
00435     }
00436 
00437   err = 0;
00438   if (verbose) cout << "XXXXX Testing random_A (Test3) ";
00439 
00440   // Next test that the first element on each processor of the first column of Rand1_A is different from all others
00441   // First we will gather them all to PE 0
00442   
00443 
00444   Epetra_Map RandstartsMap(-1, 1, 0, Comm); // This Map has a single element on each PE
00445   int itmp = 0;
00446   int nproc = Comm.NumProc();
00447   if (MyPID==0) itmp = nproc;
00448   Epetra_Map AllrandstartsMap(nproc, itmp, 0, Comm); // Map has NumProc elements on PE 0, none elsewhere
00449   Epetra_MultiVector Randstarts(RandstartsMap, NumVectors);
00450   Epetra_MultiVector Allrandstarts(AllrandstartsMap, NumVectors);
00451   for (i=0; i< NumVectors; i++) Randstarts[i][0] = Rand1_A[i][0]; // Load first value of local multivector
00452 
00453   Epetra_Import Randimporter(AllrandstartsMap,RandstartsMap);
00454   EPETRA_TEST_ERR(Allrandstarts.Import(Randstarts,Randimporter,Insert),err);
00455   // cout << "Randstarts = " << Randstarts << endl << "Allrandstarts = " << Allrandstarts << endl;
00456   // Allrandstarts now contains the first values for each local section of Rand1_A.
00457   // Next test that this is true.
00458   randvalsdiffer = 1; // Assume they all differ
00459   if (MyPID==0) {
00460     for (i=0; i< NumVectors; i++) 
00461       for (int irand=0; irand<nproc; irand++)
00462   for (int jrand=irand+1; jrand<nproc; jrand++) 
00463     if (Allrandstarts[i][irand]==Allrandstarts[i][jrand]) randvalsdiffer = 0; // make false if equal
00464   }
00465   allrandvals = 0;
00466   Comm.MinAll(&randvalsdiffer, &allrandvals, 1); // get min of all values across all processors
00467 
00468   EPETRA_TEST_ERR(1-allrandvals, err); // If allrandvals is anything but 1, this will cause an error 
00469   locerr = err;
00470   Comm.MinAll(&locerr, &err, 1);
00471   if (verbose)
00472     if (err==0) {
00473       cout << "\t Checked OK" << endl;
00474     } else {
00475       cout << "\t Checked Failed" << endl;
00476     }
00477 
00478   // Delete everything
00479   
00480   delete [] dotvec_AB;
00481   delete [] norm1_A;
00482   delete [] norm2_sqrtA;
00483   delete [] norminf_A;
00484   delete [] normw_A;
00485   delete [] minval_A;
00486   delete [] maxval_A;
00487   delete [] meanval_A;
00488   delete [] residual;
00489 
00490   //*******************************************************************
00491   // Post-construction modification tests
00492   //*******************************************************************
00493   
00494   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a multivector"
00495         <<endl<<endl;
00496 
00497   err = 0;
00498 
00499   Epetra_MultiVector X(Map, NumVectors);
00500   X.Random();
00501 
00502   // Pick middle range values for GID, LID and Vector Index
00503   int testGID = Map.NumGlobalElements()/2;
00504   int testVecIndex = NumVectors/2;
00505 
00506   int GIDSize = 1;
00507   int LIDOfGID = 0;
00508   int FirstEntryOfGID = 0;
00509 
00510   if (Map.MyGID(testGID)) {
00511     LIDOfGID = Map.LID(testGID);
00512     GIDSize = Map.ElementSize(LIDOfGID);
00513     FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00514   }
00515 
00516   // ========================================================================
00517   // Test int ReplaceGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00518   // ========================================================================
00519 
00520   double newGIDValue = 4.0;
00521   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00522 
00523   if (Map.MyGID(testGID)) {
00524     if (X[testVecIndex][FirstEntryOfGID]!=newGIDValue) err++;
00525     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00526           <<  X[testVecIndex][FirstEntryOfGID]
00527           << " should = " << newGIDValue << endl;
00528   }
00529   else
00530     if (locerr!=1) err++; // Test for GID out of range error (=1)
00531 
00532   // ========================================================================
00533   // Test int ReplaceGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00534   // ========================================================================
00535   newGIDValue = 8.0;
00536   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00537 
00538   if (Map.MyGID(testGID)) {
00539     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=newGIDValue) err++;
00540     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00541           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00542           << " should = " << newGIDValue << endl;
00543   }
00544   else
00545     if (locerr!=1) err++; // Test for GID out of range error (=1)
00546   
00547   // ========================================================================
00548   // Test int SumIntoGlobalValue (int GlobalRow, int VectorIndex, double ScalarValue)
00549   // ========================================================================
00550 
00551   newGIDValue = 1.0;
00552   locerr = X.ReplaceGlobalValue(testGID, testVecIndex, newGIDValue);
00553   locerr = X.SumIntoGlobalValue(testGID, testVecIndex, newGIDValue);
00554   if (Map.MyGID(testGID)) {
00555     if (X[testVecIndex][FirstEntryOfGID]!=(newGIDValue+newGIDValue)) err++;
00556     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID<<"] = "
00557           <<  X[testVecIndex][FirstEntryOfGID]
00558           << " should = " << newGIDValue << endl;
00559   }
00560   else 
00561     if (locerr!=1) err++; // Test for GID out of range error (=1)
00562     
00563   // ========================================================================
00564   // Test int SumIntoGlobalValue (int GlobalRow, intBlockRowOffset, int VectorIndex, double ScalarValue)
00565   // ========================================================================
00566     
00567   newGIDValue = 1.0;
00568   locerr = X.ReplaceGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00569   locerr = X.SumIntoGlobalValue(testGID, GIDSize-1, testVecIndex, newGIDValue);
00570 
00571   if (Map.MyGID(testGID)) {
00572     if (X[testVecIndex][FirstEntryOfGID+GIDSize-1]!=(newGIDValue+newGIDValue)) err++;
00573     if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfGID+GIDSize-1<<"] = "
00574           <<  X[testVecIndex][FirstEntryOfGID+GIDSize-1]
00575           << " should = " << newGIDValue << endl;
00576   }
00577   else
00578     if (locerr!=1) err++; // Test for GID out of range error (=1)
00579 
00580   // ========================================================================
00581   // Test Local "My" versions of same routine (less complicated)
00582   // ========================================================================
00583   
00584   // Pick middle range values for LID
00585   int testLID = Map.NumMyElements()/2;
00586 
00587   int LIDSize = Map.ElementSize(testLID);
00588   int FirstEntryOfLID = Map.FirstPointInElement(testLID);
00589 
00590 
00591   double newLIDValue = 4.0;
00592   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00593 
00594   if (X[testVecIndex][FirstEntryOfLID]!=newLIDValue) err++;
00595   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00596         <<  X[testVecIndex][FirstEntryOfLID]
00597         << " should = " << newLIDValue << endl;
00598   
00599   newLIDValue = 8.0;
00600   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00601   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=newLIDValue) err++;
00602   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00603         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00604         << " should = " << newLIDValue << endl;
00605   newLIDValue = 1.0;
00606   locerr = X.ReplaceMyValue(testLID, testVecIndex, newLIDValue);
00607   locerr = X.SumIntoMyValue(testLID, testVecIndex, newLIDValue);
00608   if (X[testVecIndex][FirstEntryOfLID]!=(newLIDValue+newLIDValue)) err++;
00609   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID<<"] = "
00610         <<  X[testVecIndex][FirstEntryOfLID]
00611         << " should = " << newLIDValue << endl;
00612   newLIDValue = 2.0;
00613   locerr = X.ReplaceMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00614   locerr = X.SumIntoMyValue(testLID, LIDSize-1, testVecIndex, newLIDValue);
00615   if (verbose) cout << "X["<<testVecIndex<<"]["<<FirstEntryOfLID+LIDSize-1<<"] = "
00616         <<  X[testVecIndex][FirstEntryOfLID+LIDSize-1]
00617         << " should = " << newLIDValue << endl;
00618   if (X[testVecIndex][FirstEntryOfLID+LIDSize-1]!=(newLIDValue+newLIDValue)) err++;
00619 
00620   ierr += err;
00621 
00622   // ========================================================================
00623   // Test Post-construction modification of an Epetra_Vector using a vector
00624   // our multivector X
00625   // ========================================================================
00626 
00627   if (verbose) cout <<  "\n\nXXXXX Testing Post-construction modification of a vector"
00628         << endl << endl;
00629 
00630   Epetra_Vector * x = X(testVecIndex);
00631 
00632   int NumEntries = 2;
00633   double * VecValues = new double[NumEntries];
00634   int * VecGIDs = new int[NumEntries];
00635   VecGIDs[0] = testGID;
00636   VecGIDs[1] = testGID+1; // Some pathological chance that these GIDs are not valid
00637 
00638   // ========================================================================
00639   // Test int ReplaceGlobalValues (int NumEntries, double *Values, int *Indices)
00640   // ========================================================================
00641 
00642   VecValues[0] = 2.0; VecValues[1] = 4.0;
00643   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00644 
00645   for (i=0; i<NumEntries; i++) {
00646     testGID = VecGIDs[i];
00647     if (Map.MyGID(testGID)) {
00648       LIDOfGID = Map.LID(testGID);
00649       GIDSize = EPETRA_MIN(GIDSize,Map.ElementSize(LIDOfGID)); // Need this value below
00650       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00651       if ((*x)[FirstEntryOfGID]!=VecValues[i]) err++;
00652       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00653       << (*x)[FirstEntryOfGID] 
00654       << " should = " << VecValues[i] << endl;
00655     }
00656     else
00657       if (locerr!=1) err++; // Test for GID out of range error (=1)
00658   }
00659 
00660 
00661   // ========================================================================
00662   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00663   // ========================================================================
00664 
00665   VecValues[0] = 4.0; VecValues[1] = 8.0;
00666   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00667 
00668   for (i=0; i<NumEntries; i++) {
00669     testGID = VecGIDs[i];
00670     if (Map.MyGID(testGID)) {
00671       LIDOfGID = Map.LID(testGID);
00672       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00673       if ((*x)[FirstEntryOfGID+GIDSize-1]!=VecValues[i]) err++;
00674       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00675       << (*x)[FirstEntryOfGID+GIDSize-1] 
00676       << " should = " << VecValues[i] << endl;
00677     }
00678     else
00679       if (locerr!=1) err++; // Test for GID out of range error (=1)
00680   }
00681 
00682   // ========================================================================
00683   // Test int SumIntoGlobalValues (int NumEntries, double *Values, int *Indices)
00684   // ========================================================================
00685 
00686   VecValues[0] = 1.0; VecValues[1] = 2.0;
00687   locerr = x->ReplaceGlobalValues(NumEntries, VecValues, VecGIDs);
00688   locerr = x->SumIntoGlobalValues(NumEntries, VecValues, VecGIDs);
00689 
00690   for (i=0; i<NumEntries; i++) {
00691     testGID = VecGIDs[i];
00692     if (Map.MyGID(testGID)) {
00693       LIDOfGID = Map.LID(testGID);
00694       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00695       if ((*x)[FirstEntryOfGID]!=(VecValues[i]+VecValues[i])) err++;
00696       if (verbose) cout << "x["<<FirstEntryOfGID<<"] = "
00697       << (*x)[FirstEntryOfGID] 
00698       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00699     }
00700     else
00701       if (locerr!=1) err++; // Test for GID out of range error (=1)
00702   }
00703   // ========================================================================
00704   // Test int ReplaceGlobalValues (int NumEntries, int BlockOffset, double *Values, int *Indices)
00705   // ========================================================================
00706 
00707   VecValues[0] = 1.0; VecValues[1] = 2.0;
00708   locerr = x->ReplaceGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00709   locerr = x->SumIntoGlobalValues(NumEntries, GIDSize-1, VecValues, VecGIDs);
00710 
00711   for (i=0; i<NumEntries; i++) {
00712     testGID = VecGIDs[i];
00713     if (Map.MyGID(testGID)) {
00714       LIDOfGID = Map.LID(testGID);
00715       FirstEntryOfGID = Map.FirstPointInElement(LIDOfGID);
00716       if ((*x)[FirstEntryOfGID+GIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00717       if (verbose) cout << "x["<<FirstEntryOfGID+GIDSize-1<<"] = "
00718       << (*x)[FirstEntryOfGID+GIDSize-1] 
00719       << " should = " << (VecValues[i]+VecValues[i]) << endl;
00720     }
00721     else
00722       if (locerr!=1) err++; // Test for GID out of range error (=1)
00723   }
00724 
00725   // ========================================================================
00726   // Test Local "My" versions of same routine (less complicated)
00727   // ========================================================================
00728   int * VecLIDs = new int[NumEntries];
00729   VecLIDs[0] = testLID;
00730   VecLIDs[1] = testLID+1; // Some pathological chance that these LIDs are not valid
00731 
00732   VecValues[0] = 2.0; VecValues[1] = 4.0;
00733   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00734 
00735   for (i=0; i<NumEntries; i++) {
00736     testLID = VecLIDs[i];
00737     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00738     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00739     if ((*x)[FirstEntryOfLID]!=VecValues[i]) err++;
00740     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00741           << (*x)[FirstEntryOfLID] 
00742           << " should = " << VecValues[i] << endl;
00743   }
00744 
00745   VecValues[0] = 4.0; VecValues[1] = 8.0;
00746   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00747 
00748   for (i=0; i<NumEntries; i++) {
00749     testLID = VecLIDs[i];
00750     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00751     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00752     if ((*x)[FirstEntryOfLID+LIDSize-1]!=VecValues[i]) err++;
00753     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00754           << (*x)[FirstEntryOfLID+LIDSize-1] 
00755           << " should = " << VecValues[i] << endl;
00756   }
00757 
00758   VecValues[0] = 1.0; VecValues[1] = 1.0;
00759   locerr = x->ReplaceMyValues(NumEntries, VecValues, VecLIDs);
00760   locerr = x->SumIntoMyValues(NumEntries, VecValues, VecLIDs);
00761 
00762   for (i=0; i<NumEntries; i++) {
00763     testLID = VecLIDs[i];
00764     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00765     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00766     if ((*x)[FirstEntryOfLID]!=(VecValues[i]+VecValues[i])) err++;
00767     if (verbose) cout << "x["<<FirstEntryOfLID<<"] = "
00768           << (*x)[FirstEntryOfLID] 
00769           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00770   }
00771 
00772   VecValues[0] = 2.0; VecValues[1] = 4.0;
00773   locerr = x->ReplaceMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00774   locerr = x->SumIntoMyValues(NumEntries, LIDSize-1, VecValues, VecLIDs);
00775 
00776   for (i=0; i<NumEntries; i++) {
00777     testLID = VecLIDs[i];
00778     LIDSize = EPETRA_MIN(LIDSize,Map.ElementSize(testLID)); // Need this value below
00779     FirstEntryOfLID = Map.FirstPointInElement(testLID);
00780     if ((*x)[FirstEntryOfLID+LIDSize-1]!=(VecValues[i]+VecValues[i])) err++;
00781     if (verbose) cout << "x["<<FirstEntryOfLID+LIDSize-1<<"] = "
00782           << (*x)[FirstEntryOfLID+LIDSize-1] 
00783           << " should = " << (VecValues[i]+VecValues[i]) << endl;
00784   }
00785 
00786     delete [] VecValues;
00787     delete [] VecGIDs;
00788     delete [] VecLIDs;
00789 
00790   return(ierr);
00791 }
00792 
00793 int BadResidual(bool verbose, double * Residual, int NumVectors)
00794 {
00795   double threshold = 5.0E-6;
00796   int ierr = 0;
00797   for (int i=0; i<NumVectors; i++) {
00798     if (Residual[i]>threshold) {
00799       ierr = 1;// Output will be more useful after returning from this method
00800       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i];
00801     }
00802   }
00803   if (verbose)
00804     if (ierr==0) cout << "\t Checked OK" << endl;
00805   
00806   return(ierr);
00807 }
00808 
00809 // This version tests to make sure residuals are large (when we want vectors to be different)
00810 int BadResidual1(bool verbose, double * Residual, int NumVectors)
00811 {
00812   double threshold = 5.0E-6;
00813   int ierr = 0;
00814   for (int i=0; i<NumVectors; i++) {
00815     if (Residual[i]<threshold) {
00816       ierr = 1;// Output will be more useful after returning from this method
00817       if (verbose) cout << endl << "     Residual[" << i <<"] = " << Residual[i] << "  Should be larger";
00818     }
00819   }
00820   if (verbose)
00821     if (ierr==0) cout << "\t Checked OK" << endl;
00822   
00823   return(ierr);
00824 }

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1