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