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