Vector/BuildTestProblems.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 "BuildTestProblems.h"
00030 
00031   int  BuildMatrixTests (Epetra_Vector & C,
00032            const char TransA, const char TransB, 
00033            const double alpha, 
00034            Epetra_Vector& A, 
00035            Epetra_Vector& B,
00036            const double beta,
00037            Epetra_Vector& C_GEMM ) {
00038 
00039     // For given values of TransA, TransB, alpha and beta, a (possibly
00040     // zero) filled Epetra_Vector C, and allocated
00041     // Epetra_Vectors A, B and C_GEMM this routine will generate values for 
00042     // Epetra_Vectors A, B and C_GEMM such that, if A, B and (this) are 
00043     // used with GEMM in this class, the results should match the results 
00044     // generated by this routine.
00045 
00046     // Test for Strided multivectors (required for GEMM ops)
00047 
00048     if (!A.ConstantStride()   ||
00049   !B.ConstantStride()      ||
00050   !C_GEMM.ConstantStride() ||
00051   !C.ConstantStride()) return(-1); // Error 
00052 
00053     int i, j;
00054     double fi, fj;  // Used for casting loop variables to floats
00055 
00056     // Get a view of the Vectors
00057 
00058     double *Ap      = 0;
00059     double *Bp      = 0;
00060     double *Cp      = 0;
00061     double *C_GEMMp = 0;
00062 
00063     int A_nrows = A.MyLength();
00064     int A_ncols = A.NumVectors();
00065     int B_nrows = B.MyLength();
00066     int B_ncols = B.NumVectors();
00067     int C_nrows = C.MyLength();
00068     int C_ncols = C.NumVectors();
00069     int A_Stride         = 0;
00070     int B_Stride         = 0;
00071     int C_Stride         = 0;
00072     int C_GEMM_Stride    = 0;
00073 
00074     A.ExtractView(&Ap);
00075     B.ExtractView(&Bp);
00076     C.ExtractView(&Cp);
00077     C_GEMM.ExtractView(&C_GEMMp);
00078 
00079       // Define some useful constants
00080 
00081 
00082     int opA_ncols = (TransA=='N') ? A.NumVectors() : A.MyLength();
00083     int opB_nrows = (TransB=='N') ? B.MyLength() : B.NumVectors();
00084   
00085     int C_global_inner_dim  = (TransA=='N') ? A.NumVectors() : A.GlobalLength();
00086   
00087 
00088     bool A_is_local = (!A.DistributedGlobal());
00089     bool B_is_local = (!B.DistributedGlobal());
00090     bool C_is_local = (!C.DistributedGlobal());
00091 
00092     int A_IndexBase = A.Map().IndexBase();
00093     int B_IndexBase = B.Map().IndexBase();
00094   
00095     // Build two new maps that we can use for defining global equation indices below
00096     Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
00097     Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());
00098 
00099     int* A_MyGlobalElements = new int[A_nrows];
00100     A_Map->MyGlobalElements(A_MyGlobalElements);
00101     int* B_MyGlobalElements = new int[B_nrows];
00102     B_Map->MyGlobalElements(B_MyGlobalElements);
00103 
00104   // Check for compatible dimensions
00105 
00106     if (C.MyLength()        != C_nrows     ||
00107   opA_ncols      != opB_nrows   ||
00108   C.NumVectors()    != C_ncols     ||
00109   C.MyLength()        != C_GEMM.MyLength()        ||
00110   C.NumVectors()    != C_GEMM.NumVectors()      ) {
00111       delete A_Map;
00112       delete B_Map;
00113       delete [] A_MyGlobalElements;
00114       delete [] B_MyGlobalElements;
00115       return(-2); // Return error
00116     }
00117 
00118     bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1 above
00119     bool Case2 = (!A_is_local && !B_is_local &&  C_is_local && TransA=='T' );// Case 2
00120     bool Case3 = (!A_is_local &&  B_is_local && !C_is_local && TransA=='N');// Case 3
00121   
00122     // Test for meaningful cases
00123 
00124     if (!(Case1 || Case2 || Case3)) {
00125       delete A_Map;
00126       delete B_Map;
00127       delete [] A_MyGlobalElements;
00128       delete [] B_MyGlobalElements;
00129       return(-3); // Meaningless case
00130     }
00131 
00132     /* Fill A, B and C with values as follows:
00133 
00134        If A_is_local is false:
00135        A(i,j) = A_MyGlobalElements[i]*j, i=1,...,numLocalEquations, j=1,...,NumVectors
00136        else
00137        A(i,j) = i*j,     i=1,...,numLocalEquations, j=1,...,NumVectors
00138        
00139        If B_is_local is false:
00140        B(i,j) = 1/(A_MyGlobalElements[i]*j), i=1,...,numLocalEquations, j=1,...,NumVectors
00141        
00142        else
00143        B(i,j) = 1/(i*j), i=1,...,numLocalEquations, j=1,...,NumVectors
00144        
00145        In addition, scale each entry by GlobalLength for A and
00146        1/GlobalLength for B--keeps the magnitude of entries in check
00147   
00148      
00149        C_GEMM will depend on A_is_local and B_is_local.  Three cases:
00150        
00151        1) A_is_local true and B_is_local true:
00152        C_GEMM will be local replicated and equal to A*B = i*NumVectors/j
00153        
00154        2) A_is_local false and B_is_local false
00155        C_GEMM will be local replicated = A(trans)*B(i,j) = i*numGlobalEquations/j
00156        
00157        3) A_is_local false B_is_local true
00158        C_GEMM will distributed global and equals A*B = A_MyGlobalElements[i]*NumVectors/j
00159        
00160     */
00161 
00162     // Define a scalar to keep magnitude of entries reasonable
00163 
00164     double sf = C_global_inner_dim;
00165     double sfinv = 1.0/sf;
00166 
00167     // Define A depending on A_is_local
00168 
00169     if (A_is_local)
00170       {
00171   for (j = 0; j <A_ncols ; j++) 
00172     for (i = 0; i<A_nrows; i++)
00173       { 
00174         fi = i+1; // Get float version of i and j, offset by 1.
00175         fj = j+1;
00176         Ap[i + A_Stride*j] = (fi*sfinv)*fj;
00177       }
00178       }
00179     else
00180       {
00181   for (j = 0; j <A_ncols ; j++) 
00182     for (i = 0; i<A_nrows; i++)
00183       { 
00184         fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00185         fj = j+1;
00186         Ap[i + A_Stride*j] = (fi*sfinv)*fj;
00187       }
00188       }
00189     
00190     // Define B depending on TransB and B_is_local
00191   
00192     if (B_is_local)
00193       {
00194   for (j = 0; j <B_ncols ; j++) 
00195     for (i = 0; i<B_nrows; i++)
00196       { 
00197         fi = i+1; // Get float version of i and j, offset by 1.
00198         fj = j+1;
00199         Bp[i + B_Stride*j] = 1.0/((fi*sfinv)*fj);
00200       }
00201       }
00202     else
00203       {
00204   for (j = 0; j <B_ncols ; j++) 
00205     for (i = 0; i<B_nrows; i++)
00206       { 
00207         fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00208         fj = j+1;
00209         Bp[i + B_Stride*j] = 1.0/((fi*sfinv)*fj);
00210       }
00211       }
00212     // Define C_GEMM depending on A_is_local and B_is_local.  C_GEMM is also a
00213     // function of alpha, beta, TransA, TransB: 
00214     
00215     //       C_GEMM = alpha*A(TransA)*B(TransB) + beta*C_GEMM
00216     
00217     if (Case1)
00218       {
00219   for (j = 0; j <C_ncols ; j++) 
00220     for (i = 0; i<C_nrows; i++)
00221       { 
00222         // Get float version of i and j, offset by 1.
00223         fi = (i+1)*C_global_inner_dim;
00224         fj = j+1;
00225         C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
00226     + beta * Cp[i + C_Stride*j];
00227       }
00228       }
00229     else if (Case2)
00230       {
00231   for (j = 0; j <C_ncols ; j++)
00232     for (i = 0; i<C_nrows; i++)
00233       { 
00234         // Get float version of i and j, offset by 1.
00235         fi = (i+1)*C_global_inner_dim;
00236         fj = j+1;
00237         C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
00238     + beta * Cp[i + C_Stride*j];
00239       }  
00240       }
00241     else
00242       {
00243   for (j = 0; j <C_ncols ; j++) 
00244     for (i = 0; i<C_nrows; i++)
00245       { 
00246         // Get float version of i and j.
00247         fi = (A_MyGlobalElements[i]+1)*C_global_inner_dim;
00248         fj = j+1;
00249         C_GEMMp[i + C_GEMM_Stride*j] = alpha * (fi/fj)
00250     + beta * Cp[i + C_Stride*j];
00251       }  
00252       }
00253     delete A_Map;
00254     delete B_Map;
00255     delete []A_MyGlobalElements;
00256     delete []B_MyGlobalElements;
00257 
00258     return(0);
00259   }
00260 int  BuildVectorTests (Epetra_Vector & C, const double alpha, 
00261         Epetra_Vector& A, 
00262         Epetra_Vector& sqrtA,
00263         Epetra_Vector& B,
00264         Epetra_Vector& C_alphaA,
00265         Epetra_Vector& C_alphaAplusB,
00266         Epetra_Vector& C_plusB,
00267         double* const dotvec_AB,
00268         double* const norm1_A,
00269         double* const norm2_sqrtA,
00270         double* const norminf_A,
00271         double* const normw_A,
00272         Epetra_Vector& Weights,
00273         double* const minval_A,
00274         double* const maxval_A,
00275         double* const meanval_A ) {
00276 
00277   // For given values alpha and a (possibly zero) filled 
00278   // Epetra_Vector (the this object), allocated double * arguments dotvec_AB, 
00279   // norm1_A, and norm2_A, and allocated Epetra_Vectors A, sqrtA,
00280   // B, C_alpha, C_alphaAplusB and C_plusB, this method will generate values for 
00281   // Epetra_Vectors A, B and all of the additional arguments on
00282   // the list above such that, if A, B and (this) are used with the methods in 
00283   // this class, the results should match the results generated by this routine.
00284   // Specifically, the results in dotvec_AB should match those from a call to 
00285   // A.dotProd (B,dotvec).  Similarly for other routines.
00286   
00287   int i,j;
00288   double fi;  // Used for casting loop variables to floats
00289   // Define some useful constants
00290   
00291   int A_nrows = A.MyLength();
00292   int A_ncols = A.NumVectors();
00293   int sqrtA_nrows = sqrtA.MyLength();
00294   int sqrtA_ncols = sqrtA.NumVectors();
00295   int B_nrows = B.MyLength();
00296   int B_ncols = B.NumVectors();
00297   
00298   double *Ap = 0;
00299   double *sqrtAp = 0;
00300   double *Bp = 0;
00301   double *Cp = 0;
00302   double *C_alphaAp = 0;
00303   double *C_alphaAplusBp = 0;
00304   double *C_plusBp = 0;
00305   double *Weightsp = 0;
00306 
00307   A.ExtractView(&Ap);
00308   sqrtA.ExtractView(&sqrtAp);
00309   B.ExtractView(&Bp);
00310   C.ExtractView(&Cp);
00311   C_alphaA.ExtractView(&C_alphaAp);
00312   C_alphaAplusB.ExtractView(&C_alphaAplusBp);
00313   C_plusB.ExtractView(&C_plusBp);
00314   Weights.ExtractView(&Weightsp);
00315   
00316   bool A_is_local = (A.MyLength() == A.GlobalLength());
00317   bool B_is_local = (B.MyLength() == B.GlobalLength());
00318   bool C_is_local = (C.MyLength()    == C.GlobalLength());
00319 
00320   int A_IndexBase = A.Map().IndexBase();
00321   int B_IndexBase = B.Map().IndexBase();
00322   
00323     // Build two new maps that we can use for defining global equation indices below
00324     Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
00325     Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());
00326 
00327     int* A_MyGlobalElements = new int[A_nrows];
00328     A_Map->MyGlobalElements(A_MyGlobalElements);
00329     int* B_MyGlobalElements = new int[B_nrows];
00330     B_Map->MyGlobalElements(B_MyGlobalElements);
00331 
00332   // Check for compatible dimensions
00333   
00334   if (C.MyLength()        != A_nrows     ||
00335       A_nrows        != B_nrows     ||
00336       C.NumVectors()    != A_ncols     ||
00337       A_ncols        != B_ncols     ||
00338       sqrtA_nrows    != A_nrows     ||
00339       sqrtA_ncols    != A_ncols     ||
00340       C.MyLength()        != C_alphaA.MyLength()     ||
00341       C.NumVectors()    != C_alphaA.NumVectors() ||
00342       C.MyLength()        != C_alphaAplusB.MyLength()     ||
00343       C.NumVectors()    != C_alphaAplusB.NumVectors() ||
00344       C.MyLength()        != C_plusB.MyLength()      ||
00345       C.NumVectors()    != C_plusB.NumVectors()     ) return(-2); // Return error
00346   
00347  
00348   bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1
00349   bool Case2 = (!A_is_local && !B_is_local && !C_is_local);// Case 2
00350   
00351   // Test for meaningful cases
00352   
00353   if (!(Case1 || Case2)) return(-3); // Meaningless case
00354   
00355   /* Fill A and B with values as follows:
00356      
00357      If A_is_local is false:
00358      A(i,j) = A_MyGlobalElements[i]*j,     i=1,...,numLocalEquations, j=1,...,NumVectors
00359      
00360      else
00361      A(i,j) = i*j,     i=1,...,numLocalEquations, j=1,...,NumVectors
00362 
00363      If B_is_local is false:
00364      B(i,j) = 1/(A_MyGlobalElements[i]*j), i=1,...,numLocalEquations, j=1,...,NumVectors
00365 
00366      else
00367      B(i,j) = 1/(i*j), i=1,...,numLocalEquations, j=1,...,NumVectors
00368 
00369      In addition, scale each entry by GlobalLength for A and
00370      1/GlobalLength for B--keeps the magnitude of entries in check
00371   */
00372 
00373   //Define scale factor
00374 
00375   double sf = A.GlobalLength();
00376   double sfinv = 1.0/sf;
00377 
00378   // Define A
00379 
00380   if (A_is_local) {
00381     for (i = 0; i<A_nrows; i++) { 
00382       fi = i+1; // Get float version of i and j, offset by 1.
00383       Ap[i] = (fi*sfinv);
00384       sqrtAp[i] = std::sqrt(Ap[i]);
00385     }
00386   }
00387   else {
00388     for (i = 0; i<A_nrows; i++) { 
00389       fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00390       Ap[i] = (fi*sfinv);
00391       sqrtAp[i] = std::sqrt(Ap[i]);
00392     }
00393   }
00394 
00395   // Define B depending on TransB and B_is_local
00396   
00397   if (B_is_local) {
00398     for (i = 0; i<B_nrows; i++) { 
00399       fi = i+1; // Get float version of i and j, offset by 1.
00400       Bp[i] = 1.0/((fi*sfinv));
00401       }
00402   }
00403   else {
00404     for (i = 0; i<B_nrows; i++) { 
00405       fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00406       Bp[i] = 1.0/((fi*sfinv));
00407     }
00408   }
00409   
00410   // Generate C_alphaA = alpha * A
00411 
00412   for (i = 0; i<A_nrows; i++) C_alphaAp[i] = alpha * Ap[i];
00413 
00414   // Generate C_alphaA = alpha * A + B
00415 
00416     for (i = 0; i<A_nrows; i++) C_alphaAplusBp[i] = alpha * Ap[i] + Bp[i];
00417   
00418   // Generate C_plusB = this + B
00419 
00420     for (i = 0; i<A_nrows; i++) C_plusBp[i] = Cp[i] + Bp[i];
00421 
00422   // Generate dotvec_AB.  Because B(i,j) = 1/A(i,j), dotvec[i] =  C.GlobalLength()
00423 
00424   for (i=0; i< A.NumVectors(); i++) dotvec_AB[i] = C.GlobalLength();
00425 
00426   // For the next two results we want to be careful how we do arithmetic 
00427   // to avoid very large numbers.
00428   // We are computing sfinv*(C.GlobalLength()*(C.GlobalLength()+1)/2)
00429 
00430       double result = C.GlobalLength();
00431       result *= sfinv;
00432       result /= 2.0;
00433       result *= (double)(C.GlobalLength()+1);
00434 
00435    // Generate norm1_A.  Can use formula for sum of first n integers.
00436 
00437   for (i=0; i< A.NumVectors(); i++) 
00438     // m1_A[i] = (i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2;
00439     norm1_A[i] = result * ((double) (i+1));
00440 
00441   // Generate norm2_sqrtA.  Can use formula for sum of first n integers. 
00442 
00443   for (i=0; i< A.NumVectors(); i++) 
00444     // norm2_sqrtA[i] = std::sqrt((double) ((i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2));
00445     norm2_sqrtA[i] = std::sqrt(result * ((double) (i+1)));
00446 
00447   // Generate norminf_A, minval_A, maxval_A, meanval_A. 
00448 
00449   for (i=0; i< A.NumVectors(); i++) 
00450     {
00451     norminf_A[i] = (double) (i+1);
00452     minval_A[i] =  (double) (i+1)/ (double) A.GlobalLength();
00453     maxval_A[i] = (double) (i+1);
00454     meanval_A[i] = norm1_A[i]/((double) (A.GlobalLength()));
00455     }
00456 
00457   // Define weights and expected weighted norm
00458       normw_A[0] = 1;
00459       for (j=0; j<A_nrows; j++) Weightsp[j] = Ap[j];
00460 
00461   delete A_Map;
00462   delete B_Map;
00463   delete [] A_MyGlobalElements;
00464   delete [] B_MyGlobalElements;
00465   
00466   return(0);
00467 }
00468 
00469 //========================================================================

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