50 #include <Teuchos_StandardCatchMacros.hpp>
65 template<
class Scalar,
class MV>
68 typedef Scalar scalar_type;
98 Array<RCP<MV> > V (numBlocks);
99 for (
int k = 0; k < numBlocks; ++k) {
105 RCP<Time> timer = TimeMonitor::getNewCounter (
"Baseline for OrthoManager benchmark");
111 TimeMonitor monitor (*timer);
112 for (
int trial = 0; trial < numTrials; ++trial) {
113 for (
int k = 0; k < numBlocks; ++k) {
114 for (
int j = 0; j < k; ++j)
155 const std::string& orthoManName,
156 const std::string& normalization,
162 std::ostream& resultStream,
163 const bool displayResultsCompactly=
false)
178 "numCols = " << numCols <<
" < 1");
179 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 1, std::invalid_argument,
180 "numBlocks = " << numBlocks <<
" < 1");
181 TEUCHOS_TEST_FOR_EXCEPTION(numTrials < 1, std::invalid_argument,
182 "numTrials = " << numTrials <<
" < 1");
184 std::ostream& debugOut = outMan->stream(
Debug);
196 Array<RCP<mat_type> > C (numBlocks);
197 for (
int k = 0; k < numBlocks; ++k) {
200 RCP<mat_type> B (
new mat_type (numCols, numCols));
206 Array<RCP<MV> > V (numBlocks);
207 for (
int k = 0; k < numBlocks; ++k) {
215 RCP<Time> firstRunTimer;
217 std::ostringstream os;
218 os <<
"OrthoManager: " << orthoManName <<
" first run";
219 firstRunTimer = TimeMonitor::getNewCounter (os.str());
223 std::ostringstream os;
224 os <<
"OrthoManager: " << orthoManName <<
" total over "
225 << numTrials <<
" trials (excluding first run above)";
226 timer = TimeMonitor::getNewCounter (os.str());
232 TimeMonitor monitor (*firstRunTimer);
234 (void) orthoMan->normalize (*V[0], B);
235 for (
int k = 1; k < numBlocks; ++k) {
241 ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
242 ArrayView<RCP<const MV> > V_0k =
243 Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
244 (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
256 debugOut <<
"Orthogonality of V[0:" << (numBlocks-1)
258 for (
int k = 0; k < numBlocks; ++k) {
260 debugOut <<
"For block V[" << k <<
"]:" << endl;
261 debugOut <<
" ||<V[" << k <<
"], V[" << k <<
"]> - I|| = "
262 << orthoMan->orthonormError(*V[k]) << endl;
264 for (
int j = 0; j < k; ++j) {
265 debugOut <<
" ||< V[" << j <<
"], V[" << k <<
"] >|| = "
266 << orthoMan->orthogError(*V[j], *V[k]) << endl;
274 TimeMonitor monitor (*timer);
276 for (
int trial = 0; trial < numTrials; ++trial) {
277 (void) orthoMan->normalize (*V[0], B);
278 for (
int k = 1; k < numBlocks; ++k) {
279 ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
280 ArrayView<RCP<const MV> > V_0k =
281 Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
282 (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
288 if (displayResultsCompactly)
296 resultStream <<
"#orthoManName"
301 <<
",firstRunTimeInSeconds"
305 resultStream << orthoManName
306 <<
"," << (orthoManName==
"Simple" ? normalization :
"N/A")
310 <<
"," << firstRunTimer->totalElapsedTime()
311 <<
"," << timer->totalElapsedTime()
316 TimeMonitor::summarize (resultStream);
324 template<
class Scalar,
class MV >
355 const bool isRankRevealing,
365 using Teuchos::rcp_dynamic_cast;
366 using Teuchos::tuple;
380 std::ostream& debugOut = MyOM->stream(
Debug);
389 debugOut <<
"Generating X1,X2 for testing... ";
392 debugOut <<
"done." << endl;
399 debugOut <<
"Filling X1 with random values... ";
401 debugOut <<
"done." << endl
402 <<
"Calling normalize() on X1... ";
407 const int initialX1Rank = OM->normalize(*X1, Teuchos::null);
410 "normalize(X1) returned rank "
411 << initialX1Rank <<
" from " << sizeX1
412 <<
" vectors. Cannot continue.");
413 debugOut <<
"done." << endl
414 <<
"Calling orthonormError() on X1... ";
415 err = OM->orthonormError(*X1);
417 "After normalize(X1), orthonormError(X1) = "
418 << err <<
" > TOL = " << TOL);
419 debugOut <<
"done: ||<X1,X1> - I|| = " << err << endl;
425 debugOut <<
"Filling X2 with random values... ";
427 debugOut <<
"done." << endl
428 <<
"Calling projectAndNormalize(X2, C, B, tuple(X1))... "
439 Array<RCP<mat_type> > C (1);
440 RCP<mat_type> B = Teuchos::null;
442 OM->projectAndNormalize (*X2, C, B, tuple<RCP<const MV> >(X1));
446 "projectAndNormalize(X2,X1) returned rank "
447 << initialX2Rank <<
" from " << sizeX2
448 <<
" vectors. Cannot continue.");
449 debugOut <<
"done." << endl
450 <<
"Calling orthonormError() on X2... ";
451 err = OM->orthonormError (*X2);
454 "projectAndNormalize(X2,X1) did not meet tolerance: "
455 "orthonormError(X2) = " << err <<
" > TOL = " << TOL);
456 debugOut <<
"done: || <X2,X2> - I || = " << err << endl
457 <<
"Calling orthogError(X2, X1)... ";
458 err = OM->orthogError (*X2, *X1);
461 "projectAndNormalize(X2,X1) did not meet tolerance: "
462 "orthogError(X2,X1) = " << err <<
" > TOL = " << TOL);
463 debugOut <<
"done: || <X2,X1> || = " << err << endl;
472 RCP<mixin_type> tsqr = rcp_dynamic_cast<mixin_type>(OM);
473 if (! tsqr.is_null())
477 <<
"=== OutOfPlaceNormalizerMixin tests ==="
487 debugOut <<
"Filling X1_in with random values... ";
489 debugOut <<
"done." << endl;
490 debugOut <<
"Filling X1_out with different random values...";
493 debugOut <<
"done." << endl
494 <<
"Calling normalizeOutOfPlace(*X1_in, *X1_out, null)... ";
495 const int initialX1Rank =
496 tsqr->normalizeOutOfPlace(*X1_in, *X1_out, Teuchos::null);
498 "normalizeOutOfPlace(*X1_in, *X1_out, null) "
499 "returned rank " << initialX1Rank <<
" from "
500 << sizeX1 <<
" vectors. Cannot continue.");
501 debugOut <<
"done." << endl
502 <<
"Calling orthonormError() on X1_out... ";
503 err = OM->orthonormError(*X1_out);
505 "After calling normalizeOutOfPlace(*X1_in, "
506 "*X1_out, null), orthonormError(X1) = "
507 << err <<
" > TOL = " << TOL);
508 debugOut <<
"done: ||<X1_out,X1_out> - I|| = " << err << endl;
519 debugOut <<
"Filling X2_in with random values... ";
521 debugOut <<
"done." << endl
522 <<
"Filling X2_out with different random values...";
525 debugOut <<
"done." << endl
526 <<
"Calling projectAndNormalizeOutOfPlace(X2_in, X2_out, "
527 <<
"C, B, X1_out)...";
530 Array<RCP<mat_type> > C (1);
531 RCP<mat_type> B = Teuchos::null;
533 tsqr->projectAndNormalizeOutOfPlace (*X2_in, *X2_out, C, B,
534 tuple<RCP<const MV> >(X1_out));
538 "projectAndNormalizeOutOfPlace(*X2_in, "
539 "*X2_out, C, B, tuple(X1_out)) returned rank "
540 << initialX2Rank <<
" from " << sizeX2
541 <<
" vectors. Cannot continue.");
542 debugOut <<
"done." << endl
543 <<
"Calling orthonormError() on X2_out... ";
544 err = OM->orthonormError (*X2_out);
546 "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
547 "C, B, tuple(X1_out)) did not meet tolerance: "
548 "orthonormError(X2_out) = "
549 << err <<
" > TOL = " << TOL);
550 debugOut <<
"done: || <X2_out,X2_out> - I || = " << err << endl
551 <<
"Calling orthogError(X2_out, X1_out)... ";
552 err = OM->orthogError (*X2_out, *X1_out);
554 "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
555 "C, B, tuple(X1_out)) did not meet tolerance: "
556 "orthogError(X2_out, X1_out) = "
557 << err <<
" > TOL = " << TOL);
558 debugOut <<
"done: || <X2_out,X1_out> || = " << err << endl;
560 <<
"=== Done with OutOfPlaceNormalizerMixin tests ==="
571 debugOut <<
"Testing project() by projecting a random multivector S "
572 "against various combinations of X1 and X2 " << endl;
573 const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
574 numFailed += thisNumFailed;
575 if (thisNumFailed > 0)
576 debugOut <<
" *** " << thisNumFailed
577 << (thisNumFailed > 1 ?
" tests" :
" test")
578 <<
" failed." << endl;
589 mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
590 Teuchos::randomSyncedMatrix(C1);
591 Teuchos::randomSyncedMatrix(C2);
597 debugOut <<
"Testing project() by projecting [X1 X2]-range multivector "
598 "against P_X1 P_X2 " << endl;
599 const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
600 numFailed += thisNumFailed;
601 if (thisNumFailed > 0)
602 debugOut <<
" *** " << thisNumFailed
603 << (thisNumFailed > 1 ?
" tests" :
" test")
604 <<
" failed." << endl;
609 if (isRankRevealing && sizeS > 2)
615 std::vector<int> ind(1);
619 debugOut <<
"Testing normalize() on a rank-deficient multivector " << endl;
620 const int thisNumFailed = testNormalize(OM,S,MyOM);
621 numFailed += thisNumFailed;
622 if (thisNumFailed > 0)
623 debugOut <<
" *** " << thisNumFailed
624 << (thisNumFailed > 1 ?
" tests" :
" test")
625 <<
" failed." << endl;
630 if (isRankRevealing && sizeS > 1)
636 Teuchos::randomSyncedMatrix(scaleS);
638 for (
int i=0; i<sizeS; i++)
640 std::vector<int> ind(1);
645 debugOut <<
"Testing normalize() on a rank-1 multivector " << endl;
646 const int thisNumFailed = testNormalize(OM,S,MyOM);
647 numFailed += thisNumFailed;
648 if (thisNumFailed > 0)
649 debugOut <<
" *** " << thisNumFailed
650 << (thisNumFailed > 1 ?
" tests" :
" test")
651 <<
" failed." << endl;
655 std::vector<int> ind(1);
658 debugOut <<
"Testing projectAndNormalize() on a random multivector " << endl;
659 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
660 numFailed += thisNumFailed;
661 if (thisNumFailed > 0)
662 debugOut <<
" *** " << thisNumFailed
663 << (thisNumFailed > 1 ?
" tests" :
" test")
664 <<
" failed." << endl;
675 mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
676 Teuchos::randomSyncedMatrix(C1);
677 Teuchos::randomSyncedMatrix(C2);
681 debugOut <<
"Testing projectAndNormalize() by projecting [X1 X2]-range "
682 "multivector against P_X1 P_X2 " << endl;
683 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
684 numFailed += thisNumFailed;
685 if (thisNumFailed > 0)
686 debugOut <<
" *** " << thisNumFailed
687 << (thisNumFailed > 1 ?
" tests" :
" test")
688 <<
" failed." << endl;
693 if (isRankRevealing && sizeS > 2)
699 std::vector<int> ind(1);
703 debugOut <<
"Testing projectAndNormalize() on a rank-deficient "
704 "multivector " << endl;
705 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
706 numFailed += thisNumFailed;
707 if (thisNumFailed > 0)
708 debugOut <<
" *** " << thisNumFailed
709 << (thisNumFailed > 1 ?
" tests" :
" test")
710 <<
" failed." << endl;
715 if (isRankRevealing && sizeS > 1)
721 Teuchos::randomSyncedMatrix(scaleS);
723 for (
int i=0; i<sizeS; i++)
725 std::vector<int> ind(1);
730 debugOut <<
"Testing projectAndNormalize() on a rank-1 multivector " << endl;
731 bool constantStride =
true;
733 debugOut <<
"-- S does not have constant stride" << endl;
734 constantStride =
false;
737 debugOut <<
"-- X1 does not have constant stride" << endl;
738 constantStride =
false;
741 debugOut <<
"-- X2 does not have constant stride" << endl;
742 constantStride =
false;
744 if (! constantStride) {
745 debugOut <<
"-- Skipping this test, since TSQR does not work on "
746 "multivectors with nonconstant stride" << endl;
749 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
750 numFailed += thisNumFailed;
751 if (thisNumFailed > 0) {
752 debugOut <<
" *** " << thisNumFailed
753 << (thisNumFailed > 1 ?
" tests" :
" test")
754 <<
" failed." << endl;
759 if (numFailed != 0) {
760 MyOM->stream(
Errors) << numFailed <<
" total test failures." << endl;
772 MVDiff (
const MV& X,
const MV& Y)
780 "MVDiff: X and Y should have the same number of columns."
781 " X has " << numCols <<
" column(s) and Y has "
788 return frobeniusNorm (*Resid);
796 frobeniusNorm (
const MV& X)
806 for (
int i = 0; i < numCols; ++i)
820 return testProjectAndNormalizeNew (OM, S, X1, X2, MyOM);
838 using Teuchos::tuple;
852 std::ostringstream sout;
888 sout <<
" || <S,X1> || before : " << err << endl;
892 sout <<
" || <S,X2> || before : " << err << endl;
895 for (
int t=0; t<numtests; t++) {
897 Array< RCP< const MV > > theX;
898 RCP<mat_type > B =
rcp(
new mat_type(sizeS,sizeS) );
899 Array<RCP<mat_type > > C;
900 if ( (t % 3) == 0 ) {
904 else if ( (t % 3) == 1 ) {
909 else if ( (t % 3) == 2 ) {
934 Array<RCP<MV> > S_outs;
935 Array<Array<RCP<mat_type > > > C_outs;
936 Array<RCP<mat_type > > B_outs;
943 Teuchos::randomSyncedMatrix(B);
944 for (size_type i=0; i<C.size(); i++) {
945 Teuchos::randomSyncedmatrix(*C[i]);
954 int ret = OM->projectAndNormalize(*Scopy,C,B,theX);
955 sout <<
"projectAndNormalize() returned rank " << ret << endl;
957 sout <<
" *** Error: returned rank is zero, cannot continue tests" << endl;
961 ret_out.push_back(ret);
971 std::vector<int> ind(ret);
972 for (
int i=0; i<ret; i++) {
979 S_outs.push_back( Scopy );
982 C_outs.push_back( Array<RCP<mat_type > >(0) );
984 C_outs.back().push_back(
rcp(
new mat_type(*C[0]) ) );
987 C_outs.back().push_back(
rcp(
new mat_type(*C[1]) ) );
991 if ( (t % 3) == 3 ) {
999 Teuchos::randomSyncedMatrix(B);
1000 for (size_type i=0; i<C.size(); i++) {
1001 Teuchos::randomSyncedMatrix(*C[i]);
1004 theX = tuple( theX[1], theX[0] );
1008 ret = OM->projectAndNormalize(*Scopy,C,B,theX);
1009 sout <<
"projectAndNormalize() returned rank " << ret << endl;
1011 sout <<
" *** Error: returned rank is zero, cannot continue tests" << endl;
1015 ret_out.push_back(ret);
1025 std::vector<int> ind(ret);
1026 for (
int i=0; i<ret; i++) {
1033 S_outs.push_back( Scopy );
1036 C_outs.push_back( Array<RCP<mat_type > >() );
1038 C_outs.back().push_back(
rcp(
new mat_type(*C[1]) ) );
1039 C_outs.back().push_back(
rcp(
new mat_type(*C[0]) ) );
1041 theX = tuple( theX[1], theX[0] );
1046 for (size_type o=0; o<S_outs.size(); o++) {
1052 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1053 <<
" total tests) failed: Tolerance exceeded! Error = "
1054 << err <<
" > TOL = " << TOL <<
"."
1058 sout <<
" || <S,S> - I || after : " << err << endl;
1064 if (C_outs[o].size() > 0) {
1066 if (C_outs[o].size() > 1) {
1071 if (err > ATOL*TOL) {
1073 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1074 <<
" total tests) failed: Tolerance exceeded! Error = "
1075 << err <<
" > ATOL*TOL = " << (ATOL*TOL) <<
"."
1079 sout <<
" " << t <<
"|| S_in - X1*C1 - X2*C2 - S_out*B || : " << err << endl;
1082 if (theX.size() > 0 && theX[0] != null) {
1086 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1087 <<
" total tests) failed: Tolerance exceeded! Error = "
1088 << err <<
" > TOL = " << TOL <<
"."
1092 sout <<
" " << t <<
"|| <X[0],S> || after : " << err << endl;
1095 if (theX.size() > 1 && theX[1] != null) {
1099 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1100 <<
" total tests) failed: Tolerance exceeded! Error = "
1101 << err <<
" > TOL = " << TOL <<
"."
1105 sout <<
" " << t <<
"|| <X[1],S> || after : " << err << endl;
1110 sout <<
" *** Error: OrthoManager threw exception: " << e.what() << endl;
1121 const int msgType = (numerr > 0) ?
1122 (static_cast<int>(
Debug) |
static_cast<int>(
Errors)) :
1123 static_cast<int>(
Debug);
1127 MyOM->stream(static_cast< MsgType >(msgType)) << sout.str() << endl;
1143 using Teuchos::tuple;
1146 std::ostringstream sout;
1191 sout <<
"The test matrix S has Frobenius norm " << ATOL
1192 <<
", and the relative error tolerance is TOL = "
1193 << TOL <<
"." << endl;
1195 const int numtests = 1;
1196 for (
int t = 0; t < numtests; ++t) {
1208 RCP< mat_type > B (
new mat_type (sizeS, sizeS));
1213 Teuchos::randomSyncedMatrix(B);
1215 const int reportedRank = OM->normalize (*S_copy, B);
1216 sout <<
"normalize() returned rank " << reportedRank << endl;
1217 if (reportedRank == 0) {
1218 sout <<
" *** Error: Cannot continue, since normalize() "
1219 "reports that S has rank 0" << endl;
1233 std::vector<int> indices (reportedRank);
1234 for (
int j = 0; j < reportedRank; ++j)
1251 sout <<
" *** Error: Tolerance exceeded: err = "
1252 << err <<
" > TOL = " << TOL << endl;
1255 sout <<
" || <S,S> - I || after : " << err << endl;
1273 if (err > ATOL*TOL) {
1274 sout <<
" *** Error: Tolerance exceeded: err = "
1275 << err <<
" > ATOL*TOL = " << (ATOL*TOL) << endl;
1278 sout <<
" " << t <<
"|| S - Q*B || : " << err << endl;
1282 sout <<
" *** Error: the OrthoManager's normalize() method "
1283 "threw an exception: " << e.what() << endl;
1289 const MsgType type = (numerr == 0) ?
Debug : static_cast<MsgType> (static_cast<int>(
Errors) |
static_cast<int>(
Debug));
1290 MyOM->stream(type) << sout.str();
1291 MyOM->stream(type) << endl;
1308 using Teuchos::null;
1311 using Teuchos::tuple;
1315 std::ostringstream sout;
1352 sout <<
"-- The test matrix S has Frobenius norm " << ATOL
1353 <<
", and the relative error tolerance is TOL = "
1354 << TOL <<
"." << endl;
1362 const int num_X = 2;
1363 Array< RCP< const MV > > X (num_X);
1368 RCP< mat_type > B (
new mat_type (sizeS, sizeS));
1373 Array< RCP< mat_type > > C (num_X);
1374 for (
int k = 0; k < num_X; ++k)
1377 Teuchos::randomSyncedMatrix(*C[k]);
1381 const int reportedRank = OM->projectAndNormalize (*Q, C, B, X);
1384 std::vector<int> indices (reportedRank);
1385 for (
int j = 0; j < reportedRank; ++j)
1393 sout <<
"-- ||Q(1:" << reportedRank <<
")^* Q(1:" << reportedRank
1394 <<
") - I||_F = " << orthoError << endl;
1395 if (orthoError > TOL)
1397 sout <<
" *** Error: ||Q(1:" << reportedRank <<
")^* Q(1:"
1398 << reportedRank <<
") - I||_F = " << orthoError
1399 <<
" > TOL = " << TOL <<
"." << endl;
1420 for (
int k = 0; k < num_X; ++k)
1423 sout <<
"-- ||S - Q(:, 1:" << reportedRank <<
")*B(1:"
1424 << reportedRank <<
", :) - X1*C1 - X2*C2||_F = "
1425 << residErr << endl;
1426 if (residErr > ATOL * TOL)
1428 sout <<
" *** Error: ||S - Q(:, 1:" << reportedRank
1429 <<
")*B(1:" << reportedRank <<
", :) "
1430 <<
"- X1*C1 - X2*C2||_F = " << residErr
1431 <<
" > ATOL*TOL = " << (ATOL*TOL) <<
"." << endl;
1436 if (reportedRank == 0)
1438 sout <<
"-- Reported rank of Q is zero: skipping Q, X[k] "
1439 "orthogonality test." << endl;
1443 for (
int k = 0; k < num_X; ++k)
1447 sout <<
"-- ||<Q(1:" << reportedRank <<
"), X[" << k
1448 <<
"]>||_F = " << projErr << endl;
1449 if (projErr > ATOL*TOL)
1451 sout <<
" *** Error: ||<Q(1:" << reportedRank <<
"), X["
1452 << k <<
"]>||_F = " << projErr <<
" > ATOL*TOL = "
1453 << (ATOL*TOL) <<
"." << endl;
1459 sout <<
" *** Error: The OrthoManager subclass instance threw "
1460 "an exception: " << e.what() << endl;
1466 const MsgType type = (numerr == 0) ?
Debug : static_cast<MsgType> (static_cast<int>(
Errors) |
static_cast<int>(
Debug));
1467 MyOM->stream(type) << sout.str();
1468 MyOM->stream(type) << endl;
1485 using Teuchos::null;
1488 using Teuchos::tuple;
1492 std::ostringstream sout;
1529 sout <<
"The test matrix S has Frobenius norm " << ATOL
1530 <<
", and the relative error tolerance is TOL = "
1531 << TOL <<
"." << endl;
1538 const int num_X = 2;
1539 Array< RCP< const MV > > X (num_X);
1546 Array< RCP< mat_type > > C (num_X);
1547 for (
int k = 0; k < num_X; ++k)
1550 Teuchos::randomSyncedMatrix(*C[k]);
1554 OM->project(*S_copy, C, X);
1563 for (
int k = 0; k < num_X; ++k)
1566 sout <<
" ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr;
1567 if (residErr > ATOL * TOL)
1569 sout <<
" *** Error: ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr
1570 <<
" > ATOL*TOL = " << (ATOL*TOL) <<
".";
1573 for (
int k = 0; k < num_X; ++k)
1579 sout <<
" *** Error: S is not orthogonal to X[" << k
1580 <<
"] by a factor of " << projErr <<
" > TOL = "
1586 sout <<
" *** Error: The OrthoManager subclass instance threw "
1587 "an exception: " << e.what() << endl;
1593 const MsgType type = (numerr == 0) ?
Debug : static_cast<MsgType> (static_cast<int>(
Errors) |
static_cast<int>(
Debug));
1594 MyOM->stream(type) << sout.str();
1595 MyOM->stream(type) << endl;
1607 return testProjectNew (OM, S, X1, X2, MyOM);
1621 using Teuchos::null;
1624 using Teuchos::tuple;
1629 std::ostringstream sout;
1668 sout <<
"The test matrix S has Frobenius norm " << ATOL
1669 <<
", and the relative error tolerance is TOL = "
1670 << TOL <<
"." << endl;
1706 sout <<
" || <S,X1> || before : " << err << endl;
1710 sout <<
" || <S,X2> || before : " << err << endl;
1713 for (
int t = 0; t < numtests; ++t)
1715 Array< RCP< const MV > > theX;
1716 Array< RCP< mat_type > > C;
1717 if ( (t % 3) == 0 ) {
1721 else if ( (t % 3) == 1 ) {
1726 else if ( (t % 3) == 2 ) {
1733 theX = tuple(X1,X2);
1747 Array< RCP< MV > > S_outs;
1748 Array< Array< RCP< mat_type > > > C_outs;
1754 for (size_type i = 0; i < C.size(); ++i) {
1755 Teuchos::randomSyncedMatrix(*C[i]);
1760 OM->project(*Scopy,C,theX);
1763 S_outs.push_back( Scopy );
1764 C_outs.push_back( Array< RCP< mat_type > >(0) );
1766 C_outs.back().push_back(
rcp(
new mat_type(*C[0]) ) );
1769 C_outs.back().push_back(
rcp(
new mat_type(*C[1]) ) );
1773 if ( (t % 3) == 3 ) {
1777 for (size_type i = 0; i < C.size(); ++i) {
1778 Teuchos::randomSyncedMatrix(*C[i]);
1781 theX = tuple( theX[1], theX[0] );
1785 OM->project(*Scopy,C,theX);
1788 S_outs.push_back( Scopy );
1791 C_outs.push_back( Array<RCP<mat_type > >() );
1793 C_outs.back().push_back(
rcp(
new mat_type(*C[1]) ) );
1794 C_outs.back().push_back(
rcp(
new mat_type(*C[0]) ) );
1796 theX = tuple( theX[1], theX[0] );
1800 for (size_type o = 0; o < S_outs.size(); ++o) {
1804 if (C_outs[o].size() > 0) {
1806 if (C_outs[o].size() > 1) {
1811 if (err > ATOL*TOL) {
1812 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1815 sout <<
" " << t <<
"|| S_in - X1*C1 - X2*C2 - S_out || : " << err << endl;
1818 if (theX.size() > 0 && theX[0] != null) {
1821 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1824 sout <<
" " << t <<
"|| <X[0],S> || after : " << err << endl;
1827 if (theX.size() > 1 && theX[1] != null) {
1830 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1833 sout <<
" " << t <<
"|| <X[1],S> || after : " << err << endl;
1842 for (size_type o1=0; o1<S_outs.size(); o1++) {
1843 for (size_type o2=o1+1; o2<S_outs.size(); o2++) {
1851 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1859 sout <<
" ------------------------------------------- project() threw exception" << endl;
1860 sout <<
" Error: " << e.what() << endl;
1866 if (numerr>0) type =
Errors;
1867 MyOM->stream(type) << sout.str();
1868 MyOM->stream(type) << endl;
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
static void baseline(const Teuchos::RCP< const MV > &X, const int numCols, const int numBlocks, const int numTrials)
Establish baseline run time for OrthoManager benchmark.
static magnitudeType eps()
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Mixin for out-of-place orthogonalization.
MsgType
Available message types recognized by the linear solvers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
Teuchos::SerialDenseMatrix< int, scalar_type > mat_type
MultiVecTraits< scalar_type, MV > MVT
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
Teuchos::ScalarTraits< scalar_type > SCT
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
static void benchmark(const Teuchos::RCP< OrthoManager< Scalar, MV > > &orthoMan, const std::string &orthoManName, const std::string &normalization, const Teuchos::RCP< const MV > &X, const int numCols, const int numBlocks, const int numTrials, const Teuchos::RCP< OutputManager< Scalar > > &outMan, std::ostream &resultStream, const bool displayResultsCompactly=false)
Benchmark the given orthogonalization manager.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
SCT::magnitudeType magnitude_type
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Exception thrown to signal error in an orthogonalization manager method.
static int runTests(const Teuchos::RCP< OrthoManager< Scalar, MV > > &OM, const bool isRankRevealing, const Teuchos::RCP< MV > &S, const int sizeX1, const int sizeX2, const Teuchos::RCP< OutputManager< Scalar > > &MyOM)
Run all the tests.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static bool HasConstantStride(const MV &mv)
Whether the given multivector mv has constant stride.
static magnitudeType magnitude(T a)
Teuchos::ScalarTraits< magnitude_type > SMT
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Wrapper around OrthoManager test functionality.