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