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

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