Epetra Package Browser (Single Doxygen Collection) Development
Vector/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 
00048 #ifdef HAVE_EPETRA_TEUCHOS
00049 #  include "Teuchos_VerboseObject.hpp"
00050 #endif
00051 
00052   int MatrixTests(const Epetra_BlockMap & Map, const Epetra_LocalMap & LocalMap, 
00053     bool verbose)
00054   {
00055 
00056 #ifdef HAVE_EPETRA_TEUCHOS
00057   Teuchos::RCP<Teuchos::FancyOStream>
00058     fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
00059   std::ostream &out = *fancyOut;
00060 #else
00061   std::ostream &out = std::cout;
00062 #endif
00063 
00064     int NumVectors = 1;
00065     const Epetra_Comm & Comm = Map.Comm();
00066     int ierr = 0, i;
00067     int IndexBase = 0;
00068     double *residual = new double[NumVectors];
00069 
00070     /* get ID of this processor */
00071 
00072     // Test GEMM first.  7 cases:
00073     
00074     //                                       Num
00075     //     OPERATIONS                        case  Notes
00076     // 1) C(local) = A^X(local) * B^X(local)  4   (X=Trans or Not, No Comm needed) 
00077     // 2) C(local) = A^T(distr) * B  (distr)  1   (2D dot product, replicate C)
00078     // 3) C(distr) = A  (distr) * B^X(local)  2   (2D vector update, no Comm needed)
00079 
00080     // ==================================================================
00081     // Case 1 through 4 (A, B, C all local) Strided and non-strided cases
00082     // ==================================================================
00083 
00084     // Construct Vectors
00085 
00086   {
00087     Epetra_Vector A(LocalMap);
00088     Epetra_Vector B(LocalMap);
00089     Epetra_LocalMap  Map2d(NumVectors, IndexBase, Comm);
00090     Epetra_Vector C(Map2d);
00091     Epetra_Vector C_GEMM(Map2d);
00092 
00093     double **App, **Bpp, **Cpp;
00094     
00095     Epetra_Vector *Ap, *Bp, *Cp;
00096 
00097     // For testing non-strided mode, create Vectors that are scattered throughout memory
00098 
00099     App = new double *[NumVectors];
00100     Bpp = new double *[NumVectors];
00101     Cpp = new double *[NumVectors];
00102     for (i=0; i<NumVectors; i++) App[i] = new double[A.MyLength()+i];
00103     for (i=0; i<NumVectors; i++) Bpp[i] = new double[B.MyLength()+i];
00104     for (i=0; i<NumVectors; i++) Cpp[i] = new double[C.MyLength()+i];
00105     
00106     Epetra_Vector A1(View, LocalMap, App[0]);
00107     Epetra_Vector B1(View, LocalMap, Bpp[0]);
00108     Epetra_Vector C1(View, Map2d, Cpp[0]);
00109 
00110     for (int strided = 0; strided<2; strided++){
00111     // Loop through all trans cases using a variety of values for alpha and beta
00112     for (i=0; i<4; i++){
00113   int localierr = 0;
00114   char transa = 'N'; if (i>1) transa = 'T';
00115   char transb = 'N'; if (i%2!=0) transb = 'T';
00116   double alpha = (double) i+1;
00117   double beta  = (double) (i/2);
00118   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00119   localierr = BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00120         if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
00121     if (strided)
00122       {
00123         Ap = &A; Bp = &B; Cp = &C;
00124       }
00125     else
00126       {
00127         A.ExtractCopy(App[0]); Ap = &A1;
00128         B.ExtractCopy(Bpp[0]); Bp = &B1;
00129         C.ExtractCopy(Cpp[0]); Cp = &C1;
00130       }
00131     
00132     localierr = Cp->Multiply(transa, transb, alpha, *Ap, *Bp, beta);
00133     if (localierr!=-2) { // -2 means the shapes didn't match and we skip the tests
00134       ierr += Cp->Update(-1.0, C_GEMM, 1.0);
00135       ierr += Cp->Norm2(residual);
00136       
00137       if (verbose && ierr==0)
00138         {
00139     out << "XXXXX Replicated Local Vector GEMM tests";
00140     if (strided)
00141       out << " (Strided Multivectors)" << endl;
00142     else
00143       out << " (Non-Strided Multivectors)" << endl;
00144     out << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00145         <<", transb = " << transb;
00146         }
00147       if (ierr==0 && BadResidual(verbose,residual)) return(-1);
00148     }
00149   }
00150     }
00151     }
00152     for (i=0; i<NumVectors; i++)
00153       {
00154   delete [] App[i];
00155   delete [] Bpp[i];
00156   delete [] Cpp[i];
00157       }
00158     delete [] App;
00159     delete [] Bpp;
00160     delete [] Cpp;
00161   }
00162       
00163     // ====================================
00164     // Case 5  (A, B distributed C  local)
00165     // ====================================
00166 
00167     // Construct Vectors
00168   {
00169     Epetra_Vector A(Map);
00170     Epetra_Vector B(Map);
00171     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00172     Epetra_Vector C(Map2d);
00173     Epetra_Vector C_GEMM(Map2d);
00174 
00175     char transa = 'T';
00176     char transb = 'N';
00177     double alpha = 2.0;
00178     double beta  = 1.0;
00179     EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00180     ierr += BuildMatrixTests(C, transa, transb, alpha, A, B, beta, C_GEMM );
00181     ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00182     ierr += C.Update(-1.0, C_GEMM, 1.0);
00183     ierr += C.Norm2(residual);
00184 
00185     if (verbose && ierr==0)
00186       {
00187   out << "XXXXX Generalized 2D dot product via GEMM call     " << endl;
00188   out << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00189        <<", transb = " << transb;
00190       }
00191     if (BadResidual(verbose,residual)) return(-1);
00192     
00193     
00194   }      
00195     // ====================================
00196     // Case 6-7  (A, C distributed, B local)
00197     // ====================================
00198 
00199     // Construct Vectors
00200   {
00201     Epetra_Vector A(Map);
00202     Epetra_LocalMap Map2d(NumVectors, IndexBase, Comm);
00203     Epetra_Vector B(Map2d);
00204     Epetra_Vector C(Map);
00205     Epetra_Vector C_GEMM(Map);
00206 
00207     for (i=0; i<2; i++)
00208       {
00209   char transa = 'N';
00210   char transb = 'N'; if (i>0) transb = 'T';
00211   double alpha = 2.0;
00212   double beta  = 1.1;
00213   EPETRA_TEST_ERR(C.Random(),ierr);  // Fill C with random numbers
00214   ierr += BuildMatrixTests(C,transa, transb, alpha, A, B, beta, C_GEMM );
00215   ierr += C.Multiply(transa, transb, alpha, A, B, beta);
00216   ierr += C.Update(-1.0, C_GEMM, 1.0);
00217   ierr += C.Norm2(residual);
00218   
00219   if (verbose)
00220     {
00221       out << "XXXXX Generalized 2D vector update via GEMM call     " << endl;
00222       out << "  alpha = " << alpha << ",  beta = " << beta <<", transa = "<<transa
00223      <<", transb = " << transb;
00224     }
00225   if (BadResidual(verbose,residual)) return(-1);
00226       }
00227 
00228     delete [] residual;
00229     
00230     return(ierr);
00231   }
00232   }
00233 
00234 int VectorTests(const Epetra_BlockMap & Map, bool verbose)
00235 {
00236   int NumVectors = 1;
00237   int ierr = 0;
00238   double *residual = new double[NumVectors];
00239 
00240 #ifdef HAVE_EPETRA_TEUCHOS
00241   Teuchos::RCP<Teuchos::FancyOStream>
00242     fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
00243   std::ostream &out = *fancyOut;
00244 #else
00245   std::ostream &out = std::cout;
00246 #endif
00247   
00248   Epetra_BLAS BLAS;
00249   /* get number of processors and the name of this processor */
00250   
00251   // Construct Vectors
00252   
00253   Epetra_Vector A(Map);
00254   Epetra_Vector sqrtA(Map);
00255   Epetra_Vector B(Map);
00256   Epetra_Vector C(Map);
00257   Epetra_Vector C_alphaA(Map);
00258   Epetra_Vector C_alphaAplusB(Map);
00259   Epetra_Vector C_plusB(Map);
00260   Epetra_Vector Weights(Map);
00261   
00262   // Construct double vectors
00263   double *dotvec_AB   = new double[NumVectors];
00264   double *norm1_A     = new double[NumVectors];
00265   double *norm2_sqrtA = new double[NumVectors];
00266   double *norminf_A = new double[NumVectors];
00267   double *normw_A = new double[NumVectors];
00268   double *minval_A = new double[NumVectors];
00269   double *maxval_A = new double[NumVectors];
00270   double *meanval_A = new double[NumVectors];
00271 
00272   A.SetTracebackMode(1);
00273  
00274   // Generate data 
00275 
00276   
00277   EPETRA_TEST_ERR(C.Random(),ierr); // Fill C with random numbers.
00278   double alpha = 2.0;
00279   BuildVectorTests (C,alpha, A, sqrtA, B, C_alphaA, C_alphaAplusB,
00280            C_plusB, dotvec_AB, norm1_A, norm2_sqrtA, norminf_A, 
00281            normw_A, Weights, minval_A, maxval_A, meanval_A);
00282 
00283   int err = 0;
00284 
00285   // Test range checking for operator[](int)
00286 #ifdef HAVE_EPETRA_ARRAY_BOUNDS_CHECK
00287   try {
00288     if (verbose) out << "XXXXX Testing operator[A.MyLength()] bounds check     ";
00289     A[A.MyLength()]; // Should throw!
00290     // If we get here the test failed!
00291     EPETRA_TEST_ERR( 1, ierr );
00292   }
00293   catch ( const int& err_code ) {
00294     if (verbose) out << "\t Checked OK" << endl;
00295     EPETRA_TEST_ERR( ( err_code == -99 ? 0 : 1 ), ierr );
00296   }
00297   try {
00298     if (verbose) out << "XXXXX Testing operator[-1] bounds check     ";
00299     A[-1]; // Should throw!
00300     // If we get here the test failed!
00301     EPETRA_TEST_ERR( 1, ierr );
00302   }
00303   catch ( const int& err_code ) {
00304     if (verbose) out << "\t Checked OK" << endl;
00305     EPETRA_TEST_ERR( ( err_code == -99 ? 0 : 1 ), ierr );
00306   }
00307 #endif
00308 
00309   if (verbose) out << "XXXXX Testing alpha * A     ";
00310   // Test alpha*A
00311   Epetra_Vector alphaA(A); // Copy of A
00312   EPETRA_TEST_ERR(alphaA.Scale(alpha),err);
00313   EPETRA_TEST_ERR(alphaA.Update(-1.0, C_alphaA, 1.0),err);
00314   EPETRA_TEST_ERR(alphaA.Norm2(residual),err);
00315   
00316   if (err) ierr += err;
00317   else {
00318     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00319   }
00320 
00321   err = 0;
00322   if (verbose) out << "XXXXX Testing C = alpha * A + B      ";
00323   // Test alpha*A + B
00324   Epetra_Vector alphaAplusB(A); // Copy of A
00325   EPETRA_TEST_ERR(alphaAplusB.Update(1.0, B, alpha, A, 0.0),err);
00326   EPETRA_TEST_ERR(alphaAplusB.Update(-1.0, C_alphaAplusB, 1.0),err);
00327   EPETRA_TEST_ERR(alphaAplusB.Norm2(residual),err);
00328   
00329   if (err) ierr += err;
00330   else {
00331     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00332   }
00333 
00334   err = 0;
00335   if (verbose) out << "XXXXX Testing C += B      ";
00336   // Test + B
00337   Epetra_Vector plusB(C); // Copy of C
00338   EPETRA_TEST_ERR(plusB.Update(1.0, B, 1.0),err);
00339   EPETRA_TEST_ERR(plusB.Update(-1.0, C_plusB, 1.0),err);
00340   EPETRA_TEST_ERR(plusB.Norm2(residual),err);
00341   
00342   if (err) ierr += err;
00343   else {
00344     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00345   }
00346 
00347   err = 0;
00348   if (verbose) out << "XXXXX Testing A.dotProd(B)     ";
00349   // Test A.dotvec(B)
00350   double *dotvec = residual;
00351   EPETRA_TEST_ERR(A.Dot(B,dotvec),err);
00352   BLAS.AXPY(NumVectors,-1.0,dotvec_AB,dotvec);
00353   
00354   if (err) ierr += err;
00355   else {
00356   EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00357   }
00358   
00359   err = 0;
00360   if (verbose) out << "XXXXX Testing norm1_A      ";
00361   // Test A.norm1()
00362   double *norm1 = residual;
00363   EPETRA_TEST_ERR(A.Norm1(norm1),err);
00364   BLAS.AXPY(NumVectors,-1.0,norm1_A,norm1);
00365   
00366   if (err) ierr += err;
00367   else {
00368     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00369   }
00370     
00371   err = 0;
00372   if (verbose) out << "XXXXX Testing norm2_sqrtA     ";
00373   // Test sqrtA.norm2()
00374   double *norm2 = residual;
00375   EPETRA_TEST_ERR(sqrtA.Norm2(norm2),err);
00376   BLAS.AXPY(NumVectors,-1.0,norm2_sqrtA,norm2);
00377   
00378   if (err) ierr += err;
00379   else {
00380     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00381   }
00382   
00383   err = 0;
00384   if (verbose) out << "XXXXX Testing norminf_A     ";
00385   // Test A.norminf()
00386   double *norminf = residual;
00387   EPETRA_TEST_ERR(A.NormInf(norminf),ierr);
00388   BLAS.AXPY(NumVectors,-1.0,norminf_A,norminf);
00389   
00390   if (err) ierr += err;
00391   else {
00392     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00393   }
00394   
00395   err = 0;
00396   if (verbose) out << "XXXXX Testing normw_A     ";
00397   // Test A.NormWeighted()
00398   double *normw = residual;
00399   EPETRA_TEST_ERR(A.NormWeighted(Weights, normw),err);
00400   BLAS.AXPY(NumVectors,-1.0,normw_A,normw);
00401   
00402   if (err) ierr += err;
00403   else {
00404     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00405   }
00406   
00407   err = 0;
00408   if (verbose) out << "XXXXX Testing minval_A     ";
00409   // Test A.MinValue()
00410   double *minval = residual;
00411   EPETRA_TEST_ERR(A.MinValue(minval),err);
00412   BLAS.AXPY(NumVectors,-1.0,minval_A,minval);
00413   
00414   if (err) ierr += err;
00415   else {
00416     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00417   }
00418   
00419   err = 0;
00420   if (verbose) out << "XXXXX Testing maxval_A     ";
00421   // Test A.MaxValue()
00422   double *maxval = residual;
00423   EPETRA_TEST_ERR(A.MaxValue(maxval),err);
00424   BLAS.AXPY(NumVectors,-1.0,maxval_A,maxval);
00425   
00426   if (err) ierr += err;
00427   else {
00428     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00429   }
00430 
00431   err = 0;
00432   if (verbose) out << "XXXXX Testing meanval_A     ";
00433   // Test A.MeanValue()
00434   double *meanval = residual;
00435   EPETRA_TEST_ERR(A.MeanValue(meanval),err);
00436   BLAS.AXPY(NumVectors,-1.0,meanval_A,meanval);
00437   
00438   if (err) ierr += err;
00439   else {
00440     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00441   }
00442   
00443   err = 0;
00444   if (verbose) out << "XXXXX Testing abs_A     ";
00445   // Test A.Abs()
00446   Epetra_Vector Abs_A = A;
00447   EPETRA_TEST_ERR(Abs_A.Abs(A),err);
00448   EPETRA_TEST_ERR(Abs_A.Update(1.0, A, -1.0),err); // Abs_A = A - Abs_A (should be zero since A > 0)
00449   EPETRA_TEST_ERR(Abs_A.Norm2(residual),err);
00450 
00451   if (err) ierr += err;
00452   else {
00453     EPETRA_TEST_ERR(BadResidual(verbose,residual),ierr);
00454   } 
00455   
00456   // Delete everything
00457   
00458   delete [] dotvec_AB;
00459   delete [] norm1_A;
00460   delete [] norm2_sqrtA;
00461   delete [] norminf_A;
00462   delete [] normw_A;
00463   delete [] minval_A;
00464   delete [] maxval_A;
00465   delete [] meanval_A;
00466   delete [] residual;
00467   
00468   return(ierr);
00469 }
00470 
00471 int BadResidual(bool verbose, double * Residual)
00472 {
00473 
00474 #ifdef HAVE_EPETRA_TEUCHOS
00475   Teuchos::RCP<Teuchos::FancyOStream>
00476     fancyOut = Teuchos::VerboseObjectBase::getDefaultOStream();
00477   std::ostream &out = *fancyOut;
00478 #else
00479   std::ostream &out = std::cout;
00480 #endif
00481 
00482   double threshold = 5.0E-6;
00483   int ierr = 0;
00484     if (Residual[0]>threshold) {
00485       ierr = 1;// Output will be more useful after returning from method
00486       if (verbose) out << endl << "     Residual = " << Residual[0];
00487     }
00488   if (verbose)
00489     if (ierr==0) out << "\t Checked OK" << endl;
00490   
00491   return(ierr);
00492 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines