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