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()