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

Generated on Wed May 12 21:41:06 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7