33 #ifdef HAVE_TEUCHOS_ARPREC
34 #define SType1 mp_real
35 #elif defined(HAVE_TEUCHOS_QD)
36 #define SType1 dd_real
37 #elif defined(HAVE_TEUCHOS_GNU_MP)
38 #define SType1 mpf_class
76 template<
typename TYPE>
87 template<
typename TYPE>
88 void PrintVector(TYPE* Vector,
OType Size, std::string Name,
bool Matlab = 0);
90 template<
typename TYPE>
93 template<
typename TYPE1,
typename TYPE2>
94 bool CompareScalars(TYPE1 Scalar1, TYPE2 Scalar2, TYPE2 Tolerance );
96 template<
typename TYPE1,
typename TYPE2>
99 template<
typename TYPE1,
typename TYPE2>
104 template<
typename TYPE1,
typename TYPE2>
107 #ifdef HAVE_TEUCHOS_ARPREC
112 #ifdef HAVE_TEUCHOS_QD
117 #ifdef HAVE_TEUCHOS_GNU_MP
128 int main(
int argc,
char *argv[])
130 #ifdef HAVE_TEUCHOS_ARPREC
139 bool InvalidCmdLineArgs = 0;
141 for(i = 1; i < argc; i++)
143 if(argv[i][0] ==
'-')
154 InvalidCmdLineArgs = 1;
164 InvalidCmdLineArgs = 1;
174 InvalidCmdLineArgs = 1;
178 InvalidCmdLineArgs = 1;
187 if(InvalidCmdLineArgs || (argc > 4))
189 std::cout <<
"Invalid command line arguments detected. Use the following flags:" << std::endl
190 <<
"\t -v enables verbose mode (reports number of failed/successful tests)" << std::endl
191 <<
"\t -d enables debug mode (same as verbose with output of each test, not recommended for large numbers of tests)" << std::endl
192 <<
"\t -m enables matlab-style output; only has an effect if debug mode is enabled" << std::endl;
205 SType1 SType1alpha, SType1beta;
211 SType2 SType2alpha, SType2beta;
212 SType1 SType1ASUMresult, SType1DOTresult, SType1NRM2result, SType1SINresult, SType1COSresult;
213 SType2 SType2ASUMresult, SType2DOTresult, SType2NRM2result, SType2SINresult, SType2COSresult;
215 OType SType1IAMAXresult;
216 OType SType2IAMAXresult;
217 OType TotalTestCount = 1, GoodTestSubcount, GoodTestCount = 0, M, M2, N, N2, K, P, LDA, LDB, LDC, Mx, My;
225 std::srand(time(NULL));
232 GoodTestSubcount = 0;
240 SType2COSresult =
ConvertType(SType1COSresult, convertTo);
242 SType2SINresult =
ConvertType(SType1SINresult, convertTo);
246 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
247 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
248 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
249 std::cout <<
"SType1beta = " << SType1beta << std::endl;
250 std::cout <<
"SType2beta = " << SType2beta << std::endl;
253 SType1BLAS.
ROTG(&SType1alpha, &SType1beta, &SType1COSresult, &SType1SINresult);
254 SType2BLAS.
ROTG(&SType2alpha, &SType2beta, &SType2COSresult, &SType2SINresult);
257 std::cout <<
"SType1 ROTG COS result: " << SType1COSresult << std::endl;
258 std::cout <<
"SType2 ROTG COS result: " << SType2COSresult << std::endl;
259 std::cout <<
"SType1 ROTG SIN result: " << SType1SINresult << std::endl;
260 std::cout <<
"SType2 ROTG SIN result: " << SType2SINresult << std::endl;
262 GoodTestSubcount += (
CompareScalars(SType1COSresult, SType2COSresult, TOL) &&
265 GoodTestCount += GoodTestSubcount;
266 if(verbose || debug) std::cout <<
"ROTG: " << GoodTestSubcount <<
" of " << ROTGTESTS <<
" tests were successful." << std::endl;
267 if(debug) std::cout << std::endl;
277 GoodTestSubcount = 0;
282 if (incx == 0) incx = 1;
283 if (incy == 0) incy = 1;
285 Mx = M*std::abs(incx);
286 My = M*std::abs(incy);
287 if (Mx == 0) { Mx = 1; }
288 if (My == 0) { My = 1; }
293 for(j = 0; j < Mx; j++)
298 for(j = 0; j < My; j++)
309 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
310 std::cout <<
"c1 = " << c1 <<
", s1 = " << s1 << std::endl;
311 std::cout <<
"c2 = " << c2 <<
", s2 = " << s2 << std::endl;
312 std::cout <<
"incx = " << incx <<
", incy = " << incy << std::endl;
314 PrintVector(SType1y, My,
"SType1y_before_operation", matlab);
316 PrintVector(SType2y, My,
"SType2y_before_operation", matlab);
319 SType1BLAS.
ROT(M, SType1x, incx, SType1y, incy, &c1, &s1);
320 SType2BLAS.
ROT(M, SType2x, incx, SType2y, incy, &c2, &s2);
323 PrintVector(SType1y, My,
"SType1y_after_operation", matlab);
324 PrintVector(SType2y, My,
"SType2y_after_operation", matlab);
326 GoodTestSubcount += (
CompareVectors(SType1x, SType2x, Mx, TOL) &&
333 GoodTestCount += GoodTestSubcount;
334 if(verbose || debug) std::cout <<
"ROT: " << GoodTestSubcount <<
" of " << ROTTESTS <<
" tests were successful." << std::endl;
335 if(debug) std::cout << std::endl;
343 GoodTestSubcount = 0;
352 for(j = 0; j < M2; j++)
359 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
364 SType1ASUMresult = SType1BLAS.
ASUM(M, SType1x, incx);
365 SType2ASUMresult = SType2BLAS.
ASUM(M, SType2x, incx);
368 std::cout <<
"SType1 ASUM result: " << SType1ASUMresult << std::endl;
369 std::cout <<
"SType2 ASUM result: " << SType2ASUMresult << std::endl;
371 GoodTestSubcount +=
CompareScalars(SType1ASUMresult, SType2ASUMresult, TOL);
375 GoodTestCount += GoodTestSubcount;
376 if(verbose || debug) std::cout <<
"ASUM: " << GoodTestSubcount <<
" of " << ASUMTESTS <<
" tests were successful." << std::endl;
377 if(debug) std::cout << std::endl;
386 GoodTestSubcount = 0;
392 Mx = M*std::abs(incx);
393 My = M*std::abs(incy);
394 if (Mx == 0) { Mx = 1; }
395 if (My == 0) { My = 1; }
400 for(j = 0; j < Mx; j++)
405 for(j = 0; j < My; j++)
414 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
415 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
416 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
418 PrintVector(SType1y, My,
"SType1y_before_operation", matlab);
420 PrintVector(SType2y, My,
"SType2y_before_operation", matlab);
423 SType1BLAS.
AXPY(M, SType1alpha, SType1x, incx, SType1y, incy);
424 SType2BLAS.
AXPY(M, SType2alpha, SType2x, incx, SType2y, incy);
427 PrintVector(SType1y, My,
"SType1y_after_operation", matlab);
428 PrintVector(SType2y, My,
"SType2y_after_operation", matlab);
436 GoodTestCount += GoodTestSubcount;
437 if(verbose || debug) std::cout <<
"AXPY: " << GoodTestSubcount <<
" of " << AXPYTESTS <<
" tests were successful." << std::endl;
438 if(debug) std::cout << std::endl;
446 GoodTestSubcount = 0;
452 Mx = M*std::abs(incx);
453 My = M*std::abs(incy);
454 if (Mx == 0) { Mx = 1; }
455 if (My == 0) { My = 1; }
460 for(j = 0; j < Mx; j++)
465 for(j = 0; j < My; j++)
472 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
474 PrintVector(SType1y, My,
"SType1y_before_operation", matlab);
476 PrintVector(SType2y, My,
"SType2y_before_operation", matlab);
479 SType1BLAS.
COPY(M, SType1x, incx, SType1y, incy);
480 SType2BLAS.
COPY(M, SType2x, incx, SType2y, incy);
483 PrintVector(SType1y, My,
"SType1y_after_operation", matlab);
484 PrintVector(SType2y, My,
"SType2y_after_operation", matlab);
492 GoodTestCount += GoodTestSubcount;
if(verbose || debug) std::cout <<
"COPY: " << GoodTestSubcount <<
" of " << COPYTESTS <<
" tests were successful." << std::endl;
493 if(debug) std::cout << std::endl;
501 GoodTestSubcount = 0;
507 Mx = M*std::abs(incx);
508 My = M*std::abs(incy);
509 if (Mx == 0) { Mx = 1; }
510 if (My == 0) { My = 1; }
515 for(j = 0; j < Mx; j++)
520 for(j = 0; j < My; j++)
527 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
534 SType1DOTresult = SType1BLAS.
DOT(M, SType1x, incx, SType1y, incy);
535 SType2DOTresult = SType2BLAS.
DOT(M, SType2x, incx, SType2y, incy);
538 std::cout <<
"SType1 DOT result: " << SType1DOTresult << std::endl;
539 std::cout <<
"SType2 DOT result: " << SType2DOTresult << std::endl;
541 GoodTestSubcount +=
CompareScalars(SType1DOTresult, SType2DOTresult, TOL);
547 GoodTestCount += GoodTestSubcount;
548 if(verbose || debug) std::cout <<
"DOT: " << GoodTestSubcount <<
" of " << DOTTESTS <<
" tests were successful." << std::endl;
549 if(debug) std::cout << std::endl;
557 GoodTestSubcount = 0;
565 for(j = 0; j < M2; j++)
572 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
577 SType1NRM2result = SType1BLAS.
NRM2(M, SType1x, incx);
578 SType2NRM2result = SType2BLAS.
NRM2(M, SType2x, incx);
581 std::cout <<
"SType1 NRM2 result: " << SType1NRM2result << std::endl;
582 std::cout <<
"SType2 NRM2 result: " << SType2NRM2result << std::endl;
584 GoodTestSubcount +=
CompareScalars(SType1NRM2result, SType2NRM2result, TOL);
588 GoodTestCount += GoodTestSubcount;
if(verbose || debug) std::cout <<
"NRM2: " << GoodTestSubcount <<
" of " << NRM2TESTS <<
" tests were successful." << std::endl;
589 if(debug) std::cout << std::endl;
597 GoodTestSubcount = 0;
608 for(j = 0; j < M2; j++)
617 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
618 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
619 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
620 PrintVector(SType1x, M2,
"SType1x_before_operation", matlab);
621 PrintVector(SType2x, M2,
"SType2x_before_operation", matlab);
624 SType1BLAS.
SCAL(M, SType1alpha, SType1x, incx);
625 SType2BLAS.
SCAL(M, SType2alpha, SType2x, incx);
628 PrintVector(SType1x, M2,
"SType1x_after_operation", matlab);
629 PrintVector(SType2x, M2,
"SType2x_after_operation", matlab);
635 GoodTestCount += GoodTestSubcount;
636 if(verbose || debug) std::cout <<
"SCAL: " << GoodTestSubcount <<
" of " << SCALTESTS <<
" tests were successful." << std::endl;
637 if(debug) std::cout << std::endl;
645 GoodTestSubcount = 0;
653 for(j = 0; j < M2; j++)
660 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
665 SType1IAMAXresult = SType1BLAS.
IAMAX(M, SType1x, incx);
666 SType2IAMAXresult = SType2BLAS.
IAMAX(M, SType2x, incx);
669 std::cout <<
"SType1 IAMAX result: " << SType1IAMAXresult << std::endl;
670 std::cout <<
"SType2 IAMAX result: " << SType2IAMAXresult << std::endl;
672 GoodTestSubcount += (SType1IAMAXresult == SType2IAMAXresult);
676 GoodTestCount += GoodTestSubcount;
677 if(verbose || debug) std::cout <<
"IAMAX: " << GoodTestSubcount <<
" of " << IAMAXTESTS <<
" tests were successful." << std::endl;
678 if(debug) std::cout << std::endl;
688 GoodTestSubcount = 0;
706 M2 = M*std::abs(incy);
707 N2 = N*std::abs(incx);
709 M2 = N*std::abs(incy);
710 N2 = M*std::abs(incx);
723 SType1A =
new SType1[LDA * N];
726 SType2A =
new SType2[LDA * N];
730 for(j = 0; j < LDA * N; j++)
735 for(j = 0; j < N2; j++)
740 for(j = 0; j < M2; j++)
747 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
748 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
749 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
750 std::cout <<
"SType1beta = " << SType1beta << std::endl;
751 std::cout <<
"SType2beta = " << SType2beta << std::endl;
752 PrintMatrix(SType1A, M, N, LDA,
"SType1A", matlab);
754 PrintVector(SType1y, M2,
"SType1y_before_operation", matlab);
755 PrintMatrix(SType2A, M, N, LDA,
"SType2A", matlab);
757 PrintVector(SType2y, M2,
"SType2y_before_operation", matlab);
760 SType1BLAS.
GEMV(TRANS, M, N, SType1alpha, SType1A, LDA, SType1x, incx, SType1beta, SType1y, incy);
761 SType2BLAS.
GEMV(TRANS, M, N, SType2alpha, SType2A, LDA, SType2x, incx, SType2beta, SType2y, incy);
764 PrintVector(SType1y, M2,
"SType1y_after_operation", matlab);
765 PrintVector(SType2y, M2,
"SType2y_after_operation", matlab);
775 GoodTestCount += GoodTestSubcount;
776 if(verbose || debug) std::cout <<
"GEMV: " << GoodTestSubcount <<
" of " << GEMVTESTS <<
" tests were successful." << std::endl;
777 if(debug) std::cout << std::endl;
785 GoodTestSubcount = 0;
800 N2 = N*std::abs(incx);
804 for(j = 0; j < N2; j++)
814 SType1A =
new SType1[LDA * N];
815 SType2A =
new SType2[LDA * N];
817 for(j = 0; j < N; j++)
822 for(k = 0; k < N; k++)
827 SType1A[j*LDA+k] = SType1zero;
829 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
833 while (SType1A[j*LDA+k] == SType1zero) {
836 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
838 SType1A[j*LDA+k] = SType1one;
839 SType2A[j*LDA+k] = SType2one;
846 for(k = 0; k < N; k++)
851 SType1A[j*LDA+k] = SType1zero;
853 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
857 while (SType1A[j*LDA+k] == SType1zero) {
860 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
862 SType1A[j*LDA+k] = SType1one;
863 SType2A[j*LDA+k] = SType2one;
872 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
874 PrintVector(SType1x, N2,
"SType1x_before_operation", matlab);
875 PrintMatrix(SType2A, N, N, LDA,
"SType2A", matlab);
876 PrintVector(SType2x, N2,
"SType2x_before_operation", matlab);
879 SType1BLAS.
TRMV(UPLO, TRANSA, DIAG, N, SType1A, LDA, SType1x, incx);
880 SType2BLAS.
TRMV(UPLO, TRANSA, DIAG, N, SType2A, LDA, SType2x, incx);
883 PrintVector(SType1x, N2,
"SType1x_after_operation", matlab);
884 PrintVector(SType2x, N2,
"SType2x_after_operation", matlab);
892 GoodTestCount += GoodTestSubcount;
893 if(verbose || debug) std::cout <<
"TRMV: " << GoodTestSubcount <<
" of " << TRMVTESTS <<
" tests were successful." << std::endl;
894 if(debug) std::cout << std::endl;
902 GoodTestSubcount = 0;
916 M2 = M*std::abs(incx);
917 N2 = N*std::abs(incy);
924 SType1A =
new SType1[LDA * N];
927 SType2A =
new SType2[LDA * N];
932 for(j = 0; j < LDA * N; j++)
937 for(j = 0; j < M2; j++)
942 for(j = 0; j < N2; j++)
949 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
950 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
951 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
952 PrintMatrix(SType1A, M, N, LDA,
"SType1A_before_operation", matlab);
955 PrintMatrix(SType2A, M, N, LDA,
"SType2A_before_operation", matlab);
960 SType1BLAS.
GER(M, N, SType1alpha, SType1x, incx, SType1y, incy, SType1A, LDA);
961 SType2BLAS.
GER(M, N, SType2alpha, SType2x, incx, SType2y, incy, SType2A, LDA);
964 PrintMatrix(SType1A, M, N, LDA,
"SType1A_after_operation", matlab);
965 PrintMatrix(SType2A, M, N, LDA,
"SType2A_after_operation", matlab);
975 GoodTestCount += GoodTestSubcount;
976 if(verbose || debug) std::cout <<
"GER: " << GoodTestSubcount <<
" of " << GERTESTS <<
" tests were successful." << std::endl;
977 if(debug) std::cout << std::endl;
987 GoodTestSubcount = 0;
997 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
1002 SType1A =
new SType1[LDA * P];
1003 SType2A =
new SType2[LDA * P];
1004 for(j = 0; j < LDA * P; j++)
1010 PrintMatrix(SType1A, M, P, LDA,
"SType1A", matlab);
1011 PrintMatrix(SType2A, M, P, LDA,
"SType2A", matlab);
1015 SType1A =
new SType1[LDA * M];
1016 SType2A =
new SType2[LDA * M];
1017 for(j = 0; j < LDA * M; j++)
1023 PrintMatrix(SType1A, P, M, LDA,
"SType1A", matlab);
1024 PrintMatrix(SType2A, P, M, LDA,
"SType2A", matlab);
1031 SType1B =
new SType1[LDB * N];
1032 SType2B =
new SType2[LDB * N];
1033 for(j = 0; j < LDB * N; j++)
1039 PrintMatrix(SType1B, P, N, LDB,
"SType1B", matlab);
1040 PrintMatrix(SType2B, P, N, LDB,
"SType2B", matlab);
1044 SType1B =
new SType1[LDB * P];
1045 SType2B =
new SType2[LDB * P];
1046 for(j = 0; j < LDB * P; j++)
1052 PrintMatrix(SType1B, N, P, LDB,
"SType1B", matlab);
1053 PrintMatrix(SType2B, N, P, LDB,
"SType2B", matlab);
1059 SType1C =
new SType1[LDC * N];
1060 SType2C =
new SType2[LDC * N];
1061 for(j = 0; j < LDC * N; j++) {
1067 PrintMatrix(SType1C, M, N, LDC,
"SType1C_before_operation", matlab);
1068 PrintMatrix(SType2C, M, N, LDC,
"SType2C_before_operation", matlab);
1073 SType2alpha =
ConvertType(SType1alpha, convertTo);
1077 SType1BLAS.
GEMM(TRANSA, TRANSB, M, N, P, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
1078 SType2BLAS.
GEMM(TRANSA, TRANSB, M, N, P, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
1081 PrintMatrix(SType1C, M, N, LDC,
"SType1C_after_operation", matlab);
1082 PrintMatrix(SType2C, M, N, LDC,
"SType2C_after_operation", matlab);
1084 GoodTestSubcount +=
CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
1092 GoodTestCount += GoodTestSubcount;
1093 if(verbose || debug) std::cout <<
"GEMM: " << GoodTestSubcount <<
" of " << GEMMTESTS <<
" tests were successful." << std::endl;
1094 if(debug) std::cout << std::endl;
1102 GoodTestSubcount = 0;
1113 SType1A =
new SType1[LDA * M];
1114 SType2A =
new SType2[LDA * M];
1115 for(j = 0; j < LDA * M; j++) {
1121 SType1A =
new SType1[LDA * N];
1122 SType2A =
new SType2[LDA * N];
1123 for(j = 0; j < LDA * N; j++) {
1131 SType1B =
new SType1[LDB * N];
1132 SType2B =
new SType2[LDB * N];
1133 for(j = 0; j < LDB * N; j++) {
1140 SType1C =
new SType1[LDC * N];
1141 SType2C =
new SType2[LDC * N];
1142 for(j = 0; j < LDC * N; j++) {
1149 SType2alpha =
ConvertType(SType1alpha, convertTo);
1153 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
1154 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
1155 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
1156 std::cout <<
"SType1beta = " << SType1beta << std::endl;
1157 std::cout <<
"SType2beta = " << SType2beta << std::endl;
1159 PrintMatrix(SType1A, M, M, LDA,
"SType1A", matlab);
1160 PrintMatrix(SType2A, M, M, LDA,
"SType2A", matlab);
1162 PrintMatrix(SType1A, N, N, LDA,
"SType1A", matlab);
1163 PrintMatrix(SType2A, N, N, LDA,
"SType2A", matlab);
1165 PrintMatrix(SType1B, M, N, LDB,
"SType1B", matlab);
1166 PrintMatrix(SType1C, M, N, LDC,
"SType1C_before_operation", matlab);
1167 PrintMatrix(SType2B, M, N, LDB,
"SType2B", matlab);
1168 PrintMatrix(SType2C, M, N, LDC,
"SType2C_before_operation", matlab);
1172 SType1BLAS.
SYMM(SIDE, UPLO, M, N, SType1alpha, SType1A, LDA, SType1B, LDB, SType1beta, SType1C, LDC);
1173 SType2BLAS.
SYMM(SIDE, UPLO, M, N, SType2alpha, SType2A, LDA, SType2B, LDB, SType2beta, SType2C, LDC);
1176 PrintMatrix(SType1C, M, N, LDC,
"SType1C_after_operation", matlab);
1177 PrintMatrix(SType2C, M, N, LDC,
"SType2C_after_operation", matlab);
1179 GoodTestSubcount +=
CompareMatrices(SType1C, SType2C, M, N, LDC, TOL);
1188 GoodTestCount += GoodTestSubcount;
1189 if(verbose || debug) std::cout <<
"SYMM: " << GoodTestSubcount <<
" of " << SYMMTESTS <<
" tests were successful." << std::endl;
1190 if(debug) std::cout << std::endl;
1198 GoodTestSubcount = 0;
1211 SType1A =
new SType1[LDA * K];
1212 SType2A =
new SType2[LDA * K];
1213 for(j = 0; j < LDA * K; j++) {
1219 SType1A =
new SType1[LDA * N];
1220 SType2A =
new SType2[LDA * N];
1221 for(j = 0; j < LDA * N; j++) {
1229 SType1C =
new SType1[LDC * N];
1230 SType2C =
new SType2[LDC * N];
1231 for(j = 0; j < N; j++) {
1236 for(k = 0; k < N; k++)
1241 SType1C[j*LDC+k] = SType1zero;
1243 SType2C[j*LDC+k] =
ConvertType(SType1C[j*LDC+k], convertTo);
1247 for(k = 0; k < N; k++)
1252 SType1C[j*LDC+k] = SType1zero;
1254 SType2C[j*LDC+k] =
ConvertType(SType1C[j*LDC+k], convertTo);
1261 SType2alpha =
ConvertType(SType1alpha, convertTo);
1265 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
1266 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
1267 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
1268 std::cout <<
"SType1beta = " << SType1beta << std::endl;
1269 std::cout <<
"SType2beta = " << SType2beta << std::endl;
1271 PrintMatrix(SType1A, N, K, LDA,
"SType1A", matlab);
1272 PrintMatrix(SType2A, N, K, LDA,
"SType2A", matlab);
1274 PrintMatrix(SType1A, K, N, LDA,
"SType1A", matlab);
1275 PrintMatrix(SType2A, K, N, LDA,
"SType2A", matlab);
1277 PrintMatrix(SType1C, N, N, LDC,
"SType1C_before_operation", matlab);
1278 PrintMatrix(SType2C, N, N, LDC,
"SType2C_before_operation", matlab);
1282 SType1BLAS.
SYRK(UPLO, TRANSA, N, K, SType1alpha, SType1A, LDA, SType1beta, SType1C, LDC);
1283 SType2BLAS.
SYRK(UPLO, TRANSA, N, K, SType2alpha, SType2A, LDA, SType2beta, SType2C, LDC);
1286 PrintMatrix(SType1C, N, N, LDC,
"SType1C_after_operation", matlab);
1287 PrintMatrix(SType2C, N, N, LDC,
"SType2C_after_operation", matlab);
1289 GoodTestSubcount +=
CompareMatrices(SType1C, SType2C, N, N, LDC, TOL);
1296 GoodTestCount += GoodTestSubcount;
1297 if(verbose || debug) std::cout <<
"SYRK: " << GoodTestSubcount <<
" of " << SYRKTESTS <<
" tests were successful." << std::endl;
1298 if(debug) std::cout << std::endl;
1306 GoodTestSubcount = 0;
1317 SType1B =
new SType1[LDB * N];
1318 SType2B =
new SType2[LDB * N];
1332 SType1A =
new SType1[LDA * M];
1333 SType2A =
new SType2[LDA * M];
1335 for(j = 0; j < M; j++)
1340 for(k = 0; k < M; k++)
1345 SType1A[j*LDA+k] = SType1zero;
1347 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1351 while (SType1A[j*LDA+k] == SType1zero) {
1354 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1356 SType1A[j*LDA+k] = SType1one;
1357 SType2A[j*LDA+k] = SType2one;
1364 for(k = 0; k < M; k++)
1369 SType1A[j*LDA+k] = SType1zero;
1371 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1375 while (SType1A[j*LDA+k] == SType1zero) {
1378 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1380 SType1A[j*LDA+k] = SType1one;
1381 SType2A[j*LDA+k] = SType2one;
1395 SType1A =
new SType1[LDA * N];
1396 SType2A =
new SType2[LDA * N];
1398 for(j = 0; j < N; j++)
1403 for(k = 0; k < N; k++)
1408 SType1A[j*LDA+k] = SType1zero;
1410 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1414 while (SType1A[j*LDA+k] == SType1zero) {
1417 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1419 SType1A[j*LDA+k] = SType1one;
1420 SType2A[j*LDA+k] = SType2one;
1427 for(k = 0; k < N; k++)
1432 SType1A[j*LDA+k] = SType1zero;
1434 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1438 while (SType1A[j*LDA+k] == SType1zero) {
1441 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1443 SType1A[j*LDA+k] = SType1one;
1444 SType2A[j*LDA+k] = SType2one;
1453 for(j = 0; j < N; j++) {
1454 for(k = 0; k < M; k++) {
1456 SType2B[j*LDB+k] =
ConvertType(SType1B[j*LDB+k], convertTo);
1460 SType2alpha =
ConvertType(SType1alpha, convertTo);
1463 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
1464 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
1465 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
1467 PrintMatrix(SType1A, M, M, LDA,
"SType1A", matlab);
1468 PrintMatrix(SType2A, M, M, LDA,
"SType2A", matlab);
1470 PrintMatrix(SType1A, N, N, LDA,
"SType1A", matlab);
1471 PrintMatrix(SType2A, N, N, LDA,
"SType2A", matlab);
1473 PrintMatrix(SType1B, M, N, LDB,
"SType1B_before_operation", matlab);
1474 PrintMatrix(SType2B, M, N, LDB,
"SType2B_before_operation", matlab);
1477 SType1BLAS.
TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
1478 SType2BLAS.
TRMM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
1481 PrintMatrix(SType1B, M, N, LDB,
"SType1B_after_operation", matlab);
1482 PrintMatrix(SType2B, M, N, LDB,
"SType2B_after_operation", matlab);
1484 GoodTestSubcount +=
CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
1490 GoodTestCount += GoodTestSubcount;
1491 if(verbose || debug) std::cout <<
"TRMM: " << GoodTestSubcount <<
" of " << TRMMTESTS <<
" tests were successful." << std::endl;
1492 if(debug) std::cout << std::endl;
1500 GoodTestSubcount = 0;
1511 SType1B =
new SType1[LDB * N];
1512 SType2B =
new SType2[LDB * N];
1529 SType1A =
new SType1[LDA * M];
1530 SType2A =
new SType2[LDA * M];
1532 for(j = 0; j < M; j++)
1537 for(k = 0; k < M; k++)
1542 SType1A[j*LDA+k] = SType1zero;
1544 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1548 while (SType1A[j*LDA+k] == SType1zero) {
1551 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1553 SType1A[j*LDA+k] = SType1one;
1554 SType2A[j*LDA+k] = SType2one;
1561 for(k = 0; k < M; k++)
1566 SType1A[j*LDA+k] = SType1zero;
1568 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1572 while (SType1A[j*LDA+k] == SType1zero) {
1575 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1577 SType1A[j*LDA+k] = SType1one;
1578 SType2A[j*LDA+k] = SType2one;
1592 SType1A =
new SType1[LDA * N];
1593 SType2A =
new SType2[LDA * N];
1595 for(j = 0; j < N; j++)
1600 for(k = 0; k < N; k++)
1605 SType1A[j*LDA+k] = SType1zero;
1607 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1611 while (SType1A[j*LDA+k] == SType1zero) {
1614 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1616 SType1A[j*LDA+k] = SType1one;
1617 SType2A[j*LDA+k] = SType2one;
1624 for(k = 0; k < N; k++)
1629 SType1A[j*LDA+k] = SType1zero;
1631 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1635 while (SType1A[j*LDA+k] == SType1zero) {
1638 SType2A[j*LDA+k] =
ConvertType(SType1A[j*LDA+k], convertTo);
1640 SType1A[j*LDA+k] = SType1one;
1641 SType2A[j*LDA+k] = SType2one;
1650 for(j = 0; j < N; j++)
1652 for(k = 0; k < M; k++)
1655 SType2B[j*LDB+k] =
ConvertType(SType1B[j*LDB+k], convertTo);
1660 SType2alpha =
ConvertType(SType1alpha, convertTo);
1664 std::cout <<
"Test #" << TotalTestCount <<
" --" << std::endl;
1669 std::cout <<
"M="<<M <<
"\t" <<
"N="<<N <<
"\t" <<
"LDA="<<LDA <<
"\t" <<
"LDB="<<LDB << std::endl;
1670 std::cout <<
"SType1alpha = " << SType1alpha << std::endl;
1671 std::cout <<
"SType2alpha = " << SType2alpha << std::endl;
1673 PrintMatrix(SType1A, M, M, LDA,
"SType1A", matlab);
1674 PrintMatrix(SType2A, M, M, LDA,
"SType2A", matlab);
1676 PrintMatrix(SType1A, N, N, LDA,
"SType1A", matlab);
1677 PrintMatrix(SType2A, N, N, LDA,
"SType2A", matlab);
1679 PrintMatrix(SType1B, M, N, LDB,
"SType1B_before_operation", matlab);
1680 PrintMatrix(SType2B, M, N, LDB,
"SType2B_before_operation", matlab);
1684 SType1BLAS.
TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType1alpha, SType1A, LDA, SType1B, LDB);
1685 SType2BLAS.
TRSM(SIDE, UPLO, TRANSA, DIAG, M, N, SType2alpha, SType2A, LDA, SType2B, LDB);
1689 PrintMatrix(SType1B, M, N, LDB,
"SType1B_after_operation", matlab);
1690 PrintMatrix(SType2B, M, N, LDB,
"SType2B_after_operation", matlab);
1694 std::cout <<
"FAILED TEST!!!!!!" << std::endl;
1695 GoodTestSubcount +=
CompareMatrices(SType1B, SType2B, M, N, LDB, TOL);
1702 GoodTestCount += GoodTestSubcount;
1703 if(verbose || debug) std::cout <<
"TRSM: " << GoodTestSubcount <<
" of " << TRSMTESTS <<
" tests were successful." << std::endl;
1704 if(debug) std::cout << std::endl;
1709 #ifdef HAVE_TEUCHOS_ARPREC
1714 if((((TotalTestCount - 1) - GoodTestCount) != 0) || (verbose) || (debug))
1716 std::cout << GoodTestCount <<
" of " << (TotalTestCount - 1) <<
" total tests were successful." << std::endl;
1719 if ((TotalTestCount-1) == GoodTestCount) {
1720 std::cout <<
"End Result: TEST PASSED" << std::endl;
1724 std::cout <<
"End Result: TEST FAILED" << std::endl;
1725 return (TotalTestCount-GoodTestCount-1);
1728 template<
typename TYPE>
1746 template<
typename TYPE>
1749 std::cout << Name <<
" =" << std::endl;
1751 if(Matlab) std::cout <<
"[";
1752 for(i = 0; i < Size; i++)
1754 std::cout << Vector[i] <<
" ";
1756 if(Matlab) std::cout <<
"]";
1759 std::cout << std::endl << std::endl;
1763 std::cout <<
";" << std::endl;
1767 template<
typename TYPE>
1772 std::cout << Name <<
" =" << std::endl;
1774 for(i = 0; i < Rows; i++)
1776 for(j = 0; j < Columns; j++)
1778 std::cout << Matrix[i + (j * LDM)] <<
" ";
1780 std::cout << std::endl;
1782 std::cout << std::endl;
1786 std::cout << Name <<
" = ";
1789 for(i = 0; i < Rows; i++)
1792 for(j = 0; j < Columns; j++)
1794 std::cout << Matrix[i + (j * LDM)] <<
" ";
1798 std::cout <<
"];" << std::endl;
1802 template<
typename TYPE1,
typename TYPE2>
1811 return( temp2 < Tolerance );
1818 template<
typename TYPE1,
typename TYPE2>
1827 for(i = 0; i < Size; i++)
1830 sum2 += temp2*temp2;
1839 if (temp2 > Tolerance)
1848 template<
typename TYPE1,
typename TYPE2>
1857 for(j = 0; j < Columns; j++)
1859 for(i = 0; i < Rows; i++)
1872 if (temp2 > Tolerance)
1879 template<
typename TYPE1,
typename TYPE2>
1882 return static_cast<TYPE2
>(T1);
1885 #ifdef HAVE_TEUCHOS_ARPREC
1893 #ifdef HAVE_TEUCHOS_QD
1897 return to_double(T1);
1901 #ifdef HAVE_TEUCHOS_GNU_MP
1943 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.
T magnitudeType
Mandatory typedef for result of magnitude.
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()