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