|
Teuchos Package Browser (Single Doxygen Collection)
Version of the Day
|
|
Teuchos Package Browser (Single Doxygen Collection)
Version of the Day
|
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 << "FAILED TEST!!!!!!" << std::endl; 00594 GoodTestSubcount += CompareScalars(OType1DOTresult, OType2DOTresult, TOL); 00595 00596 delete [] OType1x; 00597 delete [] OType1y; 00598 delete [] OType2x; 00599 delete [] OType2y; 00600 } 00601 GoodTestCount += GoodTestSubcount; 00602 if(verbose || debug) std::cout << "DOT: " << GoodTestSubcount << " of " << DOTTESTS << " tests were successful." << std::endl; 00603 if(debug) std::cout << std::endl; 00604 //-------------------------------------------------------------------------------- 00605 // End DOT Tests 00606 //-------------------------------------------------------------------------------- 00607 00608 //-------------------------------------------------------------------------------- 00609 // Begin NRM2 Tests 00610 //-------------------------------------------------------------------------------- 00611 GoodTestSubcount = 0; 00612 for(i = 0; i < NRM2TESTS; i++) 00613 { 00614 incx1 = GetRandom(1, MVMAX); 00615 incx2 = ConvertType( incx1, incx2 ); 00616 M1 = GetRandom(MVMIN, MVMAX); 00617 M2 = ConvertType( M1, M2 ); 00618 OType1x = new SType[M1*incx1]; 00619 OType2x = new SType[M2*incx2]; 00620 for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++) 00621 { 00622 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00623 OType2x[j2] = OType1x[j1]; 00624 } 00625 if(debug) 00626 { 00627 std::cout << "Test #" << TotalTestCount << " -- NRM2 -- " << std::endl; 00628 PrintVector(OType1x, M1*incx1, "OType1x", matlab); 00629 PrintVector(OType2x, M2*incx2, "OType2x", matlab); 00630 } 00631 TotalTestCount++; 00632 OType1NRM2result = OType1BLAS.NRM2(M1, OType1x, incx1); 00633 OType2NRM2result = OType2BLAS.NRM2(M2, OType2x, incx2); 00634 if(debug) 00635 { 00636 std::cout << "OType1 NRM2 result: " << OType1NRM2result << std::endl; 00637 std::cout << "OType2 NRM2 result: " << OType2NRM2result << std::endl; 00638 } 00639 if (CompareScalars(OType1NRM2result, OType2NRM2result, TOL)==0) 00640 std::cout << "FAILED TEST!!!!!!" << std::endl; 00641 GoodTestSubcount += CompareScalars(OType1NRM2result, OType2NRM2result, TOL); 00642 00643 delete [] OType1x; 00644 delete [] OType2x; 00645 } 00646 GoodTestCount += GoodTestSubcount; if(verbose || debug) std::cout << "NRM2: " << GoodTestSubcount << " of " << NRM2TESTS << " tests were successful." << std::endl; 00647 if(debug) std::cout << std::endl; 00648 //-------------------------------------------------------------------------------- 00649 // End NRM2 Tests 00650 //-------------------------------------------------------------------------------- 00651 00652 //-------------------------------------------------------------------------------- 00653 // Begin SCAL Tests 00654 //-------------------------------------------------------------------------------- 00655 GoodTestSubcount = 0; 00656 for(i = 0; i < SCALTESTS; i++) 00657 { 00658 // These will only test for the case that the increment is > 0, the 00659 // templated case can handle when incx < 0, but the blas library doesn't 00660 // seem to be able to on some machines. 00661 incx1 = GetRandom(1, MVMAX); 00662 incx2 = ConvertType( incx1, incx2 ); 00663 M1 = GetRandom(MVMIN, MVMAX); 00664 M2 = ConvertType( M1, M2 ); 00665 OType1x = new SType[M1*incx1]; 00666 OType2x = new SType[M2*incx2]; 00667 for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++) 00668 { 00669 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00670 OType2x[j2] = OType1x[j1]; 00671 } 00672 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 00673 OType2alpha = OType1alpha; 00674 if(debug) 00675 { 00676 std::cout << "Test #" << TotalTestCount << " -- SCAL -- " << std::endl; 00677 std::cout << "OType1alpha = " << OType1alpha << std::endl; 00678 std::cout << "OType2alpha = " << OType2alpha << std::endl; 00679 PrintVector(OType1x, M1*incx1, "OType1x_before_operation", matlab); 00680 PrintVector(OType2x, M2*incx2, "OType2x_before_operation", matlab); 00681 } 00682 TotalTestCount++; 00683 OType1BLAS.SCAL(M1, OType1alpha, OType1x, incx1); 00684 OType2BLAS.SCAL(M2, OType2alpha, OType2x, incx2); 00685 if(debug) 00686 { 00687 PrintVector(OType1x, M1*incx1, "OType1x_after_operation", matlab); 00688 PrintVector(OType2x, M2*incx2, "OType2x_after_operation", matlab); 00689 } 00690 if (CompareVectors(OType1x, M1*incx1, OType2x, M2*incx2, TOL)==0) 00691 std::cout << "FAILED TEST!!!!!!" << std::endl; 00692 GoodTestSubcount += CompareVectors(OType1x, M1*incx1, OType2x, M2*incx2, TOL); 00693 00694 delete [] OType1x; 00695 delete [] OType2x; 00696 } 00697 GoodTestCount += GoodTestSubcount; 00698 if(verbose || debug) std::cout << "SCAL: " << GoodTestSubcount << " of " << SCALTESTS << " tests were successful." << std::endl; 00699 if(debug) std::cout << std::endl; 00700 //-------------------------------------------------------------------------------- 00701 // End SCAL Tests 00702 //-------------------------------------------------------------------------------- 00703 00704 //-------------------------------------------------------------------------------- 00705 // Begin IAMAX Tests 00706 //-------------------------------------------------------------------------------- 00707 GoodTestSubcount = 0; 00708 for(i = 0; i < IAMAXTESTS; i++) 00709 { 00710 incx1 = GetRandom(1, MVMAX); 00711 incx2 = ConvertType( incx1, incx2 ); 00712 M1 = GetRandom(MVMIN, MVMAX); 00713 M2 = ConvertType( M1, M2 ); 00714 OType1x = new SType[M1*incx1]; 00715 OType2x = new SType[M2*incx2]; 00716 for(j1 = 0, j2 = 0; j1 < M1*incx1; j1++, j2++) 00717 { 00718 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00719 OType2x[j2] = OType1x[j1]; 00720 } 00721 if(debug) 00722 { 00723 std::cout << "Test #" << TotalTestCount << " -- IAMAX -- " << std::endl; 00724 PrintVector(OType1x, M1*incx1, "OType1x", matlab); 00725 PrintVector(OType2x, M2*incx2, "OType2x", matlab); 00726 } 00727 TotalTestCount++; 00728 OType1IAMAXresult = OType1BLAS.IAMAX(M1, OType1x, incx1); 00729 OType2IAMAXresult = OType2BLAS.IAMAX(M2, OType2x, incx2); 00730 if(debug) 00731 { 00732 std::cout << "OType1 IAMAX result: " << OType1IAMAXresult << std::endl; 00733 std::cout << "OType2 IAMAX result: " << OType2IAMAXresult << std::endl; 00734 } 00735 if (OType1IAMAXresult != OType2IAMAXresult) 00736 std::cout << "FAILED TEST!!!!!!" << std::endl; 00737 GoodTestSubcount += (OType1IAMAXresult == OType2IAMAXresult); 00738 00739 delete [] OType1x; 00740 delete [] OType2x; 00741 } 00742 GoodTestCount += GoodTestSubcount; 00743 if(verbose || debug) std::cout << "IAMAX: " << GoodTestSubcount << " of " << IAMAXTESTS << " tests were successful." << std::endl; 00744 if(debug) std::cout << std::endl; 00745 //-------------------------------------------------------------------------------- 00746 // End IAMAX Tests 00747 //-------------------------------------------------------------------------------- 00748 00749 //-------------------------------------------------------------------------------- 00750 // BEGIN LEVEL 2 BLAS TESTS 00751 //-------------------------------------------------------------------------------- 00752 // Begin GEMV Tests 00753 //-------------------------------------------------------------------------------- 00754 GoodTestSubcount = 0; 00755 for(i = 0; i < GEMVTESTS; i++) 00756 { 00757 // The parameters used to construct the test problem are chosen to be 00758 // valid parameters, so the GEMV routine won't bomb out. 00759 incx1 = GetRandom(1, MVMAX); 00760 while (incx1 == 0) { 00761 incx1 = GetRandom(1, MVMAX); 00762 } 00763 incy1 = GetRandom(1, MVMAX); 00764 while (incy1 == 0) { 00765 incy1 = GetRandom(1, MVMAX); 00766 } 00767 incx2 = ConvertType( incx1, incx2 ); 00768 incy2 = ConvertType( incy1, incy2 ); 00769 M1 = GetRandom(MVMIN, MVMAX); 00770 N1 = GetRandom(MVMIN, MVMAX); 00771 M2 = ConvertType( M1, M2 ); 00772 N2 = ConvertType( N1, N2 ); 00773 00774 TRANS = RandomTRANS(); 00775 OType1 M2_1 = 0, N2_1 = 0; 00776 if (Teuchos::ETranspChar[TRANS] == 'N') { 00777 M2_1 = M1*std::abs(incy1); 00778 N2_1 = N1*std::abs(incx1); 00779 } else { 00780 M2_1 = N1*std::abs(incy1); 00781 N2_1 = M1*std::abs(incx1); 00782 } 00783 OType2 M2_2 = ConvertType( M2_1, M2_2 ); 00784 OType2 N2_2 = ConvertType( N2_1, N2_2 ); 00785 00786 LDA1 = GetRandom(MVMIN, MVMAX); 00787 while (LDA1 < M1) { 00788 LDA1 = GetRandom(MVMIN, MVMAX); 00789 } 00790 LDA2 = ConvertType( LDA1, LDA2 ); 00791 00792 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 00793 OType1beta = GetRandom(-SCALARMAX, SCALARMAX); 00794 OType2alpha = OType1alpha; 00795 OType2beta = OType1beta; 00796 00797 OType1A = new SType[LDA1 * N1]; 00798 OType1x = new SType[N2_1]; 00799 OType1y = new SType[M2_1]; 00800 OType2A = new SType[LDA2 * N2]; 00801 OType2x = new SType[N2_2]; 00802 OType2y = new SType[M2_2]; 00803 00804 for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) 00805 { 00806 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00807 OType2A[j2] = OType1A[j1]; 00808 } 00809 for(j1 = 0, j2 = 0; j1 < N2_1; j1++, j2++) 00810 { 00811 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00812 OType2x[j2] = OType1x[j1]; 00813 } 00814 for(j1 = 0, j2 = 0; j1 < M2_1; j1++, j2++) 00815 { 00816 OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00817 OType2y[j2] = OType1y[j1]; 00818 } 00819 if(debug) 00820 { 00821 std::cout << "Test #" << TotalTestCount << " -- GEMV -- " << std::endl; 00822 std::cout << "TRANS = " << Teuchos::ETranspChar[TRANS] << std::endl; 00823 std::cout << "OType1alpha = " << OType1alpha << std::endl; 00824 std::cout << "OType2alpha = " << OType2alpha << std::endl; 00825 std::cout << "OType1beta = " << OType1beta << std::endl; 00826 std::cout << "OType2beta = " << OType2beta << std::endl; 00827 PrintMatrix(OType1A, M1, N1, LDA1, "OType1A", matlab); 00828 PrintVector(OType1x, N2_1, "OType1x", matlab); 00829 PrintVector(OType1y, M2_1, "OType1y_before_operation", matlab); 00830 PrintMatrix(OType2A, M2, N2, LDA2, "OType2A", matlab); 00831 PrintVector(OType2x, N2_2, "OType2x", matlab); 00832 PrintVector(OType2y, M2_2, "OType2y_before_operation", matlab); 00833 } 00834 TotalTestCount++; 00835 OType1BLAS.GEMV(TRANS, M1, N1, OType1alpha, OType1A, LDA1, OType1x, incx1, OType1beta, OType1y, incy1); 00836 OType2BLAS.GEMV(TRANS, M2, N2, OType2alpha, OType2A, LDA2, OType2x, incx2, OType2beta, OType2y, incy2); 00837 if(debug) 00838 { 00839 PrintVector(OType1y, M2_1, "OType1y_after_operation", matlab); 00840 PrintVector(OType2y, M2_2, "OType2y_after_operation", matlab); 00841 } 00842 if (CompareVectors(OType1y, M2_1, OType2y, M2_2, TOL)==0) 00843 std::cout << "FAILED TEST!!!!!!" << std::endl; 00844 GoodTestSubcount += CompareVectors(OType1y, M2_1, OType2y, M2_2, TOL); 00845 00846 delete [] OType1A; 00847 delete [] OType1x; 00848 delete [] OType1y; 00849 delete [] OType2A; 00850 delete [] OType2x; 00851 delete [] OType2y; 00852 } 00853 GoodTestCount += GoodTestSubcount; 00854 if(verbose || debug) std::cout << "GEMV: " << GoodTestSubcount << " of " << GEMVTESTS << " tests were successful." << std::endl; 00855 if(debug) std::cout << std::endl; 00856 //-------------------------------------------------------------------------------- 00857 // End GEMV Tests 00858 //-------------------------------------------------------------------------------- 00859 00860 //-------------------------------------------------------------------------------- 00861 // Begin TRMV Tests 00862 //-------------------------------------------------------------------------------- 00863 GoodTestSubcount = 0; 00864 for(i = 0; i < TRMVTESTS; i++) 00865 { 00866 UPLO = RandomUPLO(); 00867 TRANSA = RandomTRANS(); 00868 00869 // Since the entries are integers, we don't want to use the unit diagonal feature, 00870 // this creates ill-conditioned, nearly-singular matrices. 00871 //DIAG = RandomDIAG(); 00872 DIAG = Teuchos::NON_UNIT_DIAG; 00873 00874 N1 = GetRandom(MVMIN, MVMAX); 00875 N2 = ConvertType( N1, N2 ); 00876 incx1 = GetRandom(1, MVMAX); 00877 while (incx1 == 0) { 00878 incx1 = GetRandom(1, MVMAX); 00879 } 00880 incx2 = ConvertType( incx1, incx2 ); 00881 OType1x = new SType[N1*std::abs(incx1)]; 00882 OType2x = new SType[N2*std::abs(incx2)]; 00883 00884 for(j1 = 0, j2 = 0; j1 < N1*std::abs(incx1); j1++, j2++) 00885 { 00886 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 00887 OType2x[j2] = OType1x[j1]; 00888 } 00889 00890 LDA1 = GetRandom(MVMIN, MVMAX); 00891 while (LDA1 < N1) { 00892 LDA1 = GetRandom(MVMIN, MVMAX); 00893 } 00894 LDA2 = ConvertType( LDA1, LDA2 ); 00895 OType1A = new SType[LDA1 * N1]; 00896 OType2A = new SType[LDA2 * N2]; 00897 00898 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) 00899 { 00900 if(Teuchos::EUploChar[UPLO] == 'U') { 00901 // The operator is upper triangular, make sure that the entries are 00902 // only in the upper triangular part of A and the diagonal is non-zero. 00903 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 00904 { 00905 if(k1 < j1) { 00906 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00907 } else { 00908 OType1A[j1*LDA1+k1] = STypezero; 00909 } 00910 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 00911 if(k1 == j1) { 00912 if (Teuchos::EDiagChar[DIAG] == 'N') { 00913 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00914 while (OType1A[j1*LDA1+k1] == STypezero) { 00915 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00916 } 00917 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 00918 } else { 00919 OType1A[j1*LDA1+k1] = STypeone; 00920 OType2A[j2*LDA2+k2] = STypeone; 00921 } 00922 } 00923 } 00924 } else { 00925 // The operator is lower triangular, make sure that the entries are 00926 // only in the lower triangular part of A and the diagonal is non-zero. 00927 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 00928 { 00929 if(k1 > j1) { 00930 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00931 } else { 00932 OType1A[j1*LDA1+k1] = STypezero; 00933 } 00934 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 00935 if(k1 == j1) { 00936 if (Teuchos::EDiagChar[DIAG] == 'N') { 00937 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00938 while (OType1A[j1*LDA1+k1] == STypezero) { 00939 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 00940 } 00941 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 00942 } else { 00943 OType1A[j1*LDA1+k1] = STypeone; 00944 OType2A[j2*LDA2+k2] = STypeone; 00945 } 00946 } 00947 } // end for(k=0 ... 00948 } // end if(UPLO == 'U') ... 00949 } // end for(j=0 ... for(j = 0; j < N*N; j++) 00950 00951 if(debug) 00952 { 00953 std::cout << "Test #" << TotalTestCount << " -- TRMV -- " << std::endl; 00954 std::cout << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t" 00955 << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t" 00956 << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl; 00957 PrintMatrix(OType1A, N1, N1, LDA1,"OType1A", matlab); 00958 PrintVector(OType1x, N1*incx1, "OType1x_before_operation", matlab); 00959 PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab); 00960 PrintVector(OType2x, N2*incx2, "OType2x_before_operation", matlab); 00961 } 00962 TotalTestCount++; 00963 OType1BLAS.TRMV(UPLO, TRANSA, DIAG, N1, OType1A, LDA1, OType1x, incx1); 00964 OType2BLAS.TRMV(UPLO, TRANSA, DIAG, N2, OType2A, LDA2, OType2x, incx2); 00965 if(debug) 00966 { 00967 PrintVector(OType1x, N1*incx1, "OType1x_after_operation", matlab); 00968 PrintVector(OType2x, N2*incx2, "OType2x_after_operation", matlab); 00969 } 00970 if (CompareVectors(OType1x, std::abs(N1*incx1), OType2x, std::abs(N2*incx2), TOL)==0) 00971 std::cout << "FAILED TEST!!!!!!" << std::endl; 00972 GoodTestSubcount += CompareVectors(OType1x, std::abs(N1*incx1), OType2x, std::abs(N2*incx2), TOL); 00973 00974 delete [] OType1A; 00975 delete [] OType1x; 00976 delete [] OType2A; 00977 delete [] OType2x; 00978 } 00979 GoodTestCount += GoodTestSubcount; 00980 if(verbose || debug) std::cout << "TRMV: " << GoodTestSubcount << " of " << TRMVTESTS << " tests were successful." << std::endl; 00981 if(debug) std::cout << std::endl; 00982 //-------------------------------------------------------------------------------- 00983 // End TRMV Tests 00984 //-------------------------------------------------------------------------------- 00985 00986 //-------------------------------------------------------------------------------- 00987 // Begin GER Tests 00988 //-------------------------------------------------------------------------------- 00989 GoodTestSubcount = 0; 00990 for(i = 0; i < GERTESTS; i++) 00991 { 00992 incx1 = GetRandom(1, MVMAX); 00993 while (incx1 == 0) { 00994 incx1 = GetRandom(1, MVMAX); 00995 } 00996 incy1 = GetRandom(1, MVMAX); 00997 while (incy1 == 0) { 00998 incy1 = GetRandom(1, MVMAX); 00999 } 01000 incx2 = ConvertType( incx1, incx2 ); 01001 incy2 = ConvertType( incy1, incy2 ); 01002 M1 = GetRandom(MVMIN, MVMAX); 01003 N1 = GetRandom(MVMIN, MVMAX); 01004 M2 = ConvertType( M1, M2 ); 01005 N2 = ConvertType( N1, N2 ); 01006 01007 LDA1 = GetRandom(MVMIN, MVMAX); 01008 while (LDA1 < M1) { 01009 LDA1 = GetRandom(MVMIN, MVMAX); 01010 } 01011 LDA2 = ConvertType( LDA1, LDA2 ); 01012 01013 OType1A = new SType[LDA1 * N1]; 01014 OType1x = new SType[M1*std::abs(incx1)]; 01015 OType1y = new SType[N1*std::abs(incy1)]; 01016 OType2A = new SType[LDA2 * N2]; 01017 OType2x = new SType[M2*std::abs(incx2)]; 01018 OType2y = new SType[N2*std::abs(incy2)]; 01019 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01020 OType2alpha = OType1alpha; 01021 for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) 01022 { 01023 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01024 OType2A[j2] = OType1A[j1]; 01025 } 01026 for(j1 = 0, j2 = 0; j1 < std::abs(M1*incx1); j1++, j2++) 01027 { 01028 OType1x[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01029 OType2x[j2] = OType1x[j1]; 01030 } 01031 for(j1 = 0, j2 = 0; j1 < std::abs(N1*incy1); j1++, j2++) 01032 { 01033 OType1y[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01034 OType2y[j2] = OType1y[j1]; 01035 } 01036 if(debug) 01037 { 01038 std::cout << "Test #" << TotalTestCount << " -- GER -- " << std::endl; 01039 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01040 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01041 PrintMatrix(OType1A, M1, N1, LDA1,"OType1A_before_operation", matlab); 01042 PrintVector(OType1x, std::abs(M1*incx1), "OType1x", matlab); 01043 PrintVector(OType1y, std::abs(N1*incy1), "OType1y", matlab); 01044 PrintMatrix(OType2A, M2, N2, LDA2,"OType2A_before_operation", matlab); 01045 PrintVector(OType2x, std::abs(M2*incx2), "OType2x", matlab); 01046 PrintVector(OType2y, std::abs(N2*incy2), "OType2y", matlab); 01047 } 01048 TotalTestCount++; 01049 OType1BLAS.GER(M1, N1, OType1alpha, OType1x, incx1, OType1y, incy1, OType1A, LDA1); 01050 OType2BLAS.GER(M2, N2, OType2alpha, OType2x, incx2, OType2y, incy2, OType2A, LDA2); 01051 if(debug) 01052 { 01053 PrintMatrix(OType1A, M1, N1, LDA1, "OType1A_after_operation", matlab); 01054 PrintMatrix(OType2A, M2, N2, LDA2, "OType2A_after_operation", matlab); 01055 } 01056 if (CompareMatrices(OType1A, M1, N1, LDA1, OType2A, M2, N2, LDA2, TOL)==0) 01057 std::cout << "FAILED TEST!!!!!!" << std::endl; 01058 GoodTestSubcount += CompareMatrices(OType1A, M1, N1, LDA1, OType2A, M2, N2, LDA2, TOL); 01059 01060 delete [] OType1A; 01061 delete [] OType1x; 01062 delete [] OType1y; 01063 delete [] OType2A; 01064 delete [] OType2x; 01065 delete [] OType2y; 01066 } 01067 GoodTestCount += GoodTestSubcount; 01068 if(verbose || debug) std::cout << "GER: " << GoodTestSubcount << " of " << GERTESTS << " tests were successful." << std::endl; 01069 if(debug) std::cout << std::endl; 01070 //-------------------------------------------------------------------------------- 01071 // End GER Tests 01072 //-------------------------------------------------------------------------------- 01073 01074 //-------------------------------------------------------------------------------- 01075 // BEGIN LEVEL 3 BLAS TESTS 01076 //-------------------------------------------------------------------------------- 01077 // Begin GEMM Tests 01078 //-------------------------------------------------------------------------------- 01079 GoodTestSubcount = 0; 01080 for(i = 0; i < GEMMTESTS; i++) 01081 { 01082 TRANSA = RandomTRANS(); 01083 TRANSB = RandomTRANS(); 01084 M1 = GetRandom(MVMIN, MVMAX); 01085 N1 = GetRandom(MVMIN, MVMAX); 01086 P1 = GetRandom(MVMIN, MVMAX); 01087 M2 = ConvertType( M1, M2 ); 01088 N2 = ConvertType( N1, N2 ); 01089 P2 = ConvertType( P1, P2 ); 01090 01091 if(debug) { 01092 std::cout << "Test #" << TotalTestCount << " -- GEMM -- " << std::endl; 01093 std::cout << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t" 01094 << "TRANSB = " << Teuchos::ETranspChar[TRANSB] << std::endl; 01095 } 01096 LDA1 = GetRandom(MVMIN, MVMAX); 01097 if (Teuchos::ETranspChar[TRANSA] == 'N') { 01098 while (LDA1 < M1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01099 LDA2 = ConvertType( LDA1, LDA2 ); 01100 OType1A = new SType[LDA1 * P1]; 01101 OType2A = new SType[LDA2 * P2]; 01102 for(j1 = 0, j2 = 0; j1 < LDA1 * P1; j1++, j2++) 01103 { 01104 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01105 OType2A[j2] = OType1A[j1]; 01106 } 01107 if (debug) { 01108 PrintMatrix(OType1A, M1, P1, LDA1, "OType1A", matlab); 01109 PrintMatrix(OType2A, M2, P2, LDA2, "OType2A", matlab); 01110 } 01111 } else { 01112 while (LDA1 < P1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01113 LDA2 = ConvertType( LDA1, LDA2 ); 01114 OType1A = new SType[LDA1 * M1]; 01115 OType2A = new SType[LDA1 * M1]; 01116 for(j1 = 0, j2 = 0; j1 < LDA1 * M1; j1++, j2++) 01117 { 01118 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01119 OType2A[j2] = OType1A[j1]; 01120 } 01121 if (debug) { 01122 PrintMatrix(OType1A, P1, M1, LDA1, "OType1A", matlab); 01123 PrintMatrix(OType2A, P2, M2, LDA2, "OType2A", matlab); 01124 } 01125 } 01126 01127 LDB1 = GetRandom(MVMIN, MVMAX); 01128 if (Teuchos::ETranspChar[TRANSB] == 'N') { 01129 while (LDB1 < P1) { LDB1 = GetRandom(MVMIN, MVMAX); } 01130 LDB2 = ConvertType( LDB1, LDB2 ); 01131 OType1B = new SType[LDB1 * N1]; 01132 OType2B = new SType[LDB2 * N2]; 01133 for(j1 = 0, j2 = 0; j1 < LDB1 * N1; j1++, j2++) 01134 { 01135 OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01136 OType2B[j2] = OType1B[j1]; 01137 } 01138 if (debug) { 01139 PrintMatrix(OType1B, P1, N1, LDB1,"OType1B", matlab); 01140 PrintMatrix(OType2B, P2, N2, LDB2,"OType2B", matlab); 01141 } 01142 } else { 01143 while (LDB1 < N1) { LDB1 = GetRandom(MVMIN, MVMAX); } 01144 LDB2 = ConvertType( LDB1, LDB2 ); 01145 OType1B = new SType[LDB1 * P1]; 01146 OType2B = new SType[LDB2 * P2]; 01147 for(j1 = 0, j2 = 0; j1 < LDB1 * P1; j1++, j2++) 01148 { 01149 OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01150 OType2B[j2] = OType1B[j1]; 01151 } 01152 if (debug) { 01153 PrintMatrix(OType1B, N1, P1, LDB1,"OType1B", matlab); 01154 PrintMatrix(OType2B, N2, P2, LDB2,"OType2B", matlab); 01155 } 01156 } 01157 01158 LDC1 = GetRandom(MVMIN, MVMAX); 01159 while (LDC1 < M1) { LDC1 = GetRandom(MVMIN, MVMAX); } 01160 LDC2 = ConvertType( LDC1, LDC2 ); 01161 OType1C = new SType[LDC1 * N1]; 01162 OType2C = new SType[LDC2 * N2]; 01163 for(j1 = 0, j2 = 0; j1 < LDC1 * N1; j1++, j2++) { 01164 OType1C[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01165 OType2C[j2] = OType1C[j1]; 01166 } 01167 if(debug) 01168 { 01169 PrintMatrix(OType1C, M1, N1, LDC1, "OType1C_before_operation", matlab); 01170 PrintMatrix(OType2C, M2, N2, LDC2, "OType2C_before_operation", matlab); 01171 } 01172 01173 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01174 OType1beta = GetRandom(-SCALARMAX, SCALARMAX); 01175 OType2alpha = OType1alpha; 01176 OType2beta = OType1beta; 01177 01178 TotalTestCount++; 01179 OType1BLAS.GEMM(TRANSA, TRANSB, M1, N1, P1, OType1alpha, OType1A, LDA1, OType1B, LDB1, OType1beta, OType1C, LDC1); 01180 OType2BLAS.GEMM(TRANSA, TRANSB, M2, N2, P2, OType2alpha, OType2A, LDA2, OType2B, LDB2, OType2beta, OType2C, LDC2); 01181 if(debug) 01182 { 01183 std::cout << "M1="<<M1 << "\t" << "N1="<<N1 << "\t" << "P1 = " << P1 01184 << "\t" << "LDA1="<<LDA1 << "\t" << "LDB1="<<LDB1 << "\t" << "LDC1=" << LDC1 << std::endl; 01185 std::cout << "M2="<<M2 << "\t" << "N2="<<N2 << "\t" << "P2 = " << P2 01186 << "\t" << "LDA2="<<LDA2 << "\t" << "LDB2="<<LDB2 << "\t" << "LDC2=" << LDC2 << std::endl; 01187 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01188 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01189 std::cout << "OType1beta = " << OType1beta << std::endl; 01190 std::cout << "OType2beta = " << OType2beta << std::endl; 01191 PrintMatrix(OType1C, M1, N1, LDC1, "OType1C_after_operation", matlab); 01192 PrintMatrix(OType2C, M2, N2, LDC2, "OType2C_after_operation", matlab); 01193 } 01194 if (CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL)==0) 01195 std::cout << "FAILED TEST!!!!!!" << std::endl; 01196 GoodTestSubcount += CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL); 01197 01198 delete [] OType1A; 01199 delete [] OType1B; 01200 delete [] OType1C; 01201 delete [] OType2A; 01202 delete [] OType2B; 01203 delete [] OType2C; 01204 } 01205 GoodTestCount += GoodTestSubcount; 01206 if(verbose || debug) std::cout << "GEMM: " << GoodTestSubcount << " of " << GEMMTESTS << " tests were successful." << std::endl; 01207 if(debug) std::cout << std::endl; 01208 //-------------------------------------------------------------------------------- 01209 // End GEMM Tests 01210 //-------------------------------------------------------------------------------- 01211 01212 //-------------------------------------------------------------------------------- 01213 // Begin SYMM Tests 01214 //-------------------------------------------------------------------------------- 01215 GoodTestSubcount = 0; 01216 for(i = 0; i < SYMMTESTS; i++) 01217 { 01218 M1 = GetRandom(MVMIN, MVMAX); 01219 N1 = GetRandom(MVMIN, MVMAX); 01220 M2 = ConvertType( M1, M2 ); 01221 N2 = ConvertType( N1, N2 ); 01222 SIDE = RandomSIDE(); 01223 UPLO = RandomUPLO(); 01224 01225 LDA1 = GetRandom(MVMIN, MVMAX); 01226 if(Teuchos::ESideChar[SIDE] == 'L') { 01227 while (LDA1 < M1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01228 LDA2 = ConvertType( LDA1, LDA2 ); 01229 OType1A = new SType[LDA1 * M1]; 01230 OType2A = new SType[LDA2 * M2]; 01231 for(j1 = 0, j2 = 0; j1 < LDA1 * M1; j1++, j2++) { 01232 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01233 OType2A[j2] = OType1A[j1]; 01234 } 01235 } else { 01236 while (LDA1 < N1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01237 LDA2 = ConvertType( LDA1, LDA2 ); 01238 OType1A = new SType[LDA1 * N1]; 01239 OType2A = new SType[LDA2 * N2]; 01240 for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) { 01241 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01242 OType2A[j2] = OType1A[j1]; 01243 } 01244 } 01245 01246 LDB1 = GetRandom(MVMIN, MVMAX); 01247 while (LDB1 < M1) { LDB1 = GetRandom(MVMIN, MVMAX); } 01248 LDB2 = ConvertType( LDB1, LDB2 ); 01249 OType1B = new SType[LDB1 * N1]; 01250 OType2B = new SType[LDB2 * N2]; 01251 for(j1 = 0, j2 = 0; j1 < LDB1 * N1; j1++, j2++) { 01252 OType1B[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01253 OType2B[j2] = OType1B[j1]; 01254 } 01255 01256 LDC1 = GetRandom(MVMIN, MVMAX); 01257 while (LDC1 < M1) { LDC1 = GetRandom(MVMIN, MVMAX); } 01258 LDC2 = ConvertType( LDC1, LDC2 ); 01259 OType1C = new SType[LDC1 * N1]; 01260 OType2C = new SType[LDC2 * N2]; 01261 for(j1 = 0, j2 = 0; j1 < LDC1 * N1; j1++, j2++) { 01262 OType1C[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01263 OType2C[j2] = OType1C[j1]; 01264 } 01265 01266 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01267 OType1beta = GetRandom(-SCALARMAX, SCALARMAX); 01268 OType2alpha = OType1alpha; 01269 OType2beta = OType1beta; 01270 if(debug) 01271 { 01272 std::cout << "Test #" << TotalTestCount << " -- SYMM -- " << std::endl; 01273 std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t" 01274 << "UPLO = " << Teuchos::EUploChar[UPLO] << std::endl; 01275 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01276 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01277 std::cout << "OType1beta = " << OType1beta << std::endl; 01278 std::cout << "OType2beta = " << OType2beta << std::endl; 01279 if (Teuchos::ESideChar[SIDE] == 'L') { 01280 PrintMatrix(OType1A, M1, M1, LDA1,"OType1A", matlab); 01281 PrintMatrix(OType2A, M2, M2, LDA2,"OType2A", matlab); 01282 } else { 01283 PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab); 01284 PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab); 01285 } 01286 PrintMatrix(OType1B, M1, N1, LDB1,"OType1B", matlab); 01287 PrintMatrix(OType1C, M1, N1, LDC1,"OType1C_before_operation", matlab); 01288 PrintMatrix(OType2B, M2, N2, LDB2,"OType2B", matlab); 01289 PrintMatrix(OType2C, M2, N2, LDC2,"OType2C_before_operation", matlab); 01290 } 01291 TotalTestCount++; 01292 01293 OType1BLAS.SYMM(SIDE, UPLO, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1, OType1beta, OType1C, LDC1); 01294 OType2BLAS.SYMM(SIDE, UPLO, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2, OType2beta, OType2C, LDC2); 01295 if(debug) 01296 { 01297 PrintMatrix(OType1C, M1, N1, LDC1,"OType1C_after_operation", matlab); 01298 PrintMatrix(OType2C, M2, N2, LDC2,"OType2C_after_operation", matlab); 01299 } 01300 if (CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL)==0) 01301 std::cout << "FAILED TEST!!!!!!" << std::endl; 01302 GoodTestSubcount += CompareMatrices(OType1C, M1, N1, LDC1, OType2C, M2, N2, LDC2, TOL); 01303 01304 delete [] OType1A; 01305 delete [] OType1B; 01306 delete [] OType1C; 01307 delete [] OType2A; 01308 delete [] OType2B; 01309 delete [] OType2C; 01310 } 01311 GoodTestCount += GoodTestSubcount; 01312 if(verbose || debug) std::cout << "SYMM: " << GoodTestSubcount << " of " << SYMMTESTS << " tests were successful." << std::endl; 01313 if(debug) std::cout << std::endl; 01314 //-------------------------------------------------------------------------------- 01315 // End SYMM Tests 01316 //-------------------------------------------------------------------------------- 01317 01318 //-------------------------------------------------------------------------------- 01319 // Begin SYRK Tests 01320 //-------------------------------------------------------------------------------- 01321 GoodTestSubcount = 0; 01322 for(i = 0; i < SYRKTESTS; i++) 01323 { 01324 N1 = GetRandom(MVMIN, MVMAX); 01325 K1 = GetRandom(MVMIN, MVMAX); 01326 while (K1 > N1) { K1 = GetRandom(MVMIN, MVMAX); } 01327 N2 = ConvertType( N1, N2 ); 01328 K2 = ConvertType( K1, K2 ); 01329 01330 UPLO = RandomUPLO(); 01331 TRANS = RandomTRANS(); 01332 #ifdef HAVE_TEUCHOS_COMPLEX 01333 while (TRANS == Teuchos::CONJ_TRANS) { TRANS = RandomTRANS(); } 01334 #endif 01335 01336 LDA1 = GetRandom(MVMIN, MVMAX); 01337 if(Teuchos::ETranspChar[TRANS] == 'N') { 01338 while (LDA1 < N1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01339 LDA2 = ConvertType( LDA1, LDA2 ); 01340 OType1A = new SType[LDA1 * K1]; 01341 OType2A = new SType[LDA2 * K2]; 01342 for(j1 = 0, j2 = 0; j1 < LDA1 * K1; j1++, j2++) { 01343 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01344 OType2A[j2] = OType1A[j1]; 01345 } 01346 } else { 01347 while (LDA1 < K1) { LDA1 = GetRandom(MVMIN, MVMAX); } 01348 LDA2 = ConvertType( LDA1, LDA2 ); 01349 OType1A = new SType[LDA1 * N1]; 01350 OType2A = new SType[LDA2 * N2]; 01351 for(j1 = 0, j2 = 0; j1 < LDA1 * N1; j1++, j2++) { 01352 OType1A[j1] = GetRandom(-SCALARMAX, SCALARMAX); 01353 OType2A[j2] = OType1A[j1]; 01354 } 01355 } 01356 01357 LDC1 = GetRandom(MVMIN, MVMAX); 01358 while (LDC1 < N1) { LDC1 = GetRandom(MVMIN, MVMAX); } 01359 LDC2 = ConvertType( LDC1, LDC2 ); 01360 OType1C = new SType[LDC1 * N1]; 01361 OType2C = new SType[LDC2 * N2]; 01362 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) { 01363 01364 if(Teuchos::EUploChar[UPLO] == 'U') { 01365 // The operator is upper triangular, make sure that the entries are 01366 // only in the upper triangular part of C. 01367 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01368 { 01369 if(k1 <= j1) { 01370 OType1C[j1*LDC1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01371 } else { 01372 OType1C[j1*LDC1+k1] = STypezero; 01373 } 01374 OType2C[j2*LDC2+k2] = OType1C[j1*LDC1+k1]; 01375 } 01376 } 01377 else { 01378 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01379 { 01380 if(k1 >= j1) { 01381 OType1C[j1*LDC1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01382 } else { 01383 OType1C[j1*LDC1+k1] = STypezero; 01384 } 01385 OType2C[j2*LDC2+k2] = OType1C[j1*LDC1+k1]; 01386 } 01387 } 01388 } 01389 01390 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01391 OType1beta = GetRandom(-SCALARMAX, SCALARMAX); 01392 OType2alpha = OType1alpha; 01393 OType2beta = OType1beta; 01394 if(debug) 01395 { 01396 std::cout << "Test #" << TotalTestCount << " -- SYRK -- " << std::endl; 01397 std::cout << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t" 01398 << "TRANS = " << Teuchos::ETranspChar[TRANS] << std::endl; 01399 std::cout << "N1="<<N1 << "\t" << "K1 = " << K1 01400 << "\t" << "LDA1="<<LDA1 << "\t" << "LDC1=" << LDC1 << std::endl; 01401 std::cout << "N2="<<N2 << "\t" << "K2 = " << K2 01402 << "\t" << "LDA2="<<LDA2 << "\t" << "LDC2=" << LDC2 << std::endl; 01403 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01404 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01405 std::cout << "OType1beta = " << OType1beta << std::endl; 01406 std::cout << "OType2beta = " << OType2beta << std::endl; 01407 if (Teuchos::ETranspChar[TRANS] == 'N') { 01408 PrintMatrix(OType1A, N1, K1, LDA1,"OType1A", matlab); 01409 PrintMatrix(OType2A, N2, K2, LDA2,"OType2A", matlab); 01410 } else { 01411 PrintMatrix(OType1A, K1, N1, LDA1, "OType1A", matlab); 01412 PrintMatrix(OType2A, K2, N2, LDA2, "OType2A", matlab); 01413 } 01414 PrintMatrix(OType1C, N1, N1, LDC1,"OType1C_before_operation", matlab); 01415 PrintMatrix(OType2C, N2, N2, LDC2,"OType2C_before_operation", matlab); 01416 } 01417 TotalTestCount++; 01418 01419 OType1BLAS.SYRK(UPLO, TRANS, N1, K1, OType1alpha, OType1A, LDA1, OType1beta, OType1C, LDC1); 01420 OType2BLAS.SYRK(UPLO, TRANS, N2, K2, OType2alpha, OType2A, LDA2, OType2beta, OType2C, LDC2); 01421 if(debug) 01422 { 01423 PrintMatrix(OType1C, N1, N1, LDC1,"OType1C_after_operation", matlab); 01424 PrintMatrix(OType2C, N2, N2, LDC2,"OType2C_after_operation", matlab); 01425 } 01426 if (CompareMatrices(OType1C, N1, N1, LDC1, OType2C, N2, N2, LDC2, TOL)==0) 01427 std::cout << "FAILED TEST!!!!!!" << std::endl; 01428 GoodTestSubcount += CompareMatrices(OType1C, N1, N1, LDC1, OType2C, N2, N2, LDC2, TOL); 01429 01430 delete [] OType1A; 01431 delete [] OType1C; 01432 delete [] OType2A; 01433 delete [] OType2C; 01434 } 01435 GoodTestCount += GoodTestSubcount; 01436 if(verbose || debug) std::cout << "SYRK: " << GoodTestSubcount << " of " << SYRKTESTS << " tests were successful." << std::endl; 01437 if(debug) std::cout << std::endl; 01438 //-------------------------------------------------------------------------------- 01439 // End SYRK Tests 01440 //-------------------------------------------------------------------------------- 01441 01442 //-------------------------------------------------------------------------------- 01443 // Begin TRMM Tests 01444 //-------------------------------------------------------------------------------- 01445 GoodTestSubcount = 0; 01446 for(i = 0; i < TRMMTESTS; i++) 01447 { 01448 M1 = GetRandom(MVMIN, MVMAX); 01449 N1 = GetRandom(MVMIN, MVMAX); 01450 M2 = ConvertType( M1, M2 ); 01451 N2 = ConvertType( N1, N2 ); 01452 01453 LDB1 = GetRandom(MVMIN, MVMAX); 01454 while (LDB1 < M1) { 01455 LDB1 = GetRandom(MVMIN, MVMAX); 01456 } 01457 LDB2 = ConvertType( LDB1, LDB2 ); 01458 01459 OType1B = new SType[LDB1 * N1]; 01460 OType2B = new SType[LDB2 * N2]; 01461 01462 SIDE = RandomSIDE(); 01463 UPLO = RandomUPLO(); 01464 TRANSA = RandomTRANS(); 01465 DIAG = RandomDIAG(); 01466 01467 if(Teuchos::ESideChar[SIDE] == 'L') // The operator is on the left side 01468 { 01469 LDA1 = GetRandom(MVMIN, MVMAX); 01470 while (LDA1 < M1) { 01471 LDA1 = GetRandom(MVMIN, MVMAX); 01472 } 01473 LDA2 = ConvertType( LDA1, LDA2 ); 01474 01475 OType1A = new SType[LDA1 * M1]; 01476 OType2A = new SType[LDA2 * M2]; 01477 01478 for(j1 = 0, j2 = 0; j1 < M1; j1++, j2++) 01479 { 01480 if(Teuchos::EUploChar[UPLO] == 'U') { 01481 // The operator is upper triangular, make sure that the entries are 01482 // only in the upper triangular part of A and the diagonal is non-zero. 01483 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) 01484 { 01485 if(k1 < j1) { 01486 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01487 } else { 01488 OType1A[j1*LDA1+k1] = STypezero; 01489 } 01490 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01491 if(k1 == j1) { 01492 if (Teuchos::EDiagChar[DIAG] == 'N') { 01493 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01494 while (OType1A[j1*LDA1+k1] == STypezero) { 01495 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01496 } 01497 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01498 } else { 01499 OType1A[j1*LDA1+k1] = STypeone; 01500 OType2A[j2*LDA2+k2] = STypeone; 01501 } 01502 } 01503 } 01504 } else { 01505 // The operator is lower triangular, make sure that the entries are 01506 // only in the lower triangular part of A and the diagonal is non-zero. 01507 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) 01508 { 01509 if(k1 > j1) { 01510 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01511 } else { 01512 OType1A[j1*LDA1+k1] = STypezero; 01513 } 01514 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01515 if(k1 == j1) { 01516 if (Teuchos::EDiagChar[DIAG] == 'N') { 01517 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01518 while (OType1A[j1*LDA1+k1] == STypezero) { 01519 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01520 } 01521 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01522 } else { 01523 OType1A[j1*LDA1+k1] = STypeone; 01524 OType2A[j2*LDA2+k2] = STypeone; 01525 } 01526 } 01527 } // end for(k=0 ... 01528 } // end if(UPLO == 'U') ... 01529 } // end for(j=0 ... 01530 } // if(SIDE == 'L') ... 01531 else // The operator is on the right side 01532 { 01533 LDA1 = GetRandom(MVMIN, MVMAX); 01534 while (LDA1 < N1) { 01535 LDA1 = GetRandom(MVMIN, MVMAX); 01536 } 01537 LDA2 = ConvertType( LDA1, LDA2 ); 01538 01539 OType1A = new SType[LDA1 * N1]; 01540 OType2A = new SType[LDA2 * N2]; 01541 01542 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) 01543 { 01544 if(Teuchos::EUploChar[UPLO] == 'U') { 01545 // The operator is upper triangular, make sure that the entries are 01546 // only in the upper triangular part of A and the diagonal is non-zero. 01547 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01548 { 01549 if(k1 < j1) { 01550 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01551 } else { 01552 OType1A[j1*LDA1+k1] = STypezero; 01553 } 01554 OType2A[j1*LDA1+k1] = OType1A[j1*LDA1+k1]; 01555 if(k1 == j1) { 01556 if (Teuchos::EDiagChar[DIAG] == 'N') { 01557 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01558 while (OType1A[j1*LDA1+k1] == STypezero) { 01559 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01560 } 01561 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01562 } else { 01563 OType1A[j1*LDA1+k1] = STypeone; 01564 OType2A[j2*LDA2+k2] = STypeone; 01565 } 01566 } 01567 } 01568 } else { 01569 // The operator is lower triangular, make sure that the entries are 01570 // only in the lower triangular part of A and the diagonal is non-zero. 01571 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01572 { 01573 if(k1 > j1) { 01574 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01575 } else { 01576 OType1A[j1*LDA1+k1] = STypezero; 01577 } 01578 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01579 if(k1 == j1) { 01580 if (Teuchos::EDiagChar[DIAG] == 'N') { 01581 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01582 while (OType1A[j1*LDA1+k1] == STypezero) { 01583 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01584 } 01585 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01586 } else { 01587 OType1A[j1*LDA1+k1] = STypeone; 01588 OType2A[j2*LDA2+k2] = STypeone; 01589 } 01590 } 01591 } // end for(k=0 ... 01592 } // end if(UPLO == 'U') ... 01593 } // end for(j=0 ... 01594 } // end if(SIDE == 'L') ... 01595 01596 // Fill in the right hand side block B. 01597 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) { 01598 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) { 01599 OType1B[j1*LDB1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01600 OType2B[j2*LDB2+k2] = OType1B[j1*LDB1+k1]; 01601 } 01602 } 01603 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01604 OType2alpha = OType1alpha; 01605 if(debug) 01606 { 01607 std::cout << "Test #" << TotalTestCount << " -- TRMM -- " << std::endl; 01608 std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t" 01609 << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t" 01610 << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t" 01611 << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl; 01612 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01613 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01614 if(Teuchos::ESideChar[SIDE] == 'L') { 01615 PrintMatrix(OType1A, M1, M1, LDA1, "OType1A", matlab); 01616 PrintMatrix(OType2A, M2, M2, LDA2, "OType2A", matlab); 01617 } else { 01618 PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab); 01619 PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab); 01620 } 01621 PrintMatrix(OType1B, M1, N1, LDB1,"OType1B_before_operation", matlab); 01622 PrintMatrix(OType2B, M2, N2, LDB2,"OType2B_before_operation", matlab); 01623 } 01624 TotalTestCount++; 01625 OType1BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1); 01626 OType2BLAS.TRMM(SIDE, UPLO, TRANSA, DIAG, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2); 01627 if(debug) 01628 { 01629 PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_after_operation", matlab); 01630 PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_after_operation", matlab); 01631 } 01632 if (CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL)==0) 01633 std::cout << "FAILED TEST!!!!!!" << std::endl; 01634 GoodTestSubcount += CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL); 01635 delete [] OType1A; 01636 delete [] OType1B; 01637 delete [] OType2A; 01638 delete [] OType2B; 01639 } 01640 GoodTestCount += GoodTestSubcount; 01641 if(verbose || debug) std::cout << "TRMM: " << GoodTestSubcount << " of " << TRMMTESTS << " tests were successful." << std::endl; 01642 if(debug) std::cout << std::endl; 01643 //-------------------------------------------------------------------------------- 01644 // End TRMM Tests 01645 //-------------------------------------------------------------------------------- 01646 01647 //-------------------------------------------------------------------------------- 01648 // Begin TRSM Tests 01649 //-------------------------------------------------------------------------------- 01650 GoodTestSubcount = 0; 01651 for(i = 0; i < TRSMTESTS; i++) 01652 { 01653 M1 = GetRandom(MVMIN, MVMAX); 01654 N1 = GetRandom(MVMIN, MVMAX); 01655 M2 = ConvertType( M1, M2 ); 01656 N2 = ConvertType( N1, N2 ); 01657 01658 LDB1 = GetRandom(MVMIN, MVMAX); 01659 while (LDB1 < M1) { 01660 LDB1 = GetRandom(MVMIN, MVMAX); 01661 } 01662 LDB2 = ConvertType( LDB1, LDB2 ); 01663 01664 OType1B = new SType[LDB1 * N1]; 01665 OType2B = new SType[LDB2 * N2]; 01666 01667 SIDE = RandomSIDE(); 01668 UPLO = RandomUPLO(); 01669 TRANSA = RandomTRANS(); 01670 // Since the entries are integers, we don't want to use the unit diagonal feature, 01671 // this creates ill-conditioned, nearly-singular matrices. 01672 //DIAG = RandomDIAG(); 01673 DIAG = Teuchos::NON_UNIT_DIAG; 01674 01675 if(Teuchos::ESideChar[SIDE] == 'L') // The operator is on the left side 01676 { 01677 LDA1 = GetRandom(MVMIN, MVMAX); 01678 while (LDA1 < M1) { 01679 LDA1 = GetRandom(MVMIN, MVMAX); 01680 } 01681 LDA2 = ConvertType( LDA1, LDA2 ); 01682 01683 OType1A = new SType[LDA1 * M1]; 01684 OType2A = new SType[LDA2 * M2]; 01685 01686 for(j1 = 0, j2 = 0; j1 < M1; j1++, j2++) 01687 { 01688 if(Teuchos::EUploChar[UPLO] == 'U') { 01689 // The operator is upper triangular, make sure that the entries are 01690 // only in the upper triangular part of A and the diagonal is non-zero. 01691 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) 01692 { 01693 if(k1 < j1) { 01694 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01695 } else { 01696 OType1A[j1*LDA1+k1] = STypezero; 01697 } 01698 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01699 if(k1 == j1) { 01700 if (Teuchos::EDiagChar[DIAG] == 'N') { 01701 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01702 while (OType1A[j1*LDA1+k1] == STypezero) { 01703 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01704 } 01705 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01706 } else { 01707 OType1A[j1*LDA1+k1] = STypeone; 01708 OType2A[j2*LDA2+k2] = STypeone; 01709 } 01710 } 01711 } 01712 } else { 01713 // The operator is lower triangular, make sure that the entries are 01714 // only in the lower triangular part of A and the diagonal is non-zero. 01715 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) 01716 { 01717 if(k1 > j1) { 01718 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01719 } else { 01720 OType1A[j1*LDA1+k1] = STypezero; 01721 } 01722 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01723 if(k1 == j1) { 01724 if (Teuchos::EDiagChar[DIAG] == 'N') { 01725 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01726 while (OType1A[j1*LDA1+k1] == STypezero) { 01727 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01728 } 01729 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01730 } else { 01731 OType1A[j1*LDA1+k1] = STypeone; 01732 OType2A[j2*LDA2+k2] = STypeone; 01733 } 01734 } 01735 } // end for(k=0 ... 01736 } // end if(UPLO == 'U') ... 01737 } // end for(j=0 ... 01738 } // if(SIDE == 'L') ... 01739 else // The operator is on the right side 01740 { 01741 LDA1 = GetRandom(MVMIN, MVMAX); 01742 while (LDA1 < N1) { 01743 LDA1 = GetRandom(MVMIN, MVMAX); 01744 } 01745 LDA2 = ConvertType( LDA1, LDA2 ); 01746 01747 OType1A = new SType[LDA1 * N1]; 01748 OType2A = new SType[LDA2 * N2]; 01749 01750 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) 01751 { 01752 if(Teuchos::EUploChar[UPLO] == 'U') { 01753 // The operator is upper triangular, make sure that the entries are 01754 // only in the upper triangular part of A and the diagonal is non-zero. 01755 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01756 { 01757 if(k1 < j1) { 01758 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01759 } else { 01760 OType1A[j1*LDA1+k1] = STypezero; 01761 } 01762 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01763 if(k1 == j1) { 01764 if (Teuchos::EDiagChar[DIAG] == 'N') { 01765 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01766 while (OType1A[j1*LDA1+k1] == STypezero) { 01767 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01768 } 01769 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01770 } else { 01771 OType1A[j1*LDA1+k1] = STypeone; 01772 OType2A[j2*LDA2+k2] = STypeone; 01773 } 01774 } 01775 } 01776 } else { 01777 // The operator is lower triangular, make sure that the entries are 01778 // only in the lower triangular part of A and the diagonal is non-zero. 01779 for(k1 = 0, k2 = 0; k1 < N1; k1++, k2++) 01780 { 01781 if(k1 > j1) { 01782 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01783 } else { 01784 OType1A[j1*LDA1+k1] = STypezero; 01785 } 01786 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01787 if(k1 == j1) { 01788 if (Teuchos::EDiagChar[DIAG] == 'N') { 01789 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01790 while (OType1A[j1*LDA1+k1] == STypezero) { 01791 OType1A[j1*LDA1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01792 } 01793 OType2A[j2*LDA2+k2] = OType1A[j1*LDA1+k1]; 01794 } else { 01795 OType1A[j1*LDA1+k1] = STypeone; 01796 OType2A[j2*LDA2+k2] = STypeone; 01797 } 01798 } 01799 } // end for(k=0 ... 01800 } // end if(UPLO == 'U') ... 01801 } // end for(j=0 ... 01802 } // end if(SIDE == 'L') ... 01803 01804 // Fill in the right hand side block B. 01805 for(j1 = 0, j2 = 0; j1 < N1; j1++, j2++) 01806 { 01807 for(k1 = 0, k2 = 0; k1 < M1; k1++, k2++) 01808 { 01809 OType1B[j1*LDB1+k1] = GetRandom(-SCALARMAX, SCALARMAX); 01810 OType2B[j2*LDB2+k2] = OType1B[j1*LDB1+k1]; 01811 } 01812 } 01813 01814 OType1alpha = GetRandom(-SCALARMAX, SCALARMAX); 01815 OType2alpha = OType1alpha; 01816 01817 if(debug) 01818 { 01819 std::cout << "Test #" << TotalTestCount << " -- TRSM -- " << std::endl; 01820 std::cout << "SIDE = " << Teuchos::ESideChar[SIDE] << "\t" 01821 << "UPLO = " << Teuchos::EUploChar[UPLO] << "\t" 01822 << "TRANSA = " << Teuchos::ETranspChar[TRANSA] << "\t" 01823 << "DIAG = " << Teuchos::EDiagChar[DIAG] << std::endl; 01824 std::cout << "M1="<<M1 << "\t" << "N1="<<N1 << "\t" << "LDA1="<<LDA1 << "\t" << "LDB1="<<LDB1 << std::endl; 01825 std::cout << "M2="<<M2 << "\t" << "N2="<<N2 << "\t" << "LDA2="<<LDA2 << "\t" << "LDB2="<<LDB2 << std::endl; 01826 std::cout << "OType1alpha = " << OType1alpha << std::endl; 01827 std::cout << "OType2alpha = " << OType2alpha << std::endl; 01828 if (Teuchos::ESideChar[SIDE] == 'L') { 01829 PrintMatrix(OType1A, M1, M1, LDA1, "OType1A", matlab); 01830 PrintMatrix(OType2A, M2, M2, LDA2, "OType2A", matlab); 01831 } else { 01832 PrintMatrix(OType1A, N1, N1, LDA1, "OType1A", matlab); 01833 PrintMatrix(OType2A, N2, N2, LDA2, "OType2A", matlab); 01834 } 01835 PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_before_operation", matlab); 01836 PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_before_operation", matlab); 01837 } 01838 TotalTestCount++; 01839 01840 OType1BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M1, N1, OType1alpha, OType1A, LDA1, OType1B, LDB1); 01841 OType2BLAS.TRSM(SIDE, UPLO, TRANSA, DIAG, M2, N2, OType2alpha, OType2A, LDA2, OType2B, LDB2); 01842 01843 if(debug) 01844 { 01845 PrintMatrix(OType1B, M1, N1, LDB1, "OType1B_after_operation", matlab); 01846 PrintMatrix(OType2B, M2, N2, LDB2, "OType2B_after_operation", matlab); 01847 } 01848 01849 if (CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL)==0) 01850 std::cout << "FAILED TEST!!!!!!" << std::endl; 01851 GoodTestSubcount += CompareMatrices(OType1B, M1, N1, LDB1, OType2B, M2, N2, LDB2, TOL); 01852 01853 delete [] OType1A; 01854 delete [] OType1B; 01855 delete [] OType2A; 01856 delete [] OType2B; 01857 } 01858 GoodTestCount += GoodTestSubcount; 01859 if(verbose || debug) std::cout << "TRSM: " << GoodTestSubcount << " of " << TRSMTESTS << " tests were successful." << std::endl; 01860 if(debug) std::cout << std::endl; 01861 //-------------------------------------------------------------------------------- 01862 // End TRSM Tests 01863 //-------------------------------------------------------------------------------- 01864 01865 if((((TotalTestCount - 1) - GoodTestCount) != 0) || (verbose) || (debug)) 01866 { 01867 std::cout << GoodTestCount << " of " << (TotalTestCount - 1) << " total tests were successful." << std::endl; 01868 } 01869 01870 if ((TotalTestCount-1) == GoodTestCount) { 01871 std::cout << "End Result: TEST PASSED" << std::endl; 01872 return 0; 01873 } 01874 01875 std::cout << "End Result: TEST FAILED" << std::endl; 01876 return (TotalTestCount-GoodTestCount-1); 01877 } 01878 01879 template<typename TYPE> 01880 TYPE GetRandom(TYPE Low, TYPE High) 01881 { 01882 return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low); 01883 } 01884 01885 template<typename T> 01886 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High) 01887 { 01888 T lowMag = Low.real(); 01889 T highMag = High.real(); 01890 T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag; 01891 T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag; 01892 return std::complex<T>( real, imag ); 01893 } 01894 01895 template<> 01896 int GetRandom(int Low, int High) 01897 { 01898 return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low); 01899 } 01900 01901 template<> 01902 double GetRandom(double Low, double High) 01903 { 01904 return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random()); 01905 } 01906 01907 template<typename TYPE, typename OTYPE> 01908 void PrintVector(TYPE* Vector, OTYPE Size, std::string Name, bool Matlab) 01909 { 01910 std::cout << Name << " =" << std::endl; 01911 OTYPE i; 01912 if(Matlab) std::cout << "["; 01913 for(i = 0; i < Size; i++) 01914 { 01915 std::cout << Vector[i] << " "; 01916 } 01917 if(Matlab) std::cout << "]"; 01918 if(!Matlab) 01919 { 01920 std::cout << std::endl << std::endl; 01921 } 01922 else 01923 { 01924 std::cout << ";" << std::endl; 01925 } 01926 } 01927 01928 template<typename TYPE, typename OTYPE> 01929 void PrintMatrix(TYPE* Matrix, OTYPE Rows, OTYPE Columns, OTYPE LDM, std::string Name, bool Matlab) 01930 { 01931 if(!Matlab) 01932 { 01933 std::cout << Name << " =" << std::endl; 01934 OTYPE i, j; 01935 for(i = 0; i < Rows; i++) 01936 { 01937 for(j = 0; j < Columns; j++) 01938 { 01939 std::cout << Matrix[i + (j * LDM)] << " "; 01940 } 01941 std::cout << std::endl; 01942 } 01943 std::cout << std::endl; 01944 } 01945 else 01946 { 01947 std::cout << Name << " = "; 01948 OTYPE i, j; 01949 std::cout << "["; 01950 for(i = 0; i < Rows; i++) 01951 { 01952 std::cout << "["; 01953 for(j = 0; j < Columns; j++) 01954 { 01955 std::cout << Matrix[i + (j * LDM)] << " "; 01956 } 01957 std::cout << "];"; 01958 } 01959 std::cout << "];" << std::endl; 01960 } 01961 } 01962 01963 template<typename TYPE> 01964 bool CompareScalars(TYPE Scalar1, TYPE Scalar2, typename ScalarTraits<TYPE>::magnitudeType Tolerance) 01965 { 01966 typename ScalarTraits<TYPE>::magnitudeType temp = ScalarTraits<TYPE>::magnitude(Scalar2); 01967 typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<TYPE>::magnitude(Scalar1 - Scalar2); 01968 if (temp != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero()) { 01969 temp2 /= temp; 01970 } 01971 return( temp2 < Tolerance ); 01972 } 01973 01974 01975 /* Function: CompareVectors 01976 Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2 01977 */ 01978 template<typename TYPE, typename OTYPE1, typename OTYPE2> 01979 bool CompareVectors(TYPE* Vector1, OTYPE1 Size1, TYPE* Vector2, OTYPE2 Size2, typename ScalarTraits<TYPE>::magnitudeType Tolerance) 01980 { 01981 TYPE temp = ScalarTraits<TYPE>::zero(); 01982 typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 01983 typename ScalarTraits<TYPE>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 01984 typename ScalarTraits<TYPE>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 01985 typename ScalarTraits<TYPE>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 01986 OTYPE1 i1; 01987 OTYPE2 i2; 01988 for(i1 = 0, i2 = 0; i1 < Size1; i1++, i2++) 01989 { 01990 sum2 += ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(Vector2[i2])*Vector2[i2]); 01991 temp = Vector1[i1] - Vector2[i2]; 01992 sum += ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(temp)*temp); 01993 } 01994 temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum2); 01995 if (temp2 != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero()) 01996 temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum)/temp2; 01997 else 01998 temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum); 01999 if (temp3 > Tolerance ) 02000 return false; 02001 else 02002 return true; 02003 } 02004 02005 /* Function: CompareMatrices 02006 Purpose: Compares the difference between two matrices using relative frobenius-norms, i.e. ||M_1-M_2||_F/||M_2||_F 02007 */ 02008 template<typename TYPE, typename OTYPE1, typename OTYPE2> 02009 bool CompareMatrices(TYPE* Matrix1, OTYPE1 Rows1, OTYPE1 Columns1, OTYPE1 LDM1, 02010 TYPE* Matrix2, OTYPE2 Rows2, OTYPE2 Columns2, OTYPE2 LDM2, 02011 typename ScalarTraits<TYPE>::magnitudeType Tolerance) 02012 { 02013 TYPE temp = ScalarTraits<TYPE>::zero(); 02014 typename ScalarTraits<TYPE>::magnitudeType temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 02015 typename ScalarTraits<TYPE>::magnitudeType temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 02016 typename ScalarTraits<TYPE>::magnitudeType sum = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 02017 typename ScalarTraits<TYPE>::magnitudeType sum2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero(); 02018 OTYPE1 i1, j1; 02019 OTYPE2 i2, j2; 02020 for(j1 = 0, j2 = 0; j1 < Columns1; j1++, j2++) 02021 { 02022 for(i1 = 0, i2 = 0; i1 < Rows1; i1++, i2++) 02023 { 02024 sum2 = ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(Matrix2[j2*LDM2 + i2])*Matrix2[j2*LDM2 + i2]); 02025 temp = Matrix1[j1*LDM1 + i1] - Matrix2[j2*LDM2 + i2]; 02026 sum = ScalarTraits<TYPE>::magnitude(ScalarTraits<TYPE>::conjugate(temp)*temp); 02027 } 02028 } 02029 temp2 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum2); 02030 if (temp2 != ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::zero()) 02031 temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum)/temp2; 02032 else 02033 temp3 = ScalarTraits<typename ScalarTraits<TYPE>::magnitudeType>::squareroot(sum); 02034 if (temp3 > Tolerance) 02035 return false; 02036 else 02037 return true; 02038 } 02039 02040 02041 Teuchos::ESide RandomSIDE() 02042 { 02043 Teuchos::ESide result; 02044 int r = GetRandom(1, 2); 02045 if(r == 1) 02046 { 02047 result = Teuchos::LEFT_SIDE; 02048 } 02049 else 02050 { 02051 result = Teuchos::RIGHT_SIDE; 02052 } 02053 return result; 02054 } 02055 02056 Teuchos::EUplo RandomUPLO() 02057 { 02058 Teuchos::EUplo result; 02059 int r = GetRandom(1, 2); 02060 if(r == 1) 02061 { 02062 result = Teuchos::UPPER_TRI; 02063 } 02064 else 02065 { 02066 result = Teuchos::LOWER_TRI; 02067 } 02068 return result; 02069 } 02070 02071 Teuchos::ETransp RandomTRANS() 02072 { 02073 Teuchos::ETransp result; 02074 int r = GetRandom(1, 4); 02075 if(r == 1 || r == 2) 02076 { 02077 result = Teuchos::NO_TRANS; 02078 } 02079 else if(r == 3) 02080 { 02081 result = Teuchos::TRANS; 02082 } 02083 else 02084 { 02085 result = Teuchos::CONJ_TRANS; 02086 } 02087 return result; 02088 } 02089 02090 Teuchos::EDiag RandomDIAG() 02091 { 02092 Teuchos::EDiag result; 02093 int r = GetRandom(1, 2); 02094 if(r == 1) 02095 { 02096 result = Teuchos::NON_UNIT_DIAG; 02097 } 02098 else 02099 { 02100 result = Teuchos::UNIT_DIAG; 02101 } 02102 return result; 02103 } 02104
1.7.6.1
