Teuchos Package Browser (Single Doxygen Collection) Version of the Day
cxx_tmpl_main.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) 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 // 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 // Kris
00043 // 07.24.03 -- Initial checkin
00044 // 08.08.03 -- All test suites except for TRSM are finished.
00045 // 08.14.03 -- The test suite for TRSM is finished (Heidi).
00046 /*
00047 
00048 This test program is intended to check an experimental default type (e.g. mp_real) against an "officialy supported" control type (e.g. double). For each test, the program will generate the appropriate scalars and randomly-sized vectors and matrices of random numbers for the routine being tested. All of the input data for the experimental type is casted into the control type, so in theory both BLAS routines should get the same input data. Upon return, the experimental output data is casted back into the control type, and the results are compared; if they are equal (within a user-definable tolerance) the test is considered successful.
00049 
00050 The test routine for TRSM is still being developed; all of the others are more or less finalized.
00051 
00052 */
00053 
00054 #include "Teuchos_BLAS.hpp"
00055 #include "Teuchos_Time.hpp"
00056 #include "Teuchos_Version.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using Teuchos::BLAS;
00060 using Teuchos::ScalarTraits;
00061 
00062 // SType1 and SType2 define the datatypes for which BLAS output will be compared.
00063 // SType2 should generally be a control datatype "officially" supported by the BLAS; SType1 should be the experimental type being checked.
00064 
00065 // Define the first scalar type
00066 #ifdef HAVE_TEUCHOS_COMPLEX
00067 #define SType1     std::complex<float>
00068 #else
00069 #define SType1     float
00070 #endif
00071 
00072 // Define the second scalar type
00073 #ifdef HAVE_TEUCHOS_COMPLEX
00074 #define SType2     std::complex<double>
00075 #else
00076 #define SType2     double
00077 #endif
00078 
00079 // Define the ordinal type
00080 #define OType    long int 
00081 
00082 // MVMIN/MAX define the minimum and maximum dimensions of generated matrices and vectors, respectively.
00083 #define MVMIN      2
00084 #define MVMAX      20
00085 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and std::vector elements and scalars:
00086 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
00087 // Set SCALARMAX to a floating-point value (e.g. 10.0) to enable floating-point random number generation, such that
00088 // random numbers in (-SCALARMAX - 1, SCALARMAX + 1) will be generated.
00089 // Large values of SCALARMAX may cause problems with SType2 = int, as large integer values will overflow floating-point types.
00090 #define SCALARMAX  10
00091 // These define the number of tests to be run for each individual BLAS routine.
00092 #define ROTGTESTS  5
00093 #define ROTTESTS   5
00094 #define ASUMTESTS  5
00095 #define AXPYTESTS  5
00096 #define COPYTESTS  5
00097 #define DOTTESTS   5
00098 #define IAMAXTESTS 5
00099 #define NRM2TESTS  5
00100 #define SCALTESTS  5
00101 #define GEMVTESTS  5
00102 #define GERTESTS   0
00103 #define TRMVTESTS  5
00104 #define GEMMTESTS  5
00105 #define SYMMTESTS  5
00106 #define SYRKTESTS  5
00107 #define TRMMTESTS  5
00108 #define TRSMTESTS  5
00109 
00110 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
00111 template<typename TYPE>
00112 TYPE GetRandom(TYPE, TYPE);
00113 
00114 // Returns a random integer between the two input parameters, inclusive
00115 template<>
00116 int GetRandom(int, int);
00117 
00118 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
00119 template<>
00120 double GetRandom(double, double);
00121 
00122 template<typename TYPE>
00123 void PrintVector(TYPE* Vector, OType Size, std::string Name, bool Matlab = 0);
00124 
00125 template<typename TYPE>
00126 void PrintMatrix(TYPE* Matrix, OType Rows, OType Columns, OType LDM, std::string Name, bool Matlab = 0);
00127 
00128 template<typename TYPE1, typename TYPE2>
00129 bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, typename ScalarTraits<TYPE2>::magnitudeType Tolerance ); 
00130 
00131 template<typename TYPE1, typename TYPE2>
00132 bool CompareVectors(TYPE1* Vector1, TYPE2* Vector2, OType Size, typename ScalarTraits<TYPE2>::magnitudeType Tolerance ); 
00133 
00134 template<typename TYPE1, typename TYPE2>
00135 bool CompareMatrices(TYPE1* Matrix1, TYPE2* Matrix2, OType Rows, OType Columns, OType LDM, typename ScalarTraits<TYPE2>::magnitudeType Tolerance ); 
00136 
00137 // For most types, this function is just a wrapper for static_cast(), but for mp_real/double, it calls mp::dble()
00138 // The second input parameter is not used; it is only needed to determine what type to convert *to*
00139 template<typename TYPE1, typename TYPE2>
00140 TYPE2 ConvertType(TYPE1, TYPE2);
00141 
00142 // These functions return a random character appropriate for the BLAS arguments that share their names (uses GetRandom())
00143 Teuchos::ESide RandomSIDE();
00144 Teuchos::EUplo RandomUPLO();
00145 Teuchos::ETransp RandomTRANS();
00146 Teuchos::EDiag RandomDIAG();
00147 
00148 int main(int argc, char *argv[])
00149 {
00150   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00151   bool verbose = 0;
00152   bool debug = 0;
00153   bool matlab = 0;
00154   bool InvalidCmdLineArgs = 0;
00155   OType i, j, k;
00156   for(i = 1; i < argc; i++)
00157     {
00158       if(argv[i][0] == '-')
00159   {
00160     switch(argv[i][1])
00161       {
00162       case 'v':
00163         if(!verbose)
00164     {
00165       verbose = 1;
00166     }
00167         else
00168     {
00169       InvalidCmdLineArgs = 1;
00170     }
00171         break;
00172       case 'd':
00173         if(!debug)
00174     {
00175       debug = 1;
00176     }
00177         else
00178     {
00179       InvalidCmdLineArgs = 1;
00180     }
00181         break;
00182       case 'm':
00183         if(!matlab)
00184     {
00185       matlab = 1;
00186     }
00187         else
00188     {
00189       InvalidCmdLineArgs = 1;
00190     }
00191         break;
00192       default:
00193         InvalidCmdLineArgs = 1;
00194         break;
00195       }
00196   }
00197     }
00198 
00199   if (verbose)
00200     std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
00201 
00202   if(InvalidCmdLineArgs || (argc > 4))
00203     {
00204       std::cout << "Invalid command line arguments detected. Use the following flags:" << std::endl
00205      << "\t -v enables verbose mode (reports number of failed/successful tests)" << std::endl
00206      << "\t -d enables debug mode (same as verbose with output of each test, not recommended for large numbers of tests)" << std::endl
00207      << "\t -m enables matlab-style output; only has an effect if debug mode is enabled" << std::endl;
00208       return 1;
00209     }
00210   typedef ScalarTraits<SType1>::magnitudeType MType1;
00211   typedef ScalarTraits<SType2>::magnitudeType MType2; 
00212   BLAS<OType, SType1> SType1BLAS;
00213   BLAS<OType, SType2> SType2BLAS;
00214   SType1 SType1zero = ScalarTraits<SType1>::zero();
00215   SType1 SType1one = ScalarTraits<SType1>::one();
00216   SType2 SType2one = ScalarTraits<SType2>::one();
00217   SType1* SType1A;
00218   SType1* SType1B;
00219   SType1* SType1C;
00220   SType1* SType1x;
00221   SType1* SType1y;
00222   SType1 SType1alpha, SType1beta;
00223   SType2* SType2A;
00224   SType2* SType2B;
00225   SType2* SType2C;
00226   SType2* SType2x;
00227   SType2* SType2y; 
00228   SType2 SType2alpha, SType2beta;
00229   SType1 SType1ASUMresult, SType1DOTresult, SType1NRM2result, SType1SINresult;
00230   SType2 SType2ASUMresult, SType2DOTresult, SType2NRM2result, SType2SINresult;
00231   MType1 SType1COSresult;
00232   MType2 SType2COSresult;
00233   OType incx, incy;
00234   OType SType1IAMAXresult;
00235   OType SType2IAMAXresult;
00236   OType TotalTestCount = 1, GoodTestSubcount, GoodTestCount = 0, M, M2, N, N2, P, K, LDA, LDB, LDC, Mx, My;
00237   Teuchos::EUplo UPLO;
00238   Teuchos::ESide SIDE;
00239   Teuchos::ETransp TRANS, TRANSA, TRANSB;
00240   Teuchos::EDiag DIAG;
00241   SType2 convertTo = ScalarTraits<SType2>::zero();
00242   MType2 mConvertTo = ScalarTraits<MType2>::zero();
00243   MType2 TOL = 1e-5*ScalarTraits<MType2>::one();
00244   
00245   std::srand(time(NULL));
00246 
00247   //--------------------------------------------------------------------------------
00248   // BEGIN LEVEL 1 BLAS TESTS
00249   //--------------------------------------------------------------------------------
00250   // Begin ROTG Tests
00251   //--------------------------------------------------------------------------------
00252   GoodTestSubcount = 0;
00253   for(i = 0; i < ROTGTESTS; i++)
00254     {
00255       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
00256       SType2alpha = ConvertType(SType1alpha, convertTo);
00257       SType1beta = GetRandom(-SCALARMAX, SCALARMAX);
00258       SType2beta = ConvertType(SType1beta, convertTo);
00259       SType1COSresult = ScalarTraits<MType1>::zero();
00260       SType2COSresult = ConvertType(SType1COSresult, ScalarTraits<MType2>::zero());
00261       SType1SINresult = ScalarTraits<SType1>::zero();
00262       SType2SINresult = ConvertType(SType1SINresult, convertTo);
00263       
00264       if(debug)
00265   {
00266     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00267     std::cout << "SType1alpha = "  << SType1alpha << std::endl;
00268     std::cout << "SType2alpha = " << SType2alpha << std::endl;
00269     std::cout << "SType1beta = "  << SType1beta << std::endl;
00270     std::cout << "SType2beta = " << SType2beta << std::endl;
00271   }
00272       TotalTestCount++;
00273       SType1BLAS.ROTG(&SType1alpha, &SType1beta, &SType1COSresult, &SType1SINresult);
00274       SType2BLAS.ROTG(&SType2alpha, &SType2beta, &SType2COSresult, &SType2SINresult);
00275       if(debug)
00276   {
00277     std::cout << "SType1 ROTG COS result: " << SType1COSresult << std::endl;
00278     std::cout << "SType2 ROTG COS result: " << SType2COSresult << std::endl;
00279     std::cout << "SType1 ROTG SIN result: " << SType1SINresult << std::endl;
00280     std::cout << "SType2 ROTG SIN result: " << SType2SINresult << std::endl;
00281   }
00282       GoodTestSubcount += ( CompareScalars(SType1COSresult, SType2COSresult, TOL) && 
00283           CompareScalars(SType1SINresult, SType2SINresult, TOL) );
00284     }
00285   GoodTestCount += GoodTestSubcount;
00286   if(verbose || debug) std::cout << "ROTG: " << GoodTestSubcount << " of " << ROTGTESTS << " tests were successful." << std::endl;
00287   if(debug) std::cout << std::endl;
00288   //--------------------------------------------------------------------------------
00289   // End ROTG Tests
00290   //--------------------------------------------------------------------------------
00291 
00292   //--------------------------------------------------------------------------------
00293   // Begin ROT Tests
00294   //--------------------------------------------------------------------------------
00295   typedef Teuchos::ScalarTraits<SType1>::magnitudeType MType1;
00296   typedef Teuchos::ScalarTraits<SType2>::magnitudeType MType2;
00297   GoodTestSubcount = 0;
00298   for(i = 0; i < ROTTESTS; i++)
00299   {
00300     incx = GetRandom(-5,5);
00301     incy = GetRandom(-5,5);
00302     if (incx == 0) incx = 1;
00303     if (incy == 0) incy = 1;
00304     M = GetRandom(MVMIN, MVMIN+8);
00305     Mx = M*std::abs(incx);
00306     My = M*std::abs(incy);
00307     if (Mx == 0) { Mx = 1; }
00308     if (My == 0) { My = 1; }
00309     SType1x = new SType1[Mx];
00310     SType1y = new SType1[My];
00311     SType2x = new SType2[Mx];
00312     SType2y = new SType2[My];
00313     for(j = 0; j < Mx; j++)
00314     {
00315       SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00316       SType2x[j] = ConvertType(SType1x[j], convertTo);
00317     }
00318     for(j = 0; j < My; j++)
00319     {
00320       SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00321       SType2y[j] = ConvertType(SType1y[j], convertTo);
00322     }
00323     MType1 c1 = Teuchos::ScalarTraits<SType1>::magnitude(cos(static_cast<double>(GetRandom(-SCALARMAX,SCALARMAX))));
00324     MType2 c2 = ConvertType(c1, mConvertTo);
00325     SType1 s1 = sin(static_cast<double>(GetRandom(-SCALARMAX,SCALARMAX)));
00326     SType2 s2 = ConvertType(s1, convertTo);
00327     if(debug)
00328     {
00329       std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00330       std::cout << "c1 = "  << c1 << ", s1 = " << s1 << std::endl;
00331       std::cout << "c2 = " << c2 << ", s2 = " << s2 << std::endl;
00332       std::cout << "incx = " << incx << ", incy = " << incy << std::endl;
00333       PrintVector(SType1x, Mx, "SType1x", matlab);
00334       PrintVector(SType1y, My, "SType1y_before_operation", matlab);
00335       PrintVector(SType2x, Mx, "SType2x", matlab);
00336       PrintVector(SType2y, My, "SType2y_before_operation",  matlab);
00337     }
00338     TotalTestCount++;
00339     SType1BLAS.ROT(M, SType1x, incx, SType1y, incy, &c1, &s1);
00340     SType2BLAS.ROT(M, SType2x, incx, SType2y, incy, &c2, &s2);
00341     if(debug)
00342     {
00343       PrintVector(SType1y, My, "SType1y_after_operation", matlab);
00344       PrintVector(SType2y, My, "SType2y_after_operation", matlab);
00345     }
00346     GoodTestSubcount += ( CompareVectors(SType1x, SType2x, Mx, TOL) &&
00347                           CompareVectors(SType1y, SType2y, My, TOL) );
00348     delete [] SType1x;
00349     delete [] SType1y;
00350     delete [] SType2x;
00351     delete [] SType2y;
00352   }
00353   GoodTestCount += GoodTestSubcount;
00354   if(verbose || debug) std::cout << "ROT: " << GoodTestSubcount << " of " << ROTTESTS << " tests were successful." << std::endl;
00355   if(debug) std::cout << std::endl;
00356   //--------------------------------------------------------------------------------
00357   // End ROT Tests
00358   //--------------------------------------------------------------------------------
00359 
00360   //--------------------------------------------------------------------------------
00361   // Begin ASUM Tests
00362   //--------------------------------------------------------------------------------
00363   GoodTestSubcount = 0;
00364   ScalarTraits<int>::seedrandom(0);
00365   for(i = 0; i < ASUMTESTS; i++)
00366     {
00367       incx = GetRandom(1, SCALARMAX);
00368       M = GetRandom(MVMIN, MVMAX);
00369       M2 = M*incx;
00370       SType1x = new SType1[M2];
00371       SType2x = new SType2[M2];
00372       for(j = 0; j < M2; j++)
00373   {
00374     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00375     SType2x[j] = ConvertType(SType1x[j], convertTo);
00376   }
00377       if(debug)
00378   {
00379     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00380     PrintVector(SType1x, M2, "SType1x", matlab);
00381     PrintVector(SType2x, M2, "SType2x", matlab);
00382   }
00383       TotalTestCount++;
00384       SType1ASUMresult = SType1BLAS.ASUM(M, SType1x, incx);
00385       SType2ASUMresult = SType2BLAS.ASUM(M, SType2x, incx);
00386       if(debug)
00387   {
00388     std::cout << "SType1 ASUM result: " << SType1ASUMresult << std::endl;
00389     std::cout << "SType2 ASUM result: " << SType2ASUMresult << std::endl;
00390   }
00391       GoodTestSubcount += CompareScalars(SType1ASUMresult, SType2ASUMresult, TOL);
00392       delete [] SType1x;
00393       delete [] SType2x;
00394     }
00395   GoodTestCount += GoodTestSubcount;
00396   if(verbose || debug) std::cout << "ASUM: " << GoodTestSubcount << " of " << ASUMTESTS << " tests were successful." << std::endl;
00397   if(debug) std::cout << std::endl;
00398 
00399   //--------------------------------------------------------------------------------
00400   // End ASUM Tests
00401   //--------------------------------------------------------------------------------
00402 
00403   //--------------------------------------------------------------------------------
00404   // Begin AXPY Tests
00405   //--------------------------------------------------------------------------------
00406   GoodTestSubcount = 0;
00407   for(i = 0; i < AXPYTESTS; i++)
00408     {
00409       incx = GetRandom(1, MVMAX);
00410       incy = GetRandom(1, MVMAX);
00411       M = GetRandom(MVMIN, MVMAX);
00412       Mx = M*std::abs(incx);
00413       My = M*std::abs(incy);
00414       if (Mx == 0) { Mx = 1; }
00415       if (My == 0) { My = 1; }
00416       SType1x = new SType1[Mx];
00417       SType1y = new SType1[My];
00418       SType2x = new SType2[Mx];
00419       SType2y = new SType2[My]; 
00420       for(j = 0; j < Mx; j++)
00421   {
00422     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00423     SType2x[j] = ConvertType(SType1x[j], convertTo);
00424   }
00425       for(j = 0; j < My; j++)
00426   {
00427     SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00428     SType2y[j] = ConvertType(SType1y[j], convertTo);
00429   }
00430       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
00431       SType2alpha = ConvertType(SType1alpha, convertTo);
00432       if(debug)
00433   {
00434     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00435     std::cout << "SType1alpha = "  << SType1alpha << std::endl;
00436     std::cout << "SType2alpha = " << SType2alpha << std::endl;
00437     PrintVector(SType1x, Mx, "SType1x", matlab);
00438     PrintVector(SType1y, My, "SType1y_before_operation", matlab);
00439     PrintVector(SType2x, Mx, "SType2x", matlab);
00440     PrintVector(SType2y, My, "SType2y_before_operation",  matlab);
00441   }
00442       TotalTestCount++;
00443       SType1BLAS.AXPY(M, SType1alpha, SType1x, incx, SType1y, incy);
00444       SType2BLAS.AXPY(M, SType2alpha, SType2x, incx, SType2y, incy);
00445       if(debug)
00446   {
00447     PrintVector(SType1y, My, "SType1y_after_operation", matlab);
00448     PrintVector(SType2y, My, "SType2y_after_operation", matlab);
00449   }
00450       GoodTestSubcount += CompareVectors(SType1y, SType2y, My, TOL);
00451       delete [] SType1x;
00452       delete [] SType1y;
00453       delete [] SType2x;
00454       delete [] SType2y;
00455     }
00456   GoodTestCount += GoodTestSubcount;
00457   if(verbose || debug) std::cout << "AXPY: " << GoodTestSubcount << " of " << AXPYTESTS << " tests were successful." << std::endl;
00458   if(debug) std::cout << std::endl;
00459   //--------------------------------------------------------------------------------
00460   // End AXPY Tests
00461   //--------------------------------------------------------------------------------
00462 
00463   //--------------------------------------------------------------------------------
00464   // Begin COPY Tests
00465   //--------------------------------------------------------------------------------
00466   GoodTestSubcount = 0;
00467   for(i = 0; i < COPYTESTS; i++)
00468     {
00469       incx = GetRandom(1, MVMAX);
00470       incy = GetRandom(1, MVMAX);
00471       M = GetRandom(MVMIN, MVMAX);
00472       Mx = M*std::abs(incx);
00473       My = M*std::abs(incy);
00474       if (Mx == 0) { Mx = 1; }
00475       if (My == 0) { My = 1; }
00476       SType1x = new SType1[Mx];
00477       SType1y = new SType1[My];
00478       SType2x = new SType2[Mx];
00479       SType2y = new SType2[My]; 
00480       for(j = 0; j < Mx; j++)
00481   {
00482     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00483     SType2x[j] = ConvertType(SType1x[j], convertTo);
00484   }
00485       for(j = 0; j < My; j++)
00486   {
00487     SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00488     SType2y[j] = ConvertType(SType1y[j], convertTo);
00489   }
00490       if(debug)
00491   {
00492     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00493     PrintVector(SType1x, Mx, "SType1x", matlab);
00494     PrintVector(SType1y, My, "SType1y_before_operation", matlab);
00495     PrintVector(SType2x, Mx, "SType2x", matlab);
00496     PrintVector(SType2y, My, "SType2y_before_operation", matlab);
00497   }
00498       TotalTestCount++;
00499       SType1BLAS.COPY(M, SType1x, incx, SType1y, incy);
00500       SType2BLAS.COPY(M, SType2x, incx, SType2y, incy);
00501       if(debug)
00502   {
00503     PrintVector(SType1y, My, "SType1y_after_operation", matlab);
00504     PrintVector(SType2y, My, "SType2y_after_operation", matlab);
00505   }
00506       GoodTestSubcount += CompareVectors(SType1y, SType2y, My, TOL);
00507       delete [] SType1x;
00508       delete [] SType1y;
00509       delete [] SType2x;
00510       delete [] SType2y;
00511     }
00512    GoodTestCount += GoodTestSubcount; if(verbose || debug) std::cout << "COPY: " << GoodTestSubcount << " of " << COPYTESTS << " tests were successful." << std::endl;
00513   if(debug) std::cout << std::endl;
00514   //--------------------------------------------------------------------------------
00515   // End COPY Tests
00516   //--------------------------------------------------------------------------------
00517 
00518   //--------------------------------------------------------------------------------
00519   // Begin DOT Tests
00520   //--------------------------------------------------------------------------------
00521   GoodTestSubcount = 0;
00522   for(i = 0; i < DOTTESTS; i++)
00523     {
00524       incx = GetRandom(1, MVMAX);
00525       incy = GetRandom(1, MVMAX);
00526       M = GetRandom(MVMIN, MVMAX);
00527       Mx = M*std::abs(incx);
00528       My = M*std::abs(incy);
00529       if (Mx == 0) { Mx = 1; }
00530       if (My == 0) { My = 1; }
00531       SType1x = new SType1[Mx];
00532       SType1y = new SType1[My];
00533       SType2x = new SType2[Mx];
00534       SType2y = new SType2[My]; 
00535       for(j = 0; j < Mx; j++)
00536   {
00537     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00538     SType2x[j] = ConvertType(SType1x[j], convertTo);
00539   }
00540       for(j = 0; j < My; j++)
00541   {
00542     SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00543     SType2y[j] = ConvertType(SType1y[j], convertTo);
00544   }
00545       if(debug)
00546   {
00547     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00548     PrintVector(SType1x, Mx, "SType1x", matlab);
00549     PrintVector(SType1y, My, "SType1y", matlab);
00550     PrintVector(SType2x, Mx, "SType2x", matlab);
00551     PrintVector(SType2y, My, "SType2y", matlab);
00552   }
00553       TotalTestCount++;
00554       SType1DOTresult = SType1BLAS.DOT(M, SType1x, incx, SType1y, incy);
00555       SType2DOTresult = SType2BLAS.DOT(M, SType2x, incx, SType2y, incy);
00556       if(debug)
00557   {
00558     std::cout << "SType1 DOT result: " << SType1DOTresult << std::endl;
00559     std::cout << "SType2 DOT result: " << SType2DOTresult << std::endl;
00560   }
00561       GoodTestSubcount += CompareScalars(SType1DOTresult, SType2DOTresult, TOL);
00562       delete [] SType1x;
00563       delete [] SType1y;
00564       delete [] SType2x;
00565       delete [] SType2y;
00566     }
00567   GoodTestCount += GoodTestSubcount;
00568   if(verbose || debug) std::cout << "DOT: " << GoodTestSubcount << " of " << DOTTESTS << " tests were successful." << std::endl;
00569   if(debug) std::cout << std::endl;
00570   //--------------------------------------------------------------------------------
00571   // End DOT Tests
00572   //--------------------------------------------------------------------------------
00573 
00574   //--------------------------------------------------------------------------------
00575   // Begin NRM2 Tests
00576   //--------------------------------------------------------------------------------
00577   GoodTestSubcount = 0;
00578   for(i = 0; i < NRM2TESTS; i++)
00579     {
00580       incx = GetRandom(1, SCALARMAX);
00581       M = GetRandom(MVMIN, MVMAX);
00582       M2 = M*incx; 
00583       SType1x = new SType1[M2];
00584       SType2x = new SType2[M2];
00585       for(j = 0; j < M2; j++)
00586   {
00587     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00588     SType2x[j] = ConvertType(SType1x[j], convertTo);
00589   }
00590       if(debug)
00591   {
00592     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00593     PrintVector(SType1x, M2, "SType1x", matlab);
00594     PrintVector(SType2x, M2, "SType2x", matlab);
00595   }
00596       TotalTestCount++;
00597       SType1NRM2result = SType1BLAS.NRM2(M, SType1x, incx);
00598       SType2NRM2result = SType2BLAS.NRM2(M, SType2x, incx);
00599       if(debug)
00600   {
00601     std::cout << "SType1 NRM2 result: " << SType1NRM2result << std::endl;
00602     std::cout << "SType2 NRM2 result: " << SType2NRM2result << std::endl;
00603   }
00604       GoodTestSubcount += CompareScalars(SType1NRM2result, SType2NRM2result, TOL);
00605       delete [] SType1x;
00606       delete [] SType2x;
00607     }
00608    GoodTestCount += GoodTestSubcount; if(verbose || debug) std::cout << "NRM2: " << GoodTestSubcount << " of " << NRM2TESTS << " tests were successful." << std::endl;
00609   if(debug) std::cout << std::endl;
00610   //--------------------------------------------------------------------------------
00611   // End NRM2 Tests
00612   //--------------------------------------------------------------------------------
00613 
00614   //--------------------------------------------------------------------------------
00615   // Begin SCAL Tests
00616   //--------------------------------------------------------------------------------
00617   GoodTestSubcount = 0;
00618   for(i = 0; i < SCALTESTS; i++)
00619     {
00620       // These will only test for the case that the increment is > 0, the
00621       // templated case can handle when incx < 0, but the blas library doesn't 
00622       // seem to be able to on some machines.
00623       incx = GetRandom(1, SCALARMAX);
00624       M = GetRandom(MVMIN, MVMAX);
00625       M2 = M*incx;
00626       SType1x = new SType1[M2];
00627       SType2x = new SType2[M2];
00628       for(j = 0; j < M2; j++)
00629   {
00630     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00631     SType2x[j] = ConvertType(SType1x[j], convertTo);
00632   }
00633       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
00634       SType2alpha = ConvertType(SType1alpha, convertTo);
00635       if(debug)
00636   {
00637     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00638     std::cout << "SType1alpha = " << SType1alpha << std::endl;
00639     std::cout << "SType2alpha = " << SType2alpha << std::endl;
00640     PrintVector(SType1x, M2, "SType1x_before_operation", matlab);
00641     PrintVector(SType2x, M2, "SType2x_before_operation", matlab);
00642   }
00643       TotalTestCount++;
00644       SType1BLAS.SCAL(M, SType1alpha, SType1x, incx);
00645       SType2BLAS.SCAL(M, SType2alpha, SType2x, incx);
00646       if(debug)
00647   {
00648     PrintVector(SType1x, M2, "SType1x_after_operation", matlab);
00649     PrintVector(SType2x, M2, "SType2x_after_operation", matlab);
00650   }
00651       GoodTestSubcount += CompareVectors(SType1x, SType2x, M2, TOL);
00652       delete [] SType1x;
00653       delete [] SType2x;
00654     }
00655   GoodTestCount += GoodTestSubcount;
00656   if(verbose || debug) std::cout << "SCAL: " << GoodTestSubcount << " of " << SCALTESTS << " tests were successful." << std::endl;
00657   if(debug) std::cout << std::endl;
00658   //--------------------------------------------------------------------------------
00659   // End SCAL Tests
00660   //--------------------------------------------------------------------------------
00661 
00662   //--------------------------------------------------------------------------------
00663   // Begin IAMAX Tests
00664   //--------------------------------------------------------------------------------
00665   GoodTestSubcount = 0;
00666   for(i = 0; i < IAMAXTESTS; i++)
00667     {
00668       incx = GetRandom(1, SCALARMAX);
00669       M = GetRandom(MVMIN, MVMAX);
00670       M2 = M*incx;
00671       SType1x = new SType1[M2];
00672       SType2x = new SType2[M2];
00673       for(j = 0; j < M2; j++)
00674   {
00675     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00676     SType2x[j] = ConvertType(SType1x[j], convertTo);
00677   }
00678       if(debug)
00679   {
00680     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00681     PrintVector(SType1x, M2, "SType1x", matlab);
00682     PrintVector(SType2x, M2, "SType2x", matlab);
00683   }
00684       TotalTestCount++;
00685       SType1IAMAXresult = SType1BLAS.IAMAX(M, SType1x, incx);
00686       SType2IAMAXresult = SType2BLAS.IAMAX(M, SType2x, incx);
00687       if(debug)
00688   {
00689     std::cout << "SType1 IAMAX result: " << SType1IAMAXresult << std::endl;
00690     std::cout << "SType2 IAMAX result: " << SType2IAMAXresult << std::endl;
00691   }
00692       GoodTestSubcount += (SType1IAMAXresult == SType2IAMAXresult);
00693       delete [] SType1x;
00694       delete [] SType2x;
00695     }
00696   GoodTestCount += GoodTestSubcount;
00697   if(verbose || debug) std::cout << "IAMAX: " << GoodTestSubcount << " of " << IAMAXTESTS << " tests were successful." << std::endl;
00698   if(debug) std::cout << std::endl;
00699   //--------------------------------------------------------------------------------
00700   // End IAMAX Tests
00701   //--------------------------------------------------------------------------------
00702 
00703   //--------------------------------------------------------------------------------
00704   // BEGIN LEVEL 2 BLAS TESTS
00705   //--------------------------------------------------------------------------------
00706   // Begin GEMV Tests
00707   //--------------------------------------------------------------------------------
00708   GoodTestSubcount = 0;
00709   for(i = 0; i < GEMVTESTS; i++)
00710     {
00711       // The parameters used to construct the test problem are chosen to be
00712       // valid parameters, so the GEMV routine won't bomb out.
00713       incx = GetRandom(1, MVMAX);
00714       while (incx == 0) {
00715           incx = GetRandom(1, MVMAX);
00716       }
00717       incy = GetRandom(1, MVMAX);
00718       while (incy == 0) {
00719           incy = GetRandom(1, MVMAX);
00720       }   
00721       M = GetRandom(MVMIN, MVMAX);
00722       N = GetRandom(MVMIN, MVMAX);
00723 
00724       TRANS = RandomTRANS();
00725       if (Teuchos::ETranspChar[TRANS] == 'N') { 
00726         M2 = M*std::abs(incy);
00727         N2 = N*std::abs(incx);   
00728       } else {
00729   M2 = N*std::abs(incy);
00730   N2 = M*std::abs(incx);
00731       }
00732 
00733       LDA = GetRandom(MVMIN, MVMAX);
00734       while (LDA < M) {
00735           LDA = GetRandom(MVMIN, MVMAX);
00736       }   
00737 
00738       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
00739       SType1beta = GetRandom(-SCALARMAX, SCALARMAX);
00740       SType2alpha = ConvertType(SType1alpha, convertTo);
00741       SType2beta = ConvertType(SType1beta, convertTo);
00742 
00743       SType1A = new SType1[LDA * N];
00744       SType1x = new SType1[N2];
00745       SType1y = new SType1[M2];
00746       SType2A = new SType2[LDA * N];
00747       SType2x = new SType2[N2];
00748       SType2y = new SType2[M2]; 
00749 
00750       for(j = 0; j < LDA * N; j++)
00751   {
00752     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
00753     SType2A[j] = ConvertType(SType1A[j], convertTo);
00754   }
00755       for(j = 0; j < N2; j++)
00756   {
00757     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00758     SType2x[j] = ConvertType(SType1x[j], convertTo);
00759   }
00760       for(j = 0; j < M2; j++)
00761   {
00762     SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00763     SType2y[j] = ConvertType(SType1y[j], convertTo);
00764   }
00765       if(debug)
00766   {
00767     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00768     std::cout << "SType1alpha = " << SType1alpha << std::endl;
00769     std::cout << "SType2alpha = " << SType2alpha << std::endl;
00770     std::cout << "SType1beta = " << SType1beta << std::endl;
00771     std::cout << "SType2beta = " << SType2beta << std::endl;
00772     PrintMatrix(SType1A, M, N, LDA, "SType1A", matlab);
00773     PrintVector(SType1x, N2, "SType1x", matlab);
00774     PrintVector(SType1y, M2, "SType1y_before_operation", matlab);
00775     PrintMatrix(SType2A, M, N, LDA, "SType2A", matlab);
00776     PrintVector(SType2x, N2, "SType2x", matlab);
00777     PrintVector(SType2y, M2, "SType2y_before_operation", matlab);
00778   }
00779       TotalTestCount++;
00780       SType1BLAS.GEMV(TRANS, M, N, SType1alpha, SType1A, LDA, SType1x, incx, SType1beta, SType1y, incy);
00781       SType2BLAS.GEMV(TRANS, M, N, SType2alpha, SType2A, LDA, SType2x, incx, SType2beta, SType2y, incy);
00782       if(debug)
00783   {
00784     PrintVector(SType1y, M2, "SType1y_after_operation", matlab);
00785     PrintVector(SType2y, M2, "SType2y_after_operation", matlab);
00786   }
00787       GoodTestSubcount += CompareVectors(SType1y, SType2y, M2, TOL);
00788       delete [] SType1A;
00789       delete [] SType1x;
00790       delete [] SType1y;
00791       delete [] SType2A;
00792       delete [] SType2x;
00793       delete [] SType2y;
00794     }
00795   GoodTestCount += GoodTestSubcount;
00796   if(verbose || debug) std::cout << "GEMV: " << GoodTestSubcount << " of " << GEMVTESTS << " tests were successful." << std::endl;
00797   if(debug) std::cout << std::endl;
00798   //--------------------------------------------------------------------------------
00799   // End GEMV Tests
00800   //--------------------------------------------------------------------------------
00801 
00802   //--------------------------------------------------------------------------------
00803   // Begin TRMV Tests
00804   //--------------------------------------------------------------------------------
00805   GoodTestSubcount = 0;
00806   for(i = 0; i < TRMVTESTS; i++)
00807     {
00808       UPLO = RandomUPLO();
00809       TRANSA = RandomTRANS();
00810       
00811       // Since the entries are integers, we don't want to use the unit diagonal feature,
00812       // this creates ill-conditioned, nearly-singular matrices.
00813       //DIAG = RandomDIAG();  
00814       DIAG = Teuchos::NON_UNIT_DIAG;
00815 
00816       N = GetRandom(MVMIN, MVMAX);
00817       incx = GetRandom(1, MVMAX);
00818       while (incx == 0) {
00819           incx = GetRandom(1, MVMAX);
00820       }
00821       N2 = N*std::abs(incx);
00822       SType1x = new SType1[N2];
00823       SType2x = new SType2[N2];
00824 
00825       for(j = 0; j < N2; j++)
00826   {
00827     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00828     SType2x[j] = ConvertType(SType1x[j], convertTo);
00829   }
00830 
00831       LDA = GetRandom(MVMIN, MVMAX);
00832       while (LDA < N) {
00833   LDA = GetRandom(MVMIN, MVMAX);
00834       }
00835       SType1A = new SType1[LDA * N];
00836       SType2A = new SType2[LDA * N];
00837 
00838       for(j = 0; j < N; j++)
00839   {      
00840     if(Teuchos::EUploChar[UPLO] == 'U') {
00841       // The operator is upper triangular, make sure that the entries are
00842       // only in the upper triangular part of A and the diagonal is non-zero.
00843       for(k = 0; k < N; k++) 
00844       {
00845         if(k < j) {
00846     SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00847         } else {
00848     SType1A[j*LDA+k] = SType1zero;
00849         }
00850         SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
00851         if(k == j) {
00852     if (Teuchos::EDiagChar[DIAG] == 'N') {
00853       SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00854       while (SType1A[j*LDA+k] == SType1zero) {
00855         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00856       }
00857       SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
00858     } else {
00859       SType1A[j*LDA+k] = SType1one;
00860       SType2A[j*LDA+k] = SType2one;
00861     }
00862         }     
00863       }
00864     } else {
00865       // The operator is lower triangular, make sure that the entries are
00866       // only in the lower triangular part of A and the diagonal is non-zero.
00867       for(k = 0; k < N; k++) 
00868         {
00869     if(k > j) {
00870       SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00871     } else {
00872       SType1A[j*LDA+k] = SType1zero;
00873     }
00874     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
00875     if(k == j) {
00876       if (Teuchos::EDiagChar[DIAG] == 'N') {
00877         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00878         while (SType1A[j*LDA+k] == SType1zero) {
00879           SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
00880         }
00881         SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
00882       } else {
00883         SType1A[j*LDA+k] = SType1one;
00884         SType2A[j*LDA+k] = SType2one;
00885       }
00886     }     
00887         } // end for(k=0 ...    
00888     } // end if(UPLO == 'U') ...
00889   } // end for(j=0 ...      for(j = 0; j < N*N; j++)
00890       
00891       if(debug)
00892   {
00893     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00894     PrintMatrix(SType1A, N, N, LDA,"SType1A", matlab);
00895     PrintVector(SType1x, N2, "SType1x_before_operation", matlab);
00896     PrintMatrix(SType2A, N, N, LDA, "SType2A", matlab);
00897     PrintVector(SType2x, N2, "SType2x_before_operation", matlab);
00898   }
00899       TotalTestCount++;
00900       SType1BLAS.TRMV(UPLO, TRANSA, DIAG, N, SType1A, LDA, SType1x, incx);
00901       SType2BLAS.TRMV(UPLO, TRANSA, DIAG, N, SType2A, LDA, SType2x, incx);
00902       if(debug)
00903   {
00904     PrintVector(SType1x, N2, "SType1x_after_operation", matlab);
00905     PrintVector(SType2x, N2, "SType2x_after_operation", matlab);
00906   }
00907       GoodTestSubcount += CompareVectors(SType1x, SType2x, N2, TOL);
00908       delete [] SType1A;
00909       delete [] SType1x;
00910       delete [] SType2A;
00911       delete [] SType2x;
00912     }
00913   GoodTestCount += GoodTestSubcount;
00914   if(verbose || debug) std::cout << "TRMV: " << GoodTestSubcount << " of " << TRMVTESTS << " tests were successful." << std::endl;
00915   if(debug) std::cout << std::endl;
00916   //--------------------------------------------------------------------------------
00917   // End TRMV Tests
00918   //--------------------------------------------------------------------------------
00919 
00920   //--------------------------------------------------------------------------------
00921   // Begin GER Tests
00922   //--------------------------------------------------------------------------------
00923   GoodTestSubcount = 0;
00924   for(i = 0; i < GERTESTS; i++)
00925     {
00926       incx = GetRandom(1, MVMAX);
00927       while (incx == 0) {
00928           incx = GetRandom(1, MVMAX);
00929       }   
00930       incy = GetRandom(1, MVMAX);
00931       while (incy == 0) {
00932           incy = GetRandom(1, MVMAX);
00933       }   
00934       M = GetRandom(MVMIN, MVMAX);
00935       N = GetRandom(MVMIN, MVMAX);
00936 
00937       M2 = M*std::abs(incx);
00938       N2 = N*std::abs(incy);   
00939 
00940       LDA = GetRandom(MVMIN, MVMAX);
00941       while (LDA < M) {
00942           LDA = GetRandom(MVMIN, MVMAX);
00943       }   
00944 
00945       SType1A = new SType1[LDA * N];
00946       SType1x = new SType1[M2];
00947       SType1y = new SType1[N2];
00948       SType2A = new SType2[LDA * N];
00949       SType2x = new SType2[M2];
00950       SType2y = new SType2[N2];
00951       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
00952       SType2alpha = ConvertType(SType1alpha, convertTo);
00953       for(j = 0; j < LDA * N; j++)
00954   {
00955     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
00956     SType2A[j] = ConvertType(SType1A[j], convertTo);
00957   }
00958       for(j = 0; j < M2; j++)
00959   {
00960     SType1x[j] = GetRandom(-SCALARMAX, SCALARMAX);
00961     SType2x[j] = ConvertType(SType1x[j], convertTo);
00962   }
00963       for(j = 0; j < N2; j++)
00964   {
00965     SType1y[j] = GetRandom(-SCALARMAX, SCALARMAX);
00966     SType2y[j] = ConvertType(SType1y[j], convertTo);
00967   }
00968       if(debug)
00969   {
00970     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
00971     std::cout << "SType1alpha = " << SType1alpha << std::endl;
00972     std::cout << "SType2alpha = " << SType2alpha << std::endl;
00973     PrintMatrix(SType1A, M, N, LDA,"SType1A_before_operation", matlab);
00974     PrintVector(SType1x, M2, "SType1x", matlab);
00975     PrintVector(SType1y, N2, "SType1y", matlab);
00976     PrintMatrix(SType2A, M, N, LDA,"SType2A_before_operation", matlab);
00977     PrintVector(SType2x, M2, "SType2x", matlab);
00978     PrintVector(SType2y, N2, "SType2y", matlab);
00979   }
00980       TotalTestCount++;
00981       SType1BLAS.GER(M, N, SType1alpha, SType1x, incx, SType1y, incy, SType1A, LDA);
00982       SType2BLAS.GER(M, N, SType2alpha, SType2x, incx, SType2y, incy, SType2A, LDA);
00983       if(debug)
00984   {
00985     PrintMatrix(SType1A, M, N, LDA, "SType1A_after_operation", matlab);
00986     PrintMatrix(SType2A, M, N, LDA, "SType2A_after_operation", matlab);
00987   }
00988       GoodTestSubcount += CompareMatrices(SType1A, SType2A, M, N, LDA, TOL);
00989       delete [] SType1A;
00990       delete [] SType1x;
00991       delete [] SType1y;
00992       delete [] SType2A;
00993       delete [] SType2x;
00994       delete [] SType2y;
00995     }
00996   GoodTestCount += GoodTestSubcount;
00997   if(verbose || debug) std::cout << "GER: " << GoodTestSubcount << " of " << GERTESTS << " tests were successful." << std::endl;
00998   if(debug) std::cout << std::endl;
00999   //--------------------------------------------------------------------------------
01000   // End GER Tests
01001   //--------------------------------------------------------------------------------
01002 
01003   //--------------------------------------------------------------------------------
01004   // BEGIN LEVEL 3 BLAS TESTS
01005   //--------------------------------------------------------------------------------
01006   // Begin GEMM Tests
01007   //--------------------------------------------------------------------------------
01008   GoodTestSubcount = 0;
01009   for(i = 0; i < GEMMTESTS; i++)
01010     { 
01011       TRANSA = RandomTRANS();      
01012       TRANSB = RandomTRANS();
01013       M = GetRandom(MVMIN, MVMAX);
01014       N = GetRandom(MVMIN, MVMAX);
01015       P = GetRandom(MVMIN, MVMAX);
01016 
01017       if(debug) {
01018     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
01019       }
01020       LDA = GetRandom(MVMIN, MVMAX);
01021       if (Teuchos::ETranspChar[TRANSA] == 'N') {
01022   while (LDA < M) {  LDA = GetRandom(MVMIN, MVMAX); }
01023   SType1A = new SType1[LDA * P];
01024   SType2A = new SType2[LDA * P];
01025   for(j = 0; j < LDA * P; j++)
01026   {
01027     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01028     SType2A[j] = ConvertType(SType1A[j], convertTo);
01029   }
01030   if (debug) {
01031   PrintMatrix(SType1A, M, P, LDA, "SType1A", matlab);
01032   PrintMatrix(SType2A, M, P, LDA, "SType2A", matlab);
01033   }
01034       } else {
01035   while (LDA < P) {  LDA = GetRandom(MVMIN, MVMAX); }
01036   SType1A = new SType1[LDA * M];
01037   SType2A = new SType2[LDA * M];
01038   for(j = 0; j < LDA * M; j++)
01039   {
01040     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01041     SType2A[j] = ConvertType(SType1A[j], convertTo);
01042   }
01043   if (debug) {
01044   PrintMatrix(SType1A, P, M, LDA, "SType1A", matlab);
01045   PrintMatrix(SType2A, P, M, LDA, "SType2A", matlab);
01046   }
01047       }
01048 
01049       LDB = GetRandom(MVMIN, MVMAX);
01050       if (Teuchos::ETranspChar[TRANSB] == 'N') {
01051   while (LDB < P) {  LDB = GetRandom(MVMIN, MVMAX); }
01052   SType1B = new SType1[LDB * N];
01053   SType2B = new SType2[LDB * N];
01054   for(j = 0; j < LDB * N; j++)
01055   {
01056     SType1B[j] = GetRandom(-SCALARMAX, SCALARMAX);
01057     SType2B[j] = ConvertType(SType1B[j], convertTo);
01058   }
01059   if (debug) {
01060     PrintMatrix(SType1B, P, N, LDB,"SType1B", matlab);
01061     PrintMatrix(SType2B, P, N, LDB,"SType2B", matlab);
01062   }
01063       } else { 
01064   while (LDB < N) {  LDB = GetRandom(MVMIN, MVMAX); }
01065   SType1B = new SType1[LDB * P];
01066   SType2B = new SType2[LDB * P];
01067   for(j = 0; j < LDB * P; j++)
01068   {
01069     SType1B[j] = GetRandom(-SCALARMAX, SCALARMAX);
01070     SType2B[j] = ConvertType(SType1B[j], convertTo);
01071   } 
01072   if (debug) {
01073     PrintMatrix(SType1B, N, P, LDB,"SType1B", matlab);
01074     PrintMatrix(SType2B, N, P, LDB,"SType2B", matlab);
01075   }
01076       }
01077 
01078       LDC = GetRandom(MVMIN, MVMAX);
01079       while (LDC < M) {  LDC = GetRandom(MVMIN, MVMAX); }
01080       SType1C = new SType1[LDC * N];
01081       SType2C = new SType2[LDC * N];
01082       for(j = 0; j < LDC * N; j++) {
01083     SType1C[j] = GetRandom(-SCALARMAX, SCALARMAX);
01084     SType2C[j] = ConvertType(SType1C[j], convertTo);
01085       }
01086       if(debug)
01087   {
01088     PrintMatrix(SType1C, M, N, LDC, "SType1C_before_operation", matlab);
01089     PrintMatrix(SType2C, M, N, LDC, "SType2C_before_operation", matlab);
01090   }
01091   
01092       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
01093       SType1beta = GetRandom(-SCALARMAX, SCALARMAX);
01094       SType2alpha = ConvertType(SType1alpha, convertTo);
01095       SType2beta = ConvertType(SType1beta, convertTo);
01096 
01097       TotalTestCount++;
01098       SType1BLAS.GEMM(TRANSA, TRANSB, M, N, P, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
01099       SType2BLAS.GEMM(TRANSA, TRANSB, M, N, P, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
01100       if(debug)
01101   {
01102     PrintMatrix(SType1C, M, N, LDC, "SType1C_after_operation", matlab);
01103     PrintMatrix(SType2C, M, N, LDC, "SType2C_after_operation", matlab);
01104   }
01105       GoodTestSubcount += CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
01106       delete [] SType1A;
01107       delete [] SType1B;
01108       delete [] SType1C;
01109       delete [] SType2A;
01110       delete [] SType2B;
01111       delete [] SType2C;
01112     }
01113   GoodTestCount += GoodTestSubcount;
01114   if(verbose || debug) std::cout << "GEMM: " << GoodTestSubcount << " of " << GEMMTESTS << " tests were successful." << std::endl;
01115   if(debug) std::cout << std::endl;
01116   //--------------------------------------------------------------------------------
01117   // End GEMM Tests
01118   //--------------------------------------------------------------------------------
01119 
01120   //--------------------------------------------------------------------------------
01121   // Begin SYMM Tests
01122   //--------------------------------------------------------------------------------
01123   GoodTestSubcount = 0;
01124   for(i = 0; i < SYMMTESTS; i++)
01125     { 
01126       M = GetRandom(MVMIN, MVMAX);
01127       N = GetRandom(MVMIN, MVMAX);
01128       SIDE = RandomSIDE();
01129       UPLO = RandomUPLO();
01130 
01131       LDA = GetRandom(MVMIN, MVMAX);
01132       if(Teuchos::ESideChar[SIDE] == 'L') {
01133   while (LDA < M) { LDA = GetRandom(MVMIN, MVMAX); }
01134   SType1A = new SType1[LDA * M];
01135   SType2A = new SType2[LDA * M];
01136   for(j = 0; j < LDA * M; j++) {
01137     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01138     SType2A[j] = ConvertType(SType1A[j], convertTo);
01139   }
01140       } else {
01141   while (LDA < N) { LDA = GetRandom(MVMIN, MVMAX); }
01142   SType1A = new SType1[LDA * N];
01143   SType2A = new SType2[LDA * N];
01144   for(j = 0; j < LDA * N; j++) {
01145     SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01146     SType2A[j] = ConvertType(SType1A[j], convertTo);
01147   }
01148       }
01149 
01150       LDB = GetRandom(MVMIN, MVMAX);
01151       while (LDB < M) {  LDB = GetRandom(MVMIN, MVMAX); }
01152       SType1B = new SType1[LDB * N];
01153       SType2B = new SType2[LDB * N];
01154       for(j = 0; j < LDB * N; j++) {
01155     SType1B[j] = GetRandom(-SCALARMAX, SCALARMAX);
01156     SType2B[j] = ConvertType(SType1B[j], convertTo);
01157       }
01158     
01159       LDC = GetRandom(MVMIN, MVMAX);
01160       while (LDC < M) {  LDC = GetRandom(MVMIN, MVMAX); }
01161       SType1C = new SType1[LDC * N];
01162       SType2C = new SType2[LDC * N];
01163       for(j = 0; j < LDC * N; j++) {
01164     SType1C[j] = GetRandom(-SCALARMAX, SCALARMAX);
01165     SType2C[j] = ConvertType(SType1C[j], convertTo);
01166       }
01167       
01168       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
01169       SType1beta = GetRandom(-SCALARMAX, SCALARMAX);
01170       SType2alpha = ConvertType(SType1alpha, convertTo);
01171       SType2beta = ConvertType(SType1beta, convertTo);
01172       if(debug)
01173   {
01174     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
01175     std::cout << "SType1alpha = " << SType1alpha << std::endl;
01176     std::cout << "SType2alpha = " << SType2alpha << std::endl;
01177     std::cout << "SType1beta = " << SType1beta << std::endl;
01178     std::cout << "SType2beta = " << SType2beta << std::endl;
01179     if (Teuchos::ESideChar[SIDE] == 'L') {
01180         PrintMatrix(SType1A, M, M, LDA,"SType1A", matlab);
01181         PrintMatrix(SType2A, M, M, LDA,"SType2A", matlab);
01182     } else {
01183       PrintMatrix(SType1A, N, N, LDA, "SType1A", matlab);
01184       PrintMatrix(SType2A, N, N, LDA, "SType2A", matlab);
01185     }
01186     PrintMatrix(SType1B, M, N, LDB,"SType1B", matlab);
01187     PrintMatrix(SType1C, M, N, LDC,"SType1C_before_operation", matlab);
01188     PrintMatrix(SType2B, M, N, LDB,"SType2B", matlab);
01189     PrintMatrix(SType2C, M, N, LDC,"SType2C_before_operation", matlab);
01190   }
01191       TotalTestCount++;
01192 
01193       SType1BLAS.SYMM(SIDE, UPLO, M, N, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
01194       SType2BLAS.SYMM(SIDE, UPLO, M, N, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
01195       if(debug)
01196   {
01197     PrintMatrix(SType1C, M, N, LDC,"SType1C_after_operation", matlab);
01198     PrintMatrix(SType2C, M, N, LDC,"SType2C_after_operation", matlab);
01199   }
01200       GoodTestSubcount += CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
01201 
01202       delete [] SType1A;
01203       delete [] SType1B;
01204       delete [] SType1C;
01205       delete [] SType2A;
01206       delete [] SType2B;
01207       delete [] SType2C;
01208     }
01209   GoodTestCount += GoodTestSubcount;
01210   if(verbose || debug) std::cout << "SYMM: " << GoodTestSubcount << " of " << SYMMTESTS << " tests were successful." << std::endl;
01211   if(debug) std::cout << std::endl;
01212   //--------------------------------------------------------------------------------
01213   // End SYMM Tests
01214   //--------------------------------------------------------------------------------
01215 
01216   //--------------------------------------------------------------------------------
01217   // Begin SYRK Tests
01218   //--------------------------------------------------------------------------------
01219   GoodTestSubcount = 0;
01220   for(i = 0; i < SYRKTESTS; i++)
01221   {
01222     N = GetRandom(MVMIN, MVMAX);
01223     K = GetRandom(MVMIN, MVMAX);
01224     while (K > N) { K = GetRandom(MVMIN, MVMAX); }
01225 
01226     UPLO = RandomUPLO();
01227     TRANS = RandomTRANS();
01228 #ifdef HAVE_TEUCHOS_COMPLEX
01229     while (TRANS == Teuchos::CONJ_TRANS) { TRANS = RandomTRANS(); }
01230 #endif
01231 
01232     LDA = GetRandom(MVMIN, MVMAX);
01233     if(Teuchos::ETranspChar[TRANS] == 'N') {
01234       while (LDA < N) { LDA = GetRandom(MVMIN, MVMAX); }
01235       SType1A = new SType1[LDA * K];
01236       SType2A = new SType2[LDA * K];
01237       for(j = 0; j < LDA * K; j++) {
01238         SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01239         SType2A[j] = ConvertType(SType1A[j], convertTo);
01240       }
01241     } else {
01242       while (LDA < K) { LDA = GetRandom(MVMIN, MVMAX); }
01243       SType1A = new SType1[LDA * N];
01244       SType2A = new SType2[LDA * N];
01245       for(j = 0; j < LDA * N; j++) {
01246         SType1A[j] = GetRandom(-SCALARMAX, SCALARMAX);
01247         SType2A[j] = ConvertType(SType1A[j], convertTo);
01248       }
01249     }
01250 
01251     LDC = GetRandom(MVMIN, MVMAX);
01252     while (LDC < N) {  LDC = GetRandom(MVMIN, MVMAX); }
01253     SType1C = new SType1[LDC * N];
01254     SType2C = new SType2[LDC * N];
01255     for(j = 0; j < N; j++) {
01256 
01257       if(Teuchos::EUploChar[UPLO] == 'U') {
01258         // The operator is upper triangular, make sure that the entries are
01259         // only in the upper triangular part of C.
01260         for(k = 0; k < N; k++)
01261         {
01262           if(k <= j) {
01263             SType1C[j*LDC+k] = GetRandom(-SCALARMAX, SCALARMAX);
01264           } else {
01265             SType1C[j*LDC+k] = SType1zero;
01266           } 
01267           SType2C[j*LDC+k] = ConvertType(SType1C[j*LDC+k], convertTo);
01268         }
01269       }
01270       else {
01271         for(k = 0; k < N; k++)
01272         {
01273           if(k >= j) {
01274             SType1C[j*LDC+k] = GetRandom(-SCALARMAX, SCALARMAX);
01275           } else {
01276             SType1C[j*LDC+k] = SType1zero;
01277           }
01278           SType2C[j*LDC+k] = ConvertType(SType1C[j*LDC+k], convertTo);
01279         }
01280       } 
01281     }
01282        
01283     SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
01284     SType1beta = GetRandom(-SCALARMAX, SCALARMAX);
01285     SType2alpha = ConvertType(SType1alpha, convertTo);
01286     SType2beta = ConvertType(SType1beta, convertTo);
01287     if(debug)
01288     {
01289       std::cout << "Test #" << TotalTestCount << " --" << std::endl;
01290       std::cout << "SType1alpha = " << SType1alpha << std::endl;
01291       std::cout << "SType2alpha = " << SType2alpha << std::endl;
01292       std::cout << "SType1beta = " << SType1beta << std::endl;
01293       std::cout << "SType2beta = " << SType2beta << std::endl;
01294       if (Teuchos::ETranspChar[TRANS] == 'N') {
01295         PrintMatrix(SType1A, N, K, LDA,"SType1A", matlab);
01296         PrintMatrix(SType2A, N, K, LDA,"SType2A", matlab);
01297       } else {
01298         PrintMatrix(SType1A, K, N, LDA, "SType1A", matlab);
01299         PrintMatrix(SType2A, K, N, LDA, "SType2A", matlab);
01300       }
01301       PrintMatrix(SType1C, N, N, LDC,"SType1C_before_operation", matlab);
01302       PrintMatrix(SType2C, N, N, LDC,"SType2C_before_operation", matlab);
01303     }
01304     TotalTestCount++;
01305 
01306     SType1BLAS.SYRK(UPLO, TRANS, N, K, SType1alpha, SType1A, LDA, SType1beta, SType1C, LDC);
01307     SType2BLAS.SYRK(UPLO, TRANS, N, K, SType2alpha, SType2A, LDA, SType2beta, SType2C, LDC);
01308     if(debug)
01309     {
01310       PrintMatrix(SType1C, N, N, LDC,"SType1C_after_operation", matlab);
01311       PrintMatrix(SType2C, N, N, LDC,"SType2C_after_operation", matlab);
01312     }
01313     GoodTestSubcount += CompareMatrices(SType1C, SType2C, N, N, LDC, TOL);
01314 
01315     delete [] SType1A;
01316     delete [] SType1C;
01317     delete [] SType2A;
01318     delete [] SType2C;
01319   }
01320   GoodTestCount += GoodTestSubcount;
01321   if(verbose || debug) std::cout << "SYRK: " << GoodTestSubcount << " of " << SYRKTESTS << " tests were successful." << std::endl;
01322   if(debug) std::cout << std::endl;
01323   //--------------------------------------------------------------------------------
01324   // End SYRK Tests
01325   //--------------------------------------------------------------------------------
01326 
01327   //--------------------------------------------------------------------------------
01328   // Begin TRMM Tests
01329   //--------------------------------------------------------------------------------
01330   GoodTestSubcount = 0;
01331   for(i = 0; i < TRMMTESTS; i++)
01332     { 
01333       M = GetRandom(MVMIN, MVMAX);
01334       N = GetRandom(MVMIN, MVMAX);
01335 
01336       LDB = GetRandom(MVMIN, MVMAX);
01337       while (LDB < M) {
01338     LDB = GetRandom(MVMIN, MVMAX);
01339       }
01340 
01341       SType1B = new SType1[LDB * N];
01342       SType2B = new SType2[LDB * N];
01343 
01344       SIDE = RandomSIDE();
01345       UPLO = RandomUPLO();
01346       TRANSA = RandomTRANS();
01347       DIAG = RandomDIAG();
01348 
01349       if(Teuchos::ESideChar[SIDE] == 'L')  // The operator is on the left side
01350   {
01351           LDA = GetRandom(MVMIN, MVMAX);
01352           while (LDA < M) {
01353         LDA = GetRandom(MVMIN, MVMAX);
01354           }
01355 
01356     SType1A = new SType1[LDA * M];
01357     SType2A = new SType2[LDA * M];
01358 
01359     for(j = 0; j < M; j++)
01360       {      
01361         if(Teuchos::EUploChar[UPLO] == 'U') {
01362     // The operator is upper triangular, make sure that the entries are
01363     // only in the upper triangular part of A and the diagonal is non-zero.
01364     for(k = 0; k < M; k++) 
01365     {
01366         if(k < j) {
01367             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01368         } else {
01369       SType1A[j*LDA+k] = SType1zero;
01370         }
01371               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01372         if(k == j) {
01373       if (Teuchos::EDiagChar[DIAG] == 'N') {
01374                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01375           while (SType1A[j*LDA+k] == SType1zero) {
01376         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01377           }
01378                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);  
01379           } else {
01380                 SType1A[j*LDA+k] = SType1one;
01381                     SType2A[j*LDA+k] = SType2one;
01382           }     
01383         }
01384     }
01385         } else {
01386     // The operator is lower triangular, make sure that the entries are
01387     // only in the lower triangular part of A and the diagonal is non-zero.
01388     for(k = 0; k < M; k++) 
01389     {
01390         if(k > j) {
01391             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01392         } else {
01393       SType1A[j*LDA+k] = SType1zero;
01394         }
01395               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01396         if(k == j) {
01397       if (Teuchos::EDiagChar[DIAG] == 'N') {
01398                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01399           while (SType1A[j*LDA+k] == SType1zero) {
01400         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01401           }
01402                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01403           } else {
01404                 SType1A[j*LDA+k] = SType1one;
01405                     SType2A[j*LDA+k] = SType2one;
01406       }
01407         }     
01408     } // end for(k=0 ...    
01409         } // end if(UPLO == 'U') ...
01410       } // end for(j=0 ...
01411   } // if(SIDE == 'L') ...
01412       else // The operator is on the right side
01413   {
01414           LDA = GetRandom(MVMIN, MVMAX);
01415           while (LDA < N) {
01416         LDA = GetRandom(MVMIN, MVMAX);
01417           }
01418 
01419     SType1A = new SType1[LDA * N];
01420     SType2A = new SType2[LDA * N];
01421 
01422     for(j = 0; j < N; j++)
01423       {      
01424         if(Teuchos::EUploChar[UPLO] == 'U') {
01425     // The operator is upper triangular, make sure that the entries are
01426     // only in the upper triangular part of A and the diagonal is non-zero.
01427     for(k = 0; k < N; k++) 
01428     {
01429         if(k < j) {
01430             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01431         } else {
01432       SType1A[j*LDA+k] = SType1zero;
01433         }
01434               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01435         if(k == j) {
01436       if (Teuchos::EDiagChar[DIAG] == 'N') {
01437                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01438           while (SType1A[j*LDA+k] == SType1zero) {
01439         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01440           }
01441                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01442           } else {
01443                 SType1A[j*LDA+k] = SType1one;
01444                     SType2A[j*LDA+k] = SType2one;
01445       }
01446         }     
01447     }
01448         } else {
01449     // The operator is lower triangular, make sure that the entries are
01450     // only in the lower triangular part of A and the diagonal is non-zero.
01451     for(k = 0; k < N; k++) 
01452     {
01453         if(k > j) {
01454             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01455         } else {
01456       SType1A[j*LDA+k] = SType1zero;
01457         }
01458               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01459         if(k == j) {
01460       if (Teuchos::EDiagChar[DIAG] == 'N') {
01461                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01462           while (SType1A[j*LDA+k] == SType1zero) {
01463         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01464           }
01465                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01466           } else {
01467                 SType1A[j*LDA+k] = SType1one;
01468                     SType2A[j*LDA+k] = SType2one;
01469       }
01470         }     
01471     } // end for(k=0 ...    
01472         } // end if(UPLO == 'U') ...
01473       } // end for(j=0 ...
01474   } // end if(SIDE == 'L') ...
01475 
01476       // Fill in the right hand side block B.
01477       for(j = 0; j < N; j++) {
01478     for(k = 0; k < M; k++) {
01479       SType1B[j*LDB+k] = GetRandom(-SCALARMAX, SCALARMAX);
01480       SType2B[j*LDB+k] = ConvertType(SType1B[j*LDB+k], convertTo);
01481     }
01482       }
01483       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
01484       SType2alpha = ConvertType(SType1alpha, convertTo);
01485       if(debug)
01486   {
01487     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
01488     std::cout << "SType1alpha = " << SType1alpha << std::endl;
01489     std::cout << "SType2alpha = " << SType2alpha << std::endl;
01490           if(Teuchos::ESideChar[SIDE] == 'L') { 
01491       PrintMatrix(SType1A, M, M, LDA, "SType1A", matlab);
01492       PrintMatrix(SType2A, M, M, LDA, "SType2A", matlab);
01493     } else {
01494       PrintMatrix(SType1A, N, N, LDA, "SType1A", matlab);
01495       PrintMatrix(SType2A, N, N, LDA, "SType2A", matlab);
01496     }
01497     PrintMatrix(SType1B, M, N, LDB,"SType1B_before_operation", matlab);
01498     PrintMatrix(SType2B, M, N, LDB,"SType2B_before_operation", matlab);
01499   }
01500       TotalTestCount++;
01501       SType1BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
01502       SType2BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
01503       if(debug)
01504   {
01505     PrintMatrix(SType1B, M, N, LDB, "SType1B_after_operation", matlab);
01506     PrintMatrix(SType2B, M, N, LDB, "SType2B_after_operation", matlab);
01507   }
01508       GoodTestSubcount += CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
01509       delete [] SType1A;
01510       delete [] SType1B;
01511       delete [] SType2A;
01512       delete [] SType2B;
01513     }
01514   GoodTestCount += GoodTestSubcount;
01515   if(verbose || debug) std::cout << "TRMM: " << GoodTestSubcount << " of " << TRMMTESTS << " tests were successful." << std::endl;
01516   if(debug) std::cout << std::endl;
01517   //--------------------------------------------------------------------------------
01518   // End TRMM Tests
01519   //--------------------------------------------------------------------------------
01520 
01521   //--------------------------------------------------------------------------------
01522   // Begin TRSM Tests
01523   //--------------------------------------------------------------------------------
01524   GoodTestSubcount = 0;
01525   for(i = 0; i < TRSMTESTS; i++)
01526     { 
01527       M = GetRandom(MVMIN, MVMAX);
01528       N = GetRandom(MVMIN, MVMAX);
01529 
01530       LDB = GetRandom(MVMIN, MVMAX);
01531       while (LDB < M) {
01532     LDB = GetRandom(MVMIN, MVMAX);
01533       }
01534 
01535       SType1B = new SType1[LDB * N];
01536       SType2B = new SType2[LDB * N];
01537 
01538       SIDE = RandomSIDE();
01539       UPLO = RandomUPLO();
01540       TRANSA = RandomTRANS();
01541       // Since the entries are integers, we don't want to use the unit diagonal feature,
01542       // this creates ill-conditioned, nearly-singular matrices.
01543       //DIAG = RandomDIAG();  
01544       DIAG = Teuchos::NON_UNIT_DIAG;
01545 
01546       if(Teuchos::ESideChar[SIDE] == 'L')  // The operator is on the left side
01547   {
01548           LDA = GetRandom(MVMIN, MVMAX);
01549           while (LDA < M) {
01550         LDA = GetRandom(MVMIN, MVMAX);
01551           }
01552 
01553     SType1A = new SType1[LDA * M];
01554     SType2A = new SType2[LDA * M];
01555 
01556     for(j = 0; j < M; j++)
01557       {      
01558         if(Teuchos::EUploChar[UPLO] == 'U') {
01559     // The operator is upper triangular, make sure that the entries are
01560     // only in the upper triangular part of A and the diagonal is non-zero.
01561     for(k = 0; k < M; k++) 
01562     {
01563         if(k < j) {
01564             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01565         } else {
01566       SType1A[j*LDA+k] = SType1zero;
01567         }
01568               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01569         if(k == j) {
01570       if (Teuchos::EDiagChar[DIAG] == 'N') {
01571                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01572           while (SType1A[j*LDA+k] == SType1zero) {
01573         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01574           }
01575                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);  
01576           } else {
01577                 SType1A[j*LDA+k] = SType1one;
01578                     SType2A[j*LDA+k] = SType2one;
01579           }     
01580         }
01581     }
01582         } else {
01583     // The operator is lower triangular, make sure that the entries are
01584     // only in the lower triangular part of A and the diagonal is non-zero.
01585     for(k = 0; k < M; k++) 
01586     {
01587         if(k > j) {
01588             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01589         } else {
01590       SType1A[j*LDA+k] = SType1zero;
01591         }
01592               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01593         if(k == j) {
01594       if (Teuchos::EDiagChar[DIAG] == 'N') {
01595                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01596           while (SType1A[j*LDA+k] == SType1zero) {
01597         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01598           }
01599                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01600           } else {
01601                 SType1A[j*LDA+k] = SType1one;
01602                     SType2A[j*LDA+k] = SType2one;
01603       }
01604         }     
01605     } // end for(k=0 ...    
01606         } // end if(UPLO == 'U') ...
01607       } // end for(j=0 ...
01608   } // if(SIDE == 'L') ...
01609       else // The operator is on the right side
01610   {
01611           LDA = GetRandom(MVMIN, MVMAX);
01612           while (LDA < N) {
01613         LDA = GetRandom(MVMIN, MVMAX);
01614           }
01615 
01616     SType1A = new SType1[LDA * N];
01617     SType2A = new SType2[LDA * N];
01618 
01619     for(j = 0; j < N; j++)
01620       {      
01621         if(Teuchos::EUploChar[UPLO] == 'U') {
01622     // The operator is upper triangular, make sure that the entries are
01623     // only in the upper triangular part of A and the diagonal is non-zero.
01624     for(k = 0; k < N; k++) 
01625     {
01626         if(k < j) {
01627             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01628         } else {
01629       SType1A[j*LDA+k] = SType1zero;
01630         }
01631               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01632         if(k == j) {
01633       if (Teuchos::EDiagChar[DIAG] == 'N') {
01634                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01635           while (SType1A[j*LDA+k] == SType1zero) {
01636         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01637           }
01638                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01639           } else {
01640                 SType1A[j*LDA+k] = SType1one;
01641                     SType2A[j*LDA+k] = SType2one;
01642       }
01643         }     
01644     }
01645         } else {
01646     // The operator is lower triangular, make sure that the entries are
01647     // only in the lower triangular part of A and the diagonal is non-zero.
01648     for(k = 0; k < N; k++) 
01649     {
01650         if(k > j) {
01651             SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01652         } else {
01653       SType1A[j*LDA+k] = SType1zero;
01654         }
01655               SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01656         if(k == j) {
01657       if (Teuchos::EDiagChar[DIAG] == 'N') {
01658                 SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01659           while (SType1A[j*LDA+k] == SType1zero) {
01660         SType1A[j*LDA+k] = GetRandom(-SCALARMAX, SCALARMAX);
01661           }
01662                     SType2A[j*LDA+k] = ConvertType(SType1A[j*LDA+k], convertTo);
01663           } else {
01664                 SType1A[j*LDA+k] = SType1one;
01665                     SType2A[j*LDA+k] = SType2one;
01666       }
01667         }     
01668     } // end for(k=0 ...    
01669         } // end if(UPLO == 'U') ...
01670       } // end for(j=0 ...
01671   } // end if(SIDE == 'L') ...
01672 
01673       // Fill in the right hand side block B.
01674       for(j = 0; j < N; j++)
01675   {
01676     for(k = 0; k < M; k++) 
01677       {
01678       SType1B[j*LDB+k] = GetRandom(-SCALARMAX, SCALARMAX);
01679       SType2B[j*LDB+k] = ConvertType(SType1B[j*LDB+k], convertTo);
01680       }
01681   }
01682       
01683       SType1alpha = GetRandom(-SCALARMAX, SCALARMAX);
01684       SType2alpha = ConvertType(SType1alpha, convertTo);
01685       
01686       if(debug)
01687   {
01688     std::cout << "Test #" << TotalTestCount << " --" << std::endl;
01689     std::cout << Teuchos::ESideChar[SIDE] << "\t" 
01690          << Teuchos::EUploChar[UPLO] << "\t" 
01691          << Teuchos::ETranspChar[TRANSA] << "\t" 
01692          << Teuchos::EDiagChar[DIAG] << std::endl;
01693     std::cout << "M="<<M << "\t" << "N="<<N << "\t" << "LDA="<<LDA << "\t" << "LDB="<<LDB << std::endl;
01694     std::cout << "SType1alpha = " << SType1alpha << std::endl;
01695     std::cout << "SType2alpha = " << SType2alpha << std::endl;
01696     if (Teuchos::ESideChar[SIDE] == 'L') {
01697         PrintMatrix(SType1A, M, M, LDA, "SType1A", matlab);
01698         PrintMatrix(SType2A, M, M, LDA, "SType2A", matlab);
01699     } else {
01700         PrintMatrix(SType1A, N, N, LDA, "SType1A", matlab);
01701         PrintMatrix(SType2A, N, N, LDA, "SType2A", matlab);
01702     }
01703     PrintMatrix(SType1B, M, N, LDB, "SType1B_before_operation", matlab);
01704     PrintMatrix(SType2B, M, N, LDB, "SType2B_before_operation", matlab);
01705   }
01706       TotalTestCount++;
01707 
01708       SType1BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
01709       SType2BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
01710  
01711       if(debug)
01712   {
01713     PrintMatrix(SType1B, M, N, LDB, "SType1B_after_operation", matlab);
01714     PrintMatrix(SType2B, M, N, LDB, "SType2B_after_operation", matlab);
01715   }
01716 
01717       if (CompareMatrices(SType1B, SType2B, M, N, LDB, TOL)==0)
01718   std::cout << "FAILED TEST!!!!!!" << std::endl;
01719       GoodTestSubcount += CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
01720 
01721       delete [] SType1A;
01722       delete [] SType1B;
01723       delete [] SType2A;
01724       delete [] SType2B;
01725     }
01726   GoodTestCount += GoodTestSubcount; 
01727   if(verbose || debug) std::cout << "TRSM: " << GoodTestSubcount << " of " << TRSMTESTS << " tests were successful." << std::endl;
01728   if(debug) std::cout << std::endl;
01729   //--------------------------------------------------------------------------------
01730   // End TRSM Tests
01731   //--------------------------------------------------------------------------------
01732 
01733   if((((TotalTestCount - 1) - GoodTestCount) != 0) || (verbose) || (debug))
01734     {
01735       std::cout << GoodTestCount << " of " << (TotalTestCount - 1) << " total tests were successful." << std::endl;
01736     }
01737 
01738   if ((TotalTestCount-1) == GoodTestCount) {
01739     std::cout << "End Result: TEST PASSED" << std::endl;
01740     return 0;
01741   }
01742 
01743   std::cout << "End Result: TEST FAILED" << std::endl;
01744   return (TotalTestCount-GoodTestCount-1);
01745 }
01746 
01747 template<typename TYPE>
01748 TYPE GetRandom(TYPE Low, TYPE High)
01749 {
01750   return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
01751 }
01752 
01753 template<>
01754 int GetRandom(int Low, int High)
01755 {
01756   return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
01757 }
01758 
01759 template<>
01760 double GetRandom(double Low, double High)
01761 {
01762   return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
01763 }
01764 
01765 template<typename TYPE>
01766 void PrintVector(TYPE* Vector, OType Size, std::string Name, bool Matlab)
01767 {
01768   std::cout << Name << " =" << std::endl;
01769   OType i;
01770   if(Matlab) std::cout << "[";
01771   for(i = 0; i < Size; i++)
01772     {
01773       std::cout << Vector[i] << " ";
01774     }
01775   if(Matlab) std::cout << "]";
01776   if(!Matlab)
01777     {
01778       std::cout << std::endl << std::endl;
01779     }
01780   else
01781     {
01782       std::cout << ";" << std::endl;
01783     }
01784 }
01785 
01786 template<typename TYPE>
01787 void PrintMatrix(TYPE* Matrix, OType Rows, OType Columns, OType LDM, std::string Name, bool Matlab)
01788 {
01789   if(!Matlab)
01790     {
01791       std::cout << Name << " =" << std::endl;
01792       OType i, j;
01793       for(i = 0; i < Rows; i++)
01794   {
01795           for(j = 0; j < Columns; j++)
01796       {
01797         std::cout << Matrix[i + (j * LDM)] << " ";
01798       }
01799     std::cout << std::endl;
01800   }
01801       std::cout << std::endl;
01802     }
01803   else
01804     {
01805       std::cout << Name << " = ";
01806       OType i, j;
01807       std::cout << "[";
01808       for(i = 0; i < Rows; i++)
01809         {
01810     std::cout << "[";
01811           for(j = 0; j < Columns; j++)
01812       {
01813         std::cout << Matrix[i + (j * LDM)] << " ";
01814       }
01815     std::cout << "];";
01816   }
01817       std::cout << "];" << std::endl;
01818     }
01819 }
01820 
01821 template<typename TYPE1, typename TYPE2>
01822 bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, typename ScalarTraits<TYPE2>::magnitudeType Tolerance)
01823 {
01824   TYPE2 convertTo = ScalarTraits<TYPE2>::zero();
01825   typename ScalarTraits<TYPE2>::magnitudeType temp = ScalarTraits<TYPE2>::magnitude(Scalar2);
01826   typename ScalarTraits<TYPE2>::magnitudeType temp2 = ScalarTraits<TYPE2>::magnitude(ConvertType(Scalar1, convertTo) - Scalar2);
01827   if (temp != ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero()) {
01828     temp2 /= temp;
01829   }
01830   return( temp2 < Tolerance );
01831 }
01832 
01833 
01834 /*  Function:  CompareVectors
01835     Purpose:   Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
01836 */
01837 template<typename TYPE1, typename TYPE2>
01838 bool CompareVectors(TYPE1* Vector1, TYPE2* Vector2, OType Size, typename ScalarTraits<TYPE2>::magnitudeType Tolerance)
01839 {
01840   TYPE2 convertTo = ScalarTraits<TYPE2>::zero();
01841   TYPE2 temp = ScalarTraits<TYPE2>::zero();
01842   typename ScalarTraits<TYPE2>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01843   typename ScalarTraits<TYPE2>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01844   typename ScalarTraits<TYPE2>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01845   typename ScalarTraits<TYPE2>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01846   OType i;
01847   for(i = 0; i < Size; i++)
01848     {
01849       sum2 += ScalarTraits<TYPE2>::magnitude(ScalarTraits<TYPE2>::conjugate(Vector2[i])*Vector2[i]);
01850       temp = ConvertType(Vector1[i], convertTo) - Vector2[i];
01851       sum += ScalarTraits<TYPE2>::magnitude(ScalarTraits<TYPE2>::conjugate(temp)*temp);
01852     }
01853   temp2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum2);
01854   if (temp2 != ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero())
01855     temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum)/temp2;
01856   else
01857     temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum);
01858   if (temp3 > Tolerance )
01859     return false;
01860   else
01861     return true;
01862 }
01863 
01864 /*  Function:  CompareMatrices
01865     Purpose:   Compares the difference between two matrices using relative frobenius-norms, i.e. ||M_1-M_2||_F/||M_2||_F
01866 */
01867 template<typename TYPE1, typename TYPE2>
01868 bool CompareMatrices(TYPE1* Matrix1, TYPE2* Matrix2, OType Rows, OType Columns, OType LDM, typename ScalarTraits<TYPE2>::magnitudeType Tolerance)
01869 {
01870   TYPE2 convertTo = ScalarTraits<TYPE2>::zero();
01871   TYPE2 temp = ScalarTraits<TYPE2>::zero();
01872   typename ScalarTraits<TYPE2>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01873   typename ScalarTraits<TYPE2>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01874   typename ScalarTraits<TYPE2>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01875   typename ScalarTraits<TYPE2>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero();
01876   OType i,j;
01877   for(j = 0; j < Columns; j++)
01878     {
01879       for(i = 0; i < Rows; i++)
01880   {
01881     sum2 = ScalarTraits<TYPE2>::magnitude(ScalarTraits<TYPE2>::conjugate(Matrix2[j*LDM + i])*Matrix2[j*LDM + i]);
01882     temp = ConvertType(Matrix1[j*LDM + i],convertTo) - Matrix2[j*LDM + i]; 
01883     sum = ScalarTraits<TYPE2>::magnitude(ScalarTraits<TYPE2>::conjugate(temp)*temp);
01884   }
01885     }
01886   temp2 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum2);
01887   if (temp2 != ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::zero())
01888     temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum)/temp2;
01889   else
01890     temp3 = ScalarTraits<typename ScalarTraits<TYPE2>::magnitudeType>::squareroot(sum);
01891   if (temp3 > Tolerance)
01892     return false;
01893   else
01894     return true;
01895 }
01896 
01897 
01898 template<typename TYPE1, typename TYPE2>
01899 TYPE2 ConvertType(TYPE1 T1, TYPE2 T2)
01900 {
01901   return static_cast<TYPE2>(T1);
01902 }
01903 
01904 Teuchos::ESide RandomSIDE()
01905 {
01906   Teuchos::ESide result;
01907   int r = GetRandom(1, 2);
01908   if(r == 1)
01909     {
01910       result = Teuchos::LEFT_SIDE;
01911     }
01912   else
01913     {
01914       result = Teuchos::RIGHT_SIDE;
01915     }
01916   return result;
01917 }
01918 
01919 Teuchos::EUplo RandomUPLO()
01920 {
01921   Teuchos::EUplo result;
01922   int r = GetRandom(1, 2);
01923   if(r == 1)
01924     {
01925       result = Teuchos::UPPER_TRI;
01926     }
01927   else
01928     {
01929       result = Teuchos::LOWER_TRI;
01930     }
01931   return result;
01932 }
01933 
01934 Teuchos::ETransp RandomTRANS()
01935 {
01936   Teuchos::ETransp result;
01937   int r = GetRandom(1, 4);
01938   if(r == 1 || r == 2)
01939     {
01940       result = Teuchos::NO_TRANS;
01941     }
01942   else if(r == 3)
01943     {
01944       result = Teuchos::TRANS;
01945     }
01946   else
01947     {
01948       result = Teuchos::CONJ_TRANS;
01949     }
01950   return result;
01951 }
01952 
01953 Teuchos::EDiag RandomDIAG()
01954 {
01955   Teuchos::EDiag result;
01956   int r = GetRandom(1, 2);
01957   if(r == 1)
01958     {
01959       result = Teuchos::NON_UNIT_DIAG;
01960     }
01961   else
01962     {
01963       result = Teuchos::UNIT_DIAG;
01964     }
01965   return result;
01966 }
01967 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines