65 #ifdef HAVE_TEUCHOS_ARPREC 
   66 #define SType1     mp_real 
   67 #elif defined(HAVE_TEUCHOS_QD) 
   68 #define SType1     dd_real 
   69 #elif defined(HAVE_TEUCHOS_GNU_MP) 
   70 #define SType1     mpf_class 
  108 template<
typename TYPE>
 
  119 template<
typename TYPE>
 
  120 void PrintVector(TYPE* Vector, 
OType Size, std::string Name, 
bool Matlab = 0);
 
  122 template<
typename TYPE>
 
  125 template<
typename TYPE1, 
typename TYPE2>
 
  126 bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, TYPE2 Tolerance );
 
  128 template<
typename TYPE1, 
typename TYPE2>
 
  131 template<
typename TYPE1, 
typename TYPE2>
 
  136 template<
typename TYPE1, 
typename TYPE2>
 
  139 #ifdef HAVE_TEUCHOS_ARPREC 
  144 #ifdef HAVE_TEUCHOS_QD 
  149 #ifdef HAVE_TEUCHOS_GNU_MP 
  160 int main(
int argc, 
char *argv[])
 
  162 #ifdef HAVE_TEUCHOS_ARPREC 
  171   bool InvalidCmdLineArgs = 0;
 
  173   for(i = 1; i < argc; i++)
 
  175     if(argv[i][0] == 
'-')
 
  186             InvalidCmdLineArgs = 1;
 
  196             InvalidCmdLineArgs = 1;
 
  206             InvalidCmdLineArgs = 1;
 
  210           InvalidCmdLineArgs = 1;
 
  219   if(InvalidCmdLineArgs || (argc > 4))
 
  221     std::cout << 
"Invalid command line arguments detected. Use the following flags:" << std::endl
 
  222       << 
"\t -v enables verbose mode (reports number of failed/successful tests)" << std::endl
 
  223       << 
"\t -d enables debug mode (same as verbose with output of each test, not recommended for large numbers of tests)" << std::endl
 
  224       << 
"\t -m enables matlab-style output; only has an effect if debug mode is enabled" << std::endl;
 
  237   SType1 SType1alpha, SType1beta;
 
  243   SType2 SType2alpha, SType2beta;
 
  244   SType1 SType1ASUMresult, SType1DOTresult, SType1NRM2result, SType1SINresult, SType1COSresult;
 
  245   SType2 SType2ASUMresult, SType2DOTresult, SType2NRM2result, SType2SINresult, SType2COSresult;
 
  247   OType SType1IAMAXresult;
 
  248   OType SType2IAMAXresult;
 
  249   OType TotalTestCount = 1, GoodTestSubcount, GoodTestCount = 0, M, M2, N, N2, K, P, LDA, LDB, LDC, Mx, My;
 
  257   std::srand(time(NULL));
 
  264   GoodTestSubcount = 0;
 
  272     SType2COSresult = 
ConvertType(SType1COSresult, convertTo);
 
  274     SType2SINresult = 
ConvertType(SType1SINresult, convertTo);
 
  278       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  279       std::cout << 
"SType1alpha = "  << SType1alpha << std::endl;
 
  280       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
  281       std::cout << 
"SType1beta = "  << SType1beta << std::endl;
 
  282       std::cout << 
"SType2beta = " << SType2beta << std::endl;
 
  285     SType1BLAS.
ROTG(&SType1alpha, &SType1beta, &SType1COSresult, &SType1SINresult);
 
  286     SType2BLAS.
ROTG(&SType2alpha, &SType2beta, &SType2COSresult, &SType2SINresult);
 
  289       std::cout << 
"SType1 ROTG COS result: " << SType1COSresult << std::endl;
 
  290       std::cout << 
"SType2 ROTG COS result: " << SType2COSresult << std::endl;
 
  291       std::cout << 
"SType1 ROTG SIN result: " << SType1SINresult << std::endl;
 
  292       std::cout << 
"SType2 ROTG SIN result: " << SType2SINresult << std::endl;
 
  294     GoodTestSubcount += ( 
CompareScalars(SType1COSresult, SType2COSresult, TOL) &&
 
  297   GoodTestCount += GoodTestSubcount;
 
  298   if(verbose || debug) std::cout << 
"ROTG: " << GoodTestSubcount << 
" of " << ROTGTESTS << 
" tests were successful." << std::endl;
 
  299   if(debug) std::cout << std::endl;
 
  309   GoodTestSubcount = 0;
 
  314     if (incx == 0) incx = 1;
 
  315     if (incy == 0) incy = 1;
 
  317     Mx = M*std::abs(incx);
 
  318     My = M*std::abs(incy);
 
  319     if (Mx == 0) { Mx = 1; }
 
  320     if (My == 0) { My = 1; }
 
  325     for(j = 0; j < Mx; j++)
 
  330     for(j = 0; j < My; j++)
 
  341       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  342       std::cout << 
"c1 = "  << c1 << 
", s1 = " << s1 << std::endl;
 
  343       std::cout << 
"c2 = " << c2 << 
", s2 = " << s2 << std::endl;
 
  344       std::cout << 
"incx = " << incx << 
", incy = " << incy << std::endl;
 
  346       PrintVector(SType1y, My, 
"SType1y_before_operation", matlab);
 
  348       PrintVector(SType2y, My, 
"SType2y_before_operation",  matlab);
 
  351     SType1BLAS.
ROT(M, SType1x, incx, SType1y, incy, &c1, &s1);
 
  352     SType2BLAS.
ROT(M, SType2x, incx, SType2y, incy, &c2, &s2);
 
  355       PrintVector(SType1y, My, 
"SType1y_after_operation", matlab);
 
  356       PrintVector(SType2y, My, 
"SType2y_after_operation", matlab);
 
  358     GoodTestSubcount += ( 
CompareVectors(SType1x, SType2x, Mx, TOL) &&
 
  365   GoodTestCount += GoodTestSubcount;
 
  366   if(verbose || debug) std::cout << 
"ROT: " << GoodTestSubcount << 
" of " << ROTTESTS << 
" tests were successful." << std::endl;
 
  367   if(debug) std::cout << std::endl;
 
  375   GoodTestSubcount = 0;
 
  384     for(j = 0; j < M2; j++)
 
  391       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  396     SType1ASUMresult = SType1BLAS.
ASUM(M, SType1x, incx);
 
  397     SType2ASUMresult = SType2BLAS.
ASUM(M, SType2x, incx);
 
  400       std::cout << 
"SType1 ASUM result: " << SType1ASUMresult << std::endl;
 
  401       std::cout << 
"SType2 ASUM result: " << SType2ASUMresult << std::endl;
 
  403     GoodTestSubcount += 
CompareScalars(SType1ASUMresult, SType2ASUMresult, TOL);
 
  407   GoodTestCount += GoodTestSubcount;
 
  408   if(verbose || debug) std::cout << 
"ASUM: " << GoodTestSubcount << 
" of " << ASUMTESTS << 
" tests were successful." << std::endl;
 
  409   if(debug) std::cout << std::endl;
 
  418   GoodTestSubcount = 0;
 
  424     Mx = M*std::abs(incx);
 
  425     My = M*std::abs(incy);
 
  426     if (Mx == 0) { Mx = 1; }
 
  427     if (My == 0) { My = 1; }
 
  432     for(j = 0; j < Mx; j++)
 
  437     for(j = 0; j < My; j++)
 
  446       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  447       std::cout << 
"SType1alpha = "  << SType1alpha << std::endl;
 
  448       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
  450       PrintVector(SType1y, My, 
"SType1y_before_operation", matlab);
 
  452       PrintVector(SType2y, My, 
"SType2y_before_operation",  matlab);
 
  455     SType1BLAS.
AXPY(M, SType1alpha, SType1x, incx, SType1y, incy);
 
  456     SType2BLAS.
AXPY(M, SType2alpha, SType2x, incx, SType2y, incy);
 
  459       PrintVector(SType1y, My, 
"SType1y_after_operation", matlab);
 
  460       PrintVector(SType2y, My, 
"SType2y_after_operation", matlab);
 
  468   GoodTestCount += GoodTestSubcount;
 
  469   if(verbose || debug) std::cout << 
"AXPY: " << GoodTestSubcount << 
" of " << AXPYTESTS << 
" tests were successful." << std::endl;
 
  470   if(debug) std::cout << std::endl;
 
  478   GoodTestSubcount = 0;
 
  484     Mx = M*std::abs(incx);
 
  485     My = M*std::abs(incy);
 
  486     if (Mx == 0) { Mx = 1; }
 
  487     if (My == 0) { My = 1; }
 
  492     for(j = 0; j < Mx; j++)
 
  497     for(j = 0; j < My; j++)
 
  504       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  506       PrintVector(SType1y, My, 
"SType1y_before_operation", matlab);
 
  508       PrintVector(SType2y, My, 
"SType2y_before_operation", matlab);
 
  511     SType1BLAS.
COPY(M, SType1x, incx, SType1y, incy);
 
  512     SType2BLAS.
COPY(M, SType2x, incx, SType2y, incy);
 
  515       PrintVector(SType1y, My, 
"SType1y_after_operation", matlab);
 
  516       PrintVector(SType2y, My, 
"SType2y_after_operation", matlab);
 
  524   GoodTestCount += GoodTestSubcount; 
if(verbose || debug) std::cout << 
"COPY: " << GoodTestSubcount << 
" of " << COPYTESTS << 
" tests were successful." << std::endl;
 
  525   if(debug) std::cout << std::endl;
 
  533   GoodTestSubcount = 0;
 
  539     Mx = M*std::abs(incx);
 
  540     My = M*std::abs(incy);
 
  541     if (Mx == 0) { Mx = 1; }
 
  542     if (My == 0) { My = 1; }
 
  547     for(j = 0; j < Mx; j++)
 
  552     for(j = 0; j < My; j++)
 
  559       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  566     SType1DOTresult = SType1BLAS.
DOT(M, SType1x, incx, SType1y, incy);
 
  567     SType2DOTresult = SType2BLAS.
DOT(M, SType2x, incx, SType2y, incy);
 
  570       std::cout << 
"SType1 DOT result: " << SType1DOTresult << std::endl;
 
  571       std::cout << 
"SType2 DOT result: " << SType2DOTresult << std::endl;
 
  573     GoodTestSubcount += 
CompareScalars(SType1DOTresult, SType2DOTresult, TOL);
 
  579   GoodTestCount += GoodTestSubcount;
 
  580   if(verbose || debug) std::cout << 
"DOT: " << GoodTestSubcount << 
" of " << DOTTESTS << 
" tests were successful." << std::endl;
 
  581   if(debug) std::cout << std::endl;
 
  589   GoodTestSubcount = 0;
 
  597     for(j = 0; j < M2; j++)
 
  604       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  609     SType1NRM2result = SType1BLAS.
NRM2(M, SType1x, incx);
 
  610     SType2NRM2result = SType2BLAS.
NRM2(M, SType2x, incx);
 
  613       std::cout << 
"SType1 NRM2 result: " << SType1NRM2result << std::endl;
 
  614       std::cout << 
"SType2 NRM2 result: " << SType2NRM2result << std::endl;
 
  616     GoodTestSubcount += 
CompareScalars(SType1NRM2result, SType2NRM2result, TOL);
 
  620   GoodTestCount += GoodTestSubcount; 
if(verbose || debug) std::cout << 
"NRM2: " << GoodTestSubcount << 
" of " << NRM2TESTS << 
" tests were successful." << std::endl;
 
  621   if(debug) std::cout << std::endl;
 
  629   GoodTestSubcount = 0;
 
  640     for(j = 0; j < M2; j++)
 
  649       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  650       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
  651       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
  652       PrintVector(SType1x, M2, 
"SType1x_before_operation", matlab);
 
  653       PrintVector(SType2x, M2, 
"SType2x_before_operation", matlab);
 
  656     SType1BLAS.
SCAL(M, SType1alpha, SType1x, incx);
 
  657     SType2BLAS.
SCAL(M, SType2alpha, SType2x, incx);
 
  660       PrintVector(SType1x, M2, 
"SType1x_after_operation", matlab);
 
  661       PrintVector(SType2x, M2, 
"SType2x_after_operation", matlab);
 
  667   GoodTestCount += GoodTestSubcount;
 
  668   if(verbose || debug) std::cout << 
"SCAL: " << GoodTestSubcount << 
" of " << SCALTESTS << 
" tests were successful." << std::endl;
 
  669   if(debug) std::cout << std::endl;
 
  677   GoodTestSubcount = 0;
 
  685     for(j = 0; j < M2; j++)
 
  692       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  697     SType1IAMAXresult = SType1BLAS.
IAMAX(M, SType1x, incx);
 
  698     SType2IAMAXresult = SType2BLAS.
IAMAX(M, SType2x, incx);
 
  701       std::cout << 
"SType1 IAMAX result: " << SType1IAMAXresult << std::endl;
 
  702       std::cout << 
"SType2 IAMAX result: " << SType2IAMAXresult << std::endl;
 
  704     GoodTestSubcount += (SType1IAMAXresult == SType2IAMAXresult);
 
  708   GoodTestCount += GoodTestSubcount;
 
  709   if(verbose || debug) std::cout << 
"IAMAX: " << GoodTestSubcount << 
" of " << IAMAXTESTS << 
" tests were successful." << std::endl;
 
  710   if(debug) std::cout << std::endl;
 
  720   GoodTestSubcount = 0;
 
  738       M2 = M*std::abs(incy);
 
  739       N2 = N*std::abs(incx);
 
  741       M2 = N*std::abs(incy);
 
  742       N2 = M*std::abs(incx);
 
  755     SType1A = 
new SType1[LDA * N];
 
  758     SType2A = 
new SType2[LDA * N];
 
  762     for(j = 0; j < LDA * N; j++)
 
  767     for(j = 0; j < N2; j++)
 
  772     for(j = 0; j < M2; j++)
 
  779       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  780       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
  781       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
  782       std::cout << 
"SType1beta = " << SType1beta << std::endl;
 
  783       std::cout << 
"SType2beta = " << SType2beta << std::endl;
 
  784       PrintMatrix(SType1A, M, N, LDA, 
"SType1A", matlab);
 
  786       PrintVector(SType1y, M2, 
"SType1y_before_operation", matlab);
 
  787       PrintMatrix(SType2A, M, N, LDA, 
"SType2A", matlab);
 
  789       PrintVector(SType2y, M2, 
"SType2y_before_operation", matlab);
 
  792     SType1BLAS.
GEMV(TRANS, M, N, SType1alpha, SType1A, LDA, SType1x, incx, SType1beta, SType1y, incy);
 
  793     SType2BLAS.
GEMV(TRANS, M, N, SType2alpha, SType2A, LDA, SType2x, incx, SType2beta, SType2y, incy);
 
  796       PrintVector(SType1y, M2, 
"SType1y_after_operation", matlab);
 
  797       PrintVector(SType2y, M2, 
"SType2y_after_operation", matlab);
 
  807   GoodTestCount += GoodTestSubcount;
 
  808   if(verbose || debug) std::cout << 
"GEMV: " << GoodTestSubcount << 
" of " << GEMVTESTS << 
" tests were successful." << std::endl;
 
  809   if(debug) std::cout << std::endl;
 
  817   GoodTestSubcount = 0;
 
  832     N2 = N*std::abs(incx);
 
  836     for(j = 0; j < N2; j++)
 
  846     SType1A = 
new SType1[LDA * N];
 
  847     SType2A = 
new SType2[LDA * N];
 
  849     for(j = 0; j < N; j++)
 
  854         for(k = 0; k < N; k++)
 
  859             SType1A[j*LDA+k] = SType1zero;
 
  861           SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
  865               while (SType1A[j*LDA+k] == SType1zero) {
 
  868               SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
  870               SType1A[j*LDA+k] = SType1one;
 
  871               SType2A[j*LDA+k] = SType2one;
 
  878         for(k = 0; k < N; k++)
 
  883             SType1A[j*LDA+k] = SType1zero;
 
  885           SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
  889               while (SType1A[j*LDA+k] == SType1zero) {
 
  892               SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
  894               SType1A[j*LDA+k] = SType1one;
 
  895               SType2A[j*LDA+k] = SType2one;
 
  904       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  906       PrintVector(SType1x, N2, 
"SType1x_before_operation", matlab);
 
  907       PrintMatrix(SType2A, N, N, LDA, 
"SType2A", matlab);
 
  908       PrintVector(SType2x, N2, 
"SType2x_before_operation", matlab);
 
  911     SType1BLAS.
TRMV(UPLO, TRANSA, DIAG, N, SType1A, LDA, SType1x, incx);
 
  912     SType2BLAS.
TRMV(UPLO, TRANSA, DIAG, N, SType2A, LDA, SType2x, incx);
 
  915       PrintVector(SType1x, N2, 
"SType1x_after_operation", matlab);
 
  916       PrintVector(SType2x, N2, 
"SType2x_after_operation", matlab);
 
  924   GoodTestCount += GoodTestSubcount;
 
  925   if(verbose || debug) std::cout << 
"TRMV: " << GoodTestSubcount << 
" of " << TRMVTESTS << 
" tests were successful." << std::endl;
 
  926   if(debug) std::cout << std::endl;
 
  934   GoodTestSubcount = 0;
 
  948     M2 = M*std::abs(incx);
 
  949     N2 = N*std::abs(incy);
 
  956     SType1A = 
new SType1[LDA * N];
 
  959     SType2A = 
new SType2[LDA * N];
 
  964     for(j = 0; j < LDA * N; j++)
 
  969     for(j = 0; j < M2; j++)
 
  974     for(j = 0; j < N2; j++)
 
  981       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
  982       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
  983       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
  984       PrintMatrix(SType1A, M, N, LDA,
"SType1A_before_operation", matlab);
 
  987       PrintMatrix(SType2A, M, N, LDA,
"SType2A_before_operation", matlab);
 
  992     SType1BLAS.
GER(M, N, SType1alpha, SType1x, incx, SType1y, incy, SType1A, LDA);
 
  993     SType2BLAS.
GER(M, N, SType2alpha, SType2x, incx, SType2y, incy, SType2A, LDA);
 
  996       PrintMatrix(SType1A, M, N, LDA, 
"SType1A_after_operation", matlab);
 
  997       PrintMatrix(SType2A, M, N, LDA, 
"SType2A_after_operation", matlab);
 
 1007   GoodTestCount += GoodTestSubcount;
 
 1008   if(verbose || debug) std::cout << 
"GER: " << GoodTestSubcount << 
" of " << GERTESTS << 
" tests were successful." << std::endl;
 
 1009   if(debug) std::cout << std::endl;
 
 1019   GoodTestSubcount = 0;
 
 1029       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
 1034       SType1A = 
new SType1[LDA * P];
 
 1035       SType2A = 
new SType2[LDA * P];
 
 1036       for(j = 0; j < LDA * P; j++)
 
 1042         PrintMatrix(SType1A, M, P, LDA, 
"SType1A", matlab);
 
 1043         PrintMatrix(SType2A, M, P, LDA, 
"SType2A", matlab);
 
 1047       SType1A = 
new SType1[LDA * M];
 
 1048       SType2A = 
new SType2[LDA * M];
 
 1049       for(j = 0; j < LDA * M; j++)
 
 1055         PrintMatrix(SType1A, P, M, LDA, 
"SType1A", matlab);
 
 1056         PrintMatrix(SType2A, P, M, LDA, 
"SType2A", matlab);
 
 1063       SType1B = 
new SType1[LDB * N];
 
 1064       SType2B = 
new SType2[LDB * N];
 
 1065       for(j = 0; j < LDB * N; j++)
 
 1071         PrintMatrix(SType1B, P, N, LDB,
"SType1B", matlab);
 
 1072         PrintMatrix(SType2B, P, N, LDB,
"SType2B", matlab);
 
 1076       SType1B = 
new SType1[LDB * P];
 
 1077       SType2B = 
new SType2[LDB * P];
 
 1078       for(j = 0; j < LDB * P; j++)
 
 1084         PrintMatrix(SType1B, N, P, LDB,
"SType1B", matlab);
 
 1085         PrintMatrix(SType2B, N, P, LDB,
"SType2B", matlab);
 
 1091     SType1C = 
new SType1[LDC * N];
 
 1092     SType2C = 
new SType2[LDC * N];
 
 1093     for(j = 0; j < LDC * N; j++) {
 
 1099       PrintMatrix(SType1C, M, N, LDC, 
"SType1C_before_operation", matlab);
 
 1100       PrintMatrix(SType2C, M, N, LDC, 
"SType2C_before_operation", matlab);
 
 1105     SType2alpha = 
ConvertType(SType1alpha, convertTo);
 
 1109     SType1BLAS.
GEMM(TRANSA, TRANSB, M, N, P, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
 
 1110     SType2BLAS.
GEMM(TRANSA, TRANSB, M, N, P, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
 
 1113       PrintMatrix(SType1C, M, N, LDC, 
"SType1C_after_operation", matlab);
 
 1114       PrintMatrix(SType2C, M, N, LDC, 
"SType2C_after_operation", matlab);
 
 1116     GoodTestSubcount += 
CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
 
 1124   GoodTestCount += GoodTestSubcount;
 
 1125   if(verbose || debug) std::cout << 
"GEMM: " << GoodTestSubcount << 
" of " << GEMMTESTS << 
" tests were successful." << std::endl;
 
 1126   if(debug) std::cout << std::endl;
 
 1134   GoodTestSubcount = 0;
 
 1145       SType1A = 
new SType1[LDA * M];
 
 1146       SType2A = 
new SType2[LDA * M];
 
 1147       for(j = 0; j < LDA * M; j++) {
 
 1153       SType1A = 
new SType1[LDA * N];
 
 1154       SType2A = 
new SType2[LDA * N];
 
 1155       for(j = 0; j < LDA * N; j++) {
 
 1163     SType1B = 
new SType1[LDB * N];
 
 1164     SType2B = 
new SType2[LDB * N];
 
 1165     for(j = 0; j < LDB * N; j++) {
 
 1172     SType1C = 
new SType1[LDC * N];
 
 1173     SType2C = 
new SType2[LDC * N];
 
 1174     for(j = 0; j < LDC * N; j++) {
 
 1181     SType2alpha = 
ConvertType(SType1alpha, convertTo);
 
 1185       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
 1186       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
 1187       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
 1188       std::cout << 
"SType1beta = " << SType1beta << std::endl;
 
 1189       std::cout << 
"SType2beta = " << SType2beta << std::endl;
 
 1191         PrintMatrix(SType1A, M, M, LDA,
"SType1A", matlab);
 
 1192         PrintMatrix(SType2A, M, M, LDA,
"SType2A", matlab);
 
 1194         PrintMatrix(SType1A, N, N, LDA, 
"SType1A", matlab);
 
 1195         PrintMatrix(SType2A, N, N, LDA, 
"SType2A", matlab);
 
 1197       PrintMatrix(SType1B, M, N, LDB,
"SType1B", matlab);
 
 1198       PrintMatrix(SType1C, M, N, LDC,
"SType1C_before_operation", matlab);
 
 1199       PrintMatrix(SType2B, M, N, LDB,
"SType2B", matlab);
 
 1200       PrintMatrix(SType2C, M, N, LDC,
"SType2C_before_operation", matlab);
 
 1204     SType1BLAS.
SYMM(SIDE, UPLO, M, N, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
 
 1205     SType2BLAS.
SYMM(SIDE, UPLO, M, N, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
 
 1208       PrintMatrix(SType1C, M, N, LDC,
"SType1C_after_operation", matlab);
 
 1209       PrintMatrix(SType2C, M, N, LDC,
"SType2C_after_operation", matlab);
 
 1211     GoodTestSubcount += 
CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
 
 1220   GoodTestCount += GoodTestSubcount;
 
 1221   if(verbose || debug) std::cout << 
"SYMM: " << GoodTestSubcount << 
" of " << SYMMTESTS << 
" tests were successful." << std::endl;
 
 1222   if(debug) std::cout << std::endl;
 
 1230   GoodTestSubcount = 0;
 
 1243       SType1A = 
new SType1[LDA * K];
 
 1244       SType2A = 
new SType2[LDA * K];
 
 1245       for(j = 0; j < LDA * K; j++) {
 
 1251       SType1A = 
new SType1[LDA * N];
 
 1252       SType2A = 
new SType2[LDA * N];
 
 1253       for(j = 0; j < LDA * N; j++) {
 
 1261     SType1C = 
new SType1[LDC * N];
 
 1262     SType2C = 
new SType2[LDC * N];
 
 1263     for(j = 0; j < N; j++) {
 
 1268         for(k = 0; k < N; k++)
 
 1273             SType1C[j*LDC+k] = SType1zero;
 
 1275           SType2C[j*LDC+k] = 
ConvertType(SType1C[j*LDC+k], convertTo);
 
 1279         for(k = 0; k < N; k++)
 
 1284             SType1C[j*LDC+k] = SType1zero;
 
 1286           SType2C[j*LDC+k] = 
ConvertType(SType1C[j*LDC+k], convertTo);
 
 1293     SType2alpha = 
ConvertType(SType1alpha, convertTo);
 
 1297       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
 1298       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
 1299       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
 1300       std::cout << 
"SType1beta = " << SType1beta << std::endl;
 
 1301       std::cout << 
"SType2beta = " << SType2beta << std::endl;
 
 1303         PrintMatrix(SType1A, N, K, LDA,
"SType1A", matlab);
 
 1304         PrintMatrix(SType2A, N, K, LDA,
"SType2A", matlab);
 
 1306         PrintMatrix(SType1A, K, N, LDA, 
"SType1A", matlab);
 
 1307         PrintMatrix(SType2A, K, N, LDA, 
"SType2A", matlab);
 
 1309       PrintMatrix(SType1C, N, N, LDC,
"SType1C_before_operation", matlab);
 
 1310       PrintMatrix(SType2C, N, N, LDC,
"SType2C_before_operation", matlab);
 
 1314     SType1BLAS.
SYRK(UPLO, TRANSA, N, K, SType1alpha, SType1A, LDA, SType1beta, SType1C, LDC);
 
 1315     SType2BLAS.
SYRK(UPLO, TRANSA, N, K, SType2alpha, SType2A, LDA, SType2beta, SType2C, LDC);
 
 1318       PrintMatrix(SType1C, N, N, LDC,
"SType1C_after_operation", matlab);
 
 1319       PrintMatrix(SType2C, N, N, LDC,
"SType2C_after_operation", matlab);
 
 1321     GoodTestSubcount += 
CompareMatrices(SType1C, SType2C, N, N, LDC, TOL);
 
 1328   GoodTestCount += GoodTestSubcount;
 
 1329   if(verbose || debug) std::cout << 
"SYRK: " << GoodTestSubcount << 
" of " << SYRKTESTS << 
" tests were successful." << std::endl;
 
 1330   if(debug) std::cout << std::endl;
 
 1338   GoodTestSubcount = 0;
 
 1349     SType1B = 
new SType1[LDB * N];
 
 1350     SType2B = 
new SType2[LDB * N];
 
 1364       SType1A = 
new SType1[LDA * M];
 
 1365       SType2A = 
new SType2[LDA * M];
 
 1367       for(j = 0; j < M; j++)
 
 1372           for(k = 0; k < M; k++)
 
 1377               SType1A[j*LDA+k] = SType1zero;
 
 1379             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1383                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1386                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1388                 SType1A[j*LDA+k] = SType1one;
 
 1389                 SType2A[j*LDA+k] = SType2one;
 
 1396           for(k = 0; k < M; k++)
 
 1401               SType1A[j*LDA+k] = SType1zero;
 
 1403             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1407                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1410                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1412                 SType1A[j*LDA+k] = SType1one;
 
 1413                 SType2A[j*LDA+k] = SType2one;
 
 1427       SType1A = 
new SType1[LDA * N];
 
 1428       SType2A = 
new SType2[LDA * N];
 
 1430       for(j = 0; j < N; j++)
 
 1435           for(k = 0; k < N; k++)
 
 1440               SType1A[j*LDA+k] = SType1zero;
 
 1442             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1446                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1449                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1451                 SType1A[j*LDA+k] = SType1one;
 
 1452                 SType2A[j*LDA+k] = SType2one;
 
 1459           for(k = 0; k < N; k++)
 
 1464               SType1A[j*LDA+k] = SType1zero;
 
 1466             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1470                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1473                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1475                 SType1A[j*LDA+k] = SType1one;
 
 1476                 SType2A[j*LDA+k] = SType2one;
 
 1485     for(j = 0; j < N; j++) {
 
 1486       for(k = 0; k < M; k++) {
 
 1488         SType2B[j*LDB+k] = 
ConvertType(SType1B[j*LDB+k], convertTo);
 
 1492     SType2alpha = 
ConvertType(SType1alpha, convertTo);
 
 1495       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
 1496       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
 1497       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
 1499         PrintMatrix(SType1A, M, M, LDA, 
"SType1A", matlab);
 
 1500         PrintMatrix(SType2A, M, M, LDA, 
"SType2A", matlab);
 
 1502         PrintMatrix(SType1A, N, N, LDA, 
"SType1A", matlab);
 
 1503         PrintMatrix(SType2A, N, N, LDA, 
"SType2A", matlab);
 
 1505       PrintMatrix(SType1B, M, N, LDB,
"SType1B_before_operation", matlab);
 
 1506       PrintMatrix(SType2B, M, N, LDB,
"SType2B_before_operation", matlab);
 
 1509     SType1BLAS.
TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
 
 1510     SType2BLAS.
TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
 
 1513       PrintMatrix(SType1B, M, N, LDB, 
"SType1B_after_operation", matlab);
 
 1514       PrintMatrix(SType2B, M, N, LDB, 
"SType2B_after_operation", matlab);
 
 1516     GoodTestSubcount += 
CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
 
 1522   GoodTestCount += GoodTestSubcount;
 
 1523   if(verbose || debug) std::cout << 
"TRMM: " << GoodTestSubcount << 
" of " << TRMMTESTS << 
" tests were successful." << std::endl;
 
 1524   if(debug) std::cout << std::endl;
 
 1532   GoodTestSubcount = 0;
 
 1543     SType1B = 
new SType1[LDB * N];
 
 1544     SType2B = 
new SType2[LDB * N];
 
 1561       SType1A = 
new SType1[LDA * M];
 
 1562       SType2A = 
new SType2[LDA * M];
 
 1564       for(j = 0; j < M; j++)
 
 1569           for(k = 0; k < M; k++)
 
 1574               SType1A[j*LDA+k] = SType1zero;
 
 1576             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1580                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1583                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1585                 SType1A[j*LDA+k] = SType1one;
 
 1586                 SType2A[j*LDA+k] = SType2one;
 
 1593           for(k = 0; k < M; k++)
 
 1598               SType1A[j*LDA+k] = SType1zero;
 
 1600             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1604                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1607                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1609                 SType1A[j*LDA+k] = SType1one;
 
 1610                 SType2A[j*LDA+k] = SType2one;
 
 1624       SType1A = 
new SType1[LDA * N];
 
 1625       SType2A = 
new SType2[LDA * N];
 
 1627       for(j = 0; j < N; j++)
 
 1632           for(k = 0; k < N; k++)
 
 1637               SType1A[j*LDA+k] = SType1zero;
 
 1639             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1643                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1646                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1648                 SType1A[j*LDA+k] = SType1one;
 
 1649                 SType2A[j*LDA+k] = SType2one;
 
 1656           for(k = 0; k < N; k++)
 
 1661               SType1A[j*LDA+k] = SType1zero;
 
 1663             SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1667                 while (SType1A[j*LDA+k] == SType1zero) {
 
 1670                 SType2A[j*LDA+k] = 
ConvertType(SType1A[j*LDA+k], convertTo);
 
 1672                 SType1A[j*LDA+k] = SType1one;
 
 1673                 SType2A[j*LDA+k] = SType2one;
 
 1682     for(j = 0; j < N; j++)
 
 1684       for(k = 0; k < M; k++)
 
 1687         SType2B[j*LDB+k] = 
ConvertType(SType1B[j*LDB+k], convertTo);
 
 1692     SType2alpha = 
ConvertType(SType1alpha, convertTo);
 
 1696       std::cout << 
"Test #" << TotalTestCount << 
" --" << std::endl;
 
 1701       std::cout << 
"M="<<M << 
"\t" << 
"N="<<N << 
"\t" << 
"LDA="<<LDA << 
"\t" << 
"LDB="<<LDB << std::endl;
 
 1702       std::cout << 
"SType1alpha = " << SType1alpha << std::endl;
 
 1703       std::cout << 
"SType2alpha = " << SType2alpha << std::endl;
 
 1705         PrintMatrix(SType1A, M, M, LDA, 
"SType1A", matlab);
 
 1706         PrintMatrix(SType2A, M, M, LDA, 
"SType2A", matlab);
 
 1708         PrintMatrix(SType1A, N, N, LDA, 
"SType1A", matlab);
 
 1709         PrintMatrix(SType2A, N, N, LDA, 
"SType2A", matlab);
 
 1711       PrintMatrix(SType1B, M, N, LDB, 
"SType1B_before_operation", matlab);
 
 1712       PrintMatrix(SType2B, M, N, LDB, 
"SType2B_before_operation", matlab);
 
 1716     SType1BLAS.
TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
 
 1717     SType2BLAS.
TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
 
 1721       PrintMatrix(SType1B, M, N, LDB, 
"SType1B_after_operation", matlab);
 
 1722       PrintMatrix(SType2B, M, N, LDB, 
"SType2B_after_operation", matlab);
 
 1726       std::cout << 
"FAILED TEST!!!!!!" << std::endl;
 
 1727     GoodTestSubcount += 
CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
 
 1734   GoodTestCount += GoodTestSubcount;
 
 1735   if(verbose || debug) std::cout << 
"TRSM: " << GoodTestSubcount << 
" of " << TRSMTESTS << 
" tests were successful." << std::endl;
 
 1736   if(debug) std::cout << std::endl;
 
 1741 #ifdef HAVE_TEUCHOS_ARPREC 
 1746   if((((TotalTestCount - 1) - GoodTestCount) != 0) || (verbose) || (debug))
 
 1748     std::cout << GoodTestCount << 
" of " << (TotalTestCount - 1) << 
" total tests were successful." << std::endl;
 
 1751   if ((TotalTestCount-1) == GoodTestCount) {
 
 1752     std::cout << 
"End Result: TEST PASSED" << std::endl;
 
 1756   std::cout << 
"End Result: TEST FAILED" << std::endl;
 
 1757   return (TotalTestCount-GoodTestCount-1);
 
 1760 template<
typename TYPE>
 
 1778 template<
typename TYPE>
 
 1781   std::cout << Name << 
" =" << std::endl;
 
 1783   if(Matlab) std::cout << 
"[";
 
 1784   for(i = 0; i < Size; i++)
 
 1786     std::cout << Vector[i] << 
" ";
 
 1788   if(Matlab) std::cout << 
"]";
 
 1791     std::cout << std::endl << std::endl;
 
 1795     std::cout << 
";" << std::endl;
 
 1799   template<
typename TYPE>
 
 1804     std::cout << Name << 
" =" << std::endl;
 
 1806     for(i = 0; i < Rows; i++)
 
 1808       for(j = 0; j < Columns; j++)
 
 1810         std::cout << Matrix[i + (j * LDM)] << 
" ";
 
 1812       std::cout << std::endl;
 
 1814     std::cout << std::endl;
 
 1818     std::cout << Name << 
" = ";
 
 1821     for(i = 0; i < Rows; i++)
 
 1824       for(j = 0; j < Columns; j++)
 
 1826         std::cout << Matrix[i + (j * LDM)] << 
" ";
 
 1830     std::cout << 
"];" << std::endl;
 
 1834   template<
typename TYPE1, 
typename TYPE2>
 
 1843   return( temp2 < Tolerance );
 
 1850   template<
typename TYPE1, 
typename TYPE2>
 
 1859   for(i = 0; i < Size; i++)
 
 1862     sum2 += temp2*temp2;
 
 1871   if (temp2 > Tolerance)
 
 1880   template<
typename TYPE1, 
typename TYPE2>
 
 1889   for(j = 0; j < Columns; j++)
 
 1891     for(i = 0; i < Rows; i++)
 
 1904   if (temp2 > Tolerance)
 
 1911 template<
typename TYPE1, 
typename TYPE2>
 
 1914   return static_cast<TYPE2
>(T1);
 
 1917 #ifdef HAVE_TEUCHOS_ARPREC 
 1925 #ifdef HAVE_TEUCHOS_QD 
 1929   return to_double(T1);
 
 1933 #ifdef HAVE_TEUCHOS_GNU_MP 
 1975   if(r == 1 || r == 2)
 
void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const 
Solves the matrix equations: op(A)*X=alpha*B or X*op(A)=alpha*B where X and B are m by n matrices...
 
void GER(const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy, ScalarType *A, const OrdinalType &lda) const 
Performs the rank 1 operation: A <- alpha*x*y'+A. 
 
void AXPY(const OrdinalType &n, const alpha_type alpha, const x_type *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const 
Perform the operation: y <- y+alpha*x. 
 
void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType &n, const A_type *A, const OrdinalType &lda, ScalarType *x, const OrdinalType &incx) const 
Performs the matrix-vector operation: x <- A*x or x <- A'*x where A is a unit/non-unit n by n upper/low...
 
Templated interface class to BLAS routines. 
 
void GEMV(ETransp trans, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const x_type *x, const OrdinalType &incx, const beta_type beta, ScalarType *y, const OrdinalType &incy) const 
Performs the matrix-vector operation: y <- alpha*A*x+beta*y or y <- alpha*A'*x+beta*y where A is a gene...
 
void ROT(const OrdinalType &n, ScalarType *dx, const OrdinalType &incx, ScalarType *dy, const OrdinalType &incy, MagnitudeType *c, ScalarType *s) const 
Applies a Givens plane rotation. 
 
Basic wall-clock timer class. 
 
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x. 
 
ScalarTraits< ScalarType >::magnitudeType NRM2(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const 
Compute the 2-norm of the vector x. 
 
OrdinalType IAMAX(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const 
Return the index of the element of x with the maximum magnitude. 
 
void GEMM(ETransp transa, ETransp transb, const OrdinalType &m, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const 
General matrix-matrix multiply. 
 
void PrintMatrix(TYPE *Matrix, OType Rows, OType Columns, OType LDM, std::string Name, bool Matlab=0)
 
Initialize, finalize, and query the global MPI session. 
 
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EDiagChar[]
 
This structure defines some basic traits for a scalar field type. 
 
ScalarTraits< ScalarType >::magnitudeType ASUM(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx) const 
Sum the absolute values of the entries of x. 
 
Teuchos::EUplo RandomUPLO()
 
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const 
Copy the vector x to the vector y. 
 
TYPE GetRandom(TYPE, TYPE)
 
Teuchos::ESide RandomSIDE()
 
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char EUploChar[]
 
Teuchos::EDiag RandomDIAG()
 
ScalarType DOT(const OrdinalType &n, const x_type *x, const OrdinalType &incx, const y_type *y, const OrdinalType &incy) const 
Form the dot product of the vectors x and y. 
 
bool CompareVectors(TYPE1 *Vector1, TYPE2 *Vector2, OType Size, TYPE2 Tolerance)
 
bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, TYPE2 Tolerance)
 
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ESideChar[]
 
std::string Teuchos_Version()
 
void PrintVector(TYPE *Vector, OType Size, std::string Name, bool Matlab=0)
 
int main(int argc, char *argv[])
 
bool CompareMatrices(TYPE1 *Matrix1, TYPE2 *Matrix2, OType Rows, OType Columns, OType LDM, TYPE2 Tolerance)
 
TYPE2 ConvertType(TYPE1, TYPE2)
 
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a. 
 
void SYMM(ESide side, EUplo uplo, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const B_type *B, const OrdinalType &ldb, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const 
Performs the matrix-matrix operation: C <- alpha*A*B+beta*C or C <- alpha*B*A+beta*C where A is an m ...
 
void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType &m, const OrdinalType &n, const alpha_type alpha, const A_type *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb) const 
Performs the matrix-matrix operation: B <- alpha*op(A)*B or B <- alpha*B*op(A) where op(A) is an unit...
 
A MPI utilities class, providing methods for initializing, finalizing, and querying the global MPI se...
 
void SCAL(const OrdinalType &n, const ScalarType &alpha, ScalarType *x, const OrdinalType &incx) const 
Scale the vector x by the constant alpha. 
 
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const 
Computes a Givens plane rotation. 
 
TEUCHOSNUMERICS_LIB_DLL_EXPORT const char ETranspChar[]
 
void SYRK(EUplo uplo, ETransp trans, const OrdinalType &n, const OrdinalType &k, const alpha_type alpha, const A_type *A, const OrdinalType &lda, const beta_type beta, ScalarType *C, const OrdinalType &ldc) const 
Performs the symmetric rank k operation: C <- alpha*A*A'+beta*C or C <- alpha*A'*A+beta*C, where alpha and beta are scalars, C is an n by n symmetric matrix and A is an n by k matrix in the first case or k by n matrix in the second case. 
 
Teuchos::ETransp RandomTRANS()