Epetra Package Browser (Single Doxygen Collection) Development
Vector/BuildTestProblems.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: Linear Algebra Services Package
00005 //                 Copyright 2011 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_Vector & C,
00046            const char TransA, const char TransB,
00047            const double alpha,
00048            Epetra_Vector& A,
00049            Epetra_Vector& B,
00050            const double beta,
00051            Epetra_Vector& C_GEMM ) {
00052 
00053     // For given values of TransA, TransB, alpha and beta, a (possibly
00054     // zero) filled Epetra_Vector C, and allocated
00055     // Epetra_Vectors A, B and C_GEMM this routine will generate values for
00056     // Epetra_Vectors 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 Vectors
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);
00089     B.ExtractView(&Bp);
00090     C.ExtractView(&Cp);
00091     C_GEMM.ExtractView(&C_GEMMp);
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  BuildVectorTests (Epetra_Vector & C, const double alpha,
00275         Epetra_Vector& A,
00276         Epetra_Vector& sqrtA,
00277         Epetra_Vector& B,
00278         Epetra_Vector& C_alphaA,
00279         Epetra_Vector& C_alphaAplusB,
00280         Epetra_Vector& 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_Vector& 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_Vector (the this object), allocated double * arguments dotvec_AB,
00293   // norm1_A, and norm2_A, and allocated Epetra_Vectors A, sqrtA,
00294   // B, C_alpha, C_alphaAplusB and C_plusB, this method will generate values for
00295   // Epetra_Vectors 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;  // 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     for (i = 0; i<A_nrows; i++) {
00396       fi = i+1; // Get float version of i and j, offset by 1.
00397       Ap[i] = (fi*sfinv);
00398       sqrtAp[i] = std::sqrt(Ap[i]);
00399     }
00400   }
00401   else {
00402     for (i = 0; i<A_nrows; i++) {
00403       fi = A_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00404       Ap[i] = (fi*sfinv);
00405       sqrtAp[i] = std::sqrt(Ap[i]);
00406     }
00407   }
00408 
00409   // Define B depending on TransB and B_is_local
00410 
00411   if (B_is_local) {
00412     for (i = 0; i<B_nrows; i++) {
00413       fi = i+1; // Get float version of i and j, offset by 1.
00414       Bp[i] = 1.0/((fi*sfinv));
00415       }
00416   }
00417   else {
00418     for (i = 0; i<B_nrows; i++) {
00419       fi = B_MyGlobalElements[i]+1; // Get float version of i and j, offset by 1.
00420       Bp[i] = 1.0/((fi*sfinv));
00421     }
00422   }
00423 
00424   // Generate C_alphaA = alpha * A
00425 
00426   for (i = 0; i<A_nrows; i++) C_alphaAp[i] = alpha * Ap[i];
00427 
00428   // Generate C_alphaA = alpha * A + B
00429 
00430     for (i = 0; i<A_nrows; i++) C_alphaAplusBp[i] = alpha * Ap[i] + Bp[i];
00431 
00432   // Generate C_plusB = this + B
00433 
00434     for (i = 0; i<A_nrows; i++) C_plusBp[i] = Cp[i] + Bp[i];
00435 
00436   // Generate dotvec_AB.  Because B(i,j) = 1/A(i,j), dotvec[i] =  C.GlobalLength()
00437 
00438   for (i=0; i< A.NumVectors(); i++) dotvec_AB[i] = C.GlobalLength();
00439 
00440   // For the next two results we want to be careful how we do arithmetic
00441   // to avoid very large numbers.
00442   // We are computing sfinv*(C.GlobalLength()*(C.GlobalLength()+1)/2)
00443 
00444       double result = C.GlobalLength();
00445       result *= sfinv;
00446       result /= 2.0;
00447       result *= (double)(C.GlobalLength()+1);
00448 
00449    // Generate norm1_A.  Can use formula for sum of first n integers.
00450 
00451   for (i=0; i< A.NumVectors(); i++)
00452     // m1_A[i] = (i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2;
00453     norm1_A[i] = result * ((double) (i+1));
00454 
00455   // Generate norm2_sqrtA.  Can use formula for sum of first n integers.
00456 
00457   for (i=0; i< A.NumVectors(); i++)
00458     // norm2_sqrtA[i] = std::sqrt((double) ((i+1)*C.GlobalLength()*(C.GlobalLength()+1)/2));
00459     norm2_sqrtA[i] = std::sqrt(result * ((double) (i+1)));
00460 
00461   // Generate norminf_A, minval_A, maxval_A, meanval_A.
00462 
00463   for (i=0; i< A.NumVectors(); i++)
00464     {
00465     norminf_A[i] = (double) (i+1);
00466     minval_A[i] =  (double) (i+1)/ (double) A.GlobalLength();
00467     maxval_A[i] = (double) (i+1);
00468     meanval_A[i] = norm1_A[i]/((double) (A.GlobalLength()));
00469     }
00470 
00471   // Define weights and expected weighted norm
00472       normw_A[0] = 1;
00473       for (j=0; j<A_nrows; j++) Weightsp[j] = Ap[j];
00474 
00475   delete A_Map;
00476   delete B_Map;
00477   delete [] A_MyGlobalElements;
00478   delete [] B_MyGlobalElements;
00479 
00480   return(0);
00481 }
00482 
00483 //========================================================================
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines