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

Generated on Tue Oct 20 10:13:59 2009 for Teuchos Package Browser (Single Doxygen Collection) by  doxygen 1.6.1