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