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