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