50#include <Teuchos_StandardCatchMacros.hpp>
51#include <Teuchos_TimeMonitor.hpp>
52#include <Teuchos_SerialDenseHelpers.hpp>
65 template<
class Scalar,
class MV>
69 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
magnitude_type;
71 typedef Teuchos::SerialDenseMatrix<int, Scalar>
mat_type;
93 using Teuchos::TimeMonitor;
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,
157 const Teuchos::RCP<const MV>& X,
162 std::ostream& resultStream,
163 const bool displayResultsCompactly=
false)
165 using Teuchos::Array;
166 using Teuchos::ArrayView;
170 using Teuchos::TimeMonitor;
173 TEUCHOS_TEST_FOR_EXCEPTION(orthoMan.is_null(), std::invalid_argument,
175 TEUCHOS_TEST_FOR_EXCEPTION(X.is_null(), std::invalid_argument,
177 TEUCHOS_TEST_FOR_EXCEPTION(numCols < 1, std::invalid_argument,
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) {
198 C[k] = rcp (
new mat_type (numCols, numCols));
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 >
331 typedef Teuchos::ScalarTraits<scalar_type>
SCT;
333 typedef Teuchos::ScalarTraits<magnitude_type>
SMT;
335 typedef Teuchos::SerialDenseMatrix<int, scalar_type>
mat_type;
355 const bool isRankRevealing,
356 const Teuchos::RCP<MV>& S,
361 using Teuchos::Array;
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);
408 TEUCHOS_TEST_FOR_EXCEPTION(initialX1Rank != sizeX1,
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);
416 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
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));
444 TEUCHOS_TEST_FOR_EXCEPTION(initialX2Rank != sizeX2,
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);
452 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL,
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);
459 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL,
461 "projectAndNormalize(X2,X1) did not meet tolerance: "
462 "orthogError(X2,X1) = " << err <<
" > TOL = " << TOL);
463 debugOut <<
"done: || <X2,X1> || = " << err << endl;
466#ifdef HAVE_BELOS_TSQR
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);
497 TEUCHOS_TEST_FOR_EXCEPTION(initialX1Rank != sizeX1, std::runtime_error,
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);
504 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
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));
536 TEUCHOS_TEST_FOR_EXCEPTION(initialX2Rank != sizeX2,
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);
545 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
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);
553 TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
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 ==="
572 debugOut <<
"Testing project() by projecting a random multivector S "
573 "against various combinations of X1 and X2 " << endl;
574 const int thisNumFailed =
testProject(OM,S,X1,X2,MyOM);
575 numFailed += thisNumFailed;
576 if (thisNumFailed > 0)
577 debugOut <<
" *** " << thisNumFailed
578 << (thisNumFailed > 1 ?
" tests" :
" test")
579 <<
" failed." << endl;
586 debugOut <<
"Testing normalize() on bad multivectors " << endl;
588 numFailed += thisNumFailed;
599 mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
600 Teuchos::randomSyncedMatrix(C1);
601 Teuchos::randomSyncedMatrix(C2);
607 debugOut <<
"Testing project() by projecting [X1 X2]-range multivector "
608 "against P_X1 P_X2 " << endl;
609 const int thisNumFailed =
testProject(OM,S,X1,X2,MyOM);
610 numFailed += thisNumFailed;
611 if (thisNumFailed > 0)
612 debugOut <<
" *** " << thisNumFailed
613 << (thisNumFailed > 1 ?
" tests" :
" test")
614 <<
" failed." << endl;
619 if (isRankRevealing && sizeS > 2)
625 std::vector<int> ind(1);
629 debugOut <<
"Testing normalize() on a rank-deficient multivector " << endl;
631 numFailed += thisNumFailed;
632 if (thisNumFailed > 0)
633 debugOut <<
" *** " << thisNumFailed
634 << (thisNumFailed > 1 ?
" tests" :
" test")
635 <<
" failed." << endl;
640 if (isRankRevealing && sizeS > 1)
646 Teuchos::randomSyncedMatrix(scaleS);
648 for (
int i=0; i<sizeS; i++)
650 std::vector<int> ind(1);
655 debugOut <<
"Testing normalize() on a rank-1 multivector " << endl;
657 numFailed += thisNumFailed;
658 if (thisNumFailed > 0)
659 debugOut <<
" *** " << thisNumFailed
660 << (thisNumFailed > 1 ?
" tests" :
" test")
661 <<
" failed." << endl;
665 std::vector<int> ind(1);
668 debugOut <<
"Testing projectAndNormalize() on a random multivector " << endl;
670 numFailed += thisNumFailed;
671 if (thisNumFailed > 0)
672 debugOut <<
" *** " << thisNumFailed
673 << (thisNumFailed > 1 ?
" tests" :
" test")
674 <<
" failed." << endl;
685 mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
686 Teuchos::randomSyncedMatrix(C1);
687 Teuchos::randomSyncedMatrix(C2);
691 debugOut <<
"Testing projectAndNormalize() by projecting [X1 X2]-range "
692 "multivector against P_X1 P_X2 " << endl;
694 numFailed += thisNumFailed;
695 if (thisNumFailed > 0)
696 debugOut <<
" *** " << thisNumFailed
697 << (thisNumFailed > 1 ?
" tests" :
" test")
698 <<
" failed." << endl;
703 if (isRankRevealing && sizeS > 2)
709 std::vector<int> ind(1);
713 debugOut <<
"Testing projectAndNormalize() on a rank-deficient "
714 "multivector " << endl;
716 numFailed += thisNumFailed;
717 if (thisNumFailed > 0)
718 debugOut <<
" *** " << thisNumFailed
719 << (thisNumFailed > 1 ?
" tests" :
" test")
720 <<
" failed." << endl;
725 if (isRankRevealing && sizeS > 1)
731 Teuchos::randomSyncedMatrix(scaleS);
733 for (
int i=0; i<sizeS; i++)
735 std::vector<int> ind(1);
740 debugOut <<
"Testing projectAndNormalize() on a rank-1 multivector " << endl;
741 bool constantStride =
true;
743 debugOut <<
"-- S does not have constant stride" << endl;
744 constantStride =
false;
747 debugOut <<
"-- X1 does not have constant stride" << endl;
748 constantStride =
false;
751 debugOut <<
"-- X2 does not have constant stride" << endl;
752 constantStride =
false;
754 if (! constantStride) {
755 debugOut <<
"-- Skipping this test, since TSQR does not work on "
756 "multivectors with nonconstant stride" << endl;
760 numFailed += thisNumFailed;
761 if (thisNumFailed > 0) {
762 debugOut <<
" *** " << thisNumFailed
763 << (thisNumFailed > 1 ?
" tests" :
" test")
764 <<
" failed." << endl;
769 if (numFailed != 0) {
770 MyOM->stream(
Errors) << numFailed <<
" total test failures." << endl;
790 "MVDiff: X and Y should have the same number of columns."
791 " X has " << numCols <<
" column(s) and Y has "
816 for (
int i = 0; i < numCols; ++i)
817 err += SCT::magnitude (C(i,i));
819 return SCT::magnitude (SCT::squareroot (err));
825 const Teuchos::RCP< const MV >& S,
826 const Teuchos::RCP< const MV >& X1,
827 const Teuchos::RCP< const MV >& X2,
839 const Teuchos::RCP< const MV >& S,
840 const Teuchos::RCP< const MV >& X1,
841 const Teuchos::RCP< const MV >& X2,
844 using Teuchos::Array;
848 using Teuchos::tuple;
862 std::ostringstream sout;
898 sout <<
" || <S,X1> || before : " << err << endl;
902 sout <<
" || <S,X2> || before : " << err << endl;
905 for (
int t=0; t<numtests; t++) {
907 Array< RCP< const MV > > theX;
908 RCP<mat_type > B = rcp(
new mat_type(sizeS,sizeS) );
909 Array<RCP<mat_type > > C;
910 if ( (t % 3) == 0 ) {
914 else if ( (t % 3) == 1 ) {
917 C = tuple( rcp(
new mat_type(sizeX1,sizeS)) );
919 else if ( (t % 3) == 2 ) {
922 C = tuple( rcp(
new mat_type(sizeX2,sizeS)) );
927 C = tuple( rcp(
new mat_type(sizeX1,sizeS)),
944 Array<RCP<MV> > S_outs;
945 Array<Array<RCP<mat_type > > > C_outs;
946 Array<RCP<mat_type > > B_outs;
953 Teuchos::randomSyncedMatrix(*B);
955 Teuchos::randomSyncedMatrix(*C[i]);
964 int ret = OM->projectAndNormalize(*Scopy,C,B,theX);
965 sout <<
"projectAndNormalize() returned rank " << ret << endl;
967 sout <<
" *** Error: returned rank is zero, cannot continue tests" << endl;
971 ret_out.push_back(ret);
981 std::vector<int> ind(ret);
982 for (
int i=0; i<ret; i++) {
986 B_outs.push_back( rcp(
new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
989 S_outs.push_back( Scopy );
990 B_outs.push_back( rcp(
new mat_type(*B) ) );
992 C_outs.push_back( Array<RCP<mat_type > >(0) );
994 C_outs.back().push_back( rcp(
new mat_type(*C[0]) ) );
997 C_outs.back().push_back( rcp(
new mat_type(*C[1]) ) );
1001 if ( (t % 3) == 3 ) {
1009 Teuchos::randomSyncedMatrix(*B);
1011 Teuchos::randomSyncedMatrix(*C[i]);
1014 theX = tuple( theX[1], theX[0] );
1018 ret = OM->projectAndNormalize(*Scopy,C,B,theX);
1019 sout <<
"projectAndNormalize() returned rank " << ret << endl;
1021 sout <<
" *** Error: returned rank is zero, cannot continue tests" << endl;
1025 ret_out.push_back(ret);
1035 std::vector<int> ind(ret);
1036 for (
int i=0; i<ret; i++) {
1040 B_outs.push_back( rcp(
new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
1043 S_outs.push_back( Scopy );
1044 B_outs.push_back( rcp(
new mat_type(*B) ) );
1046 C_outs.push_back( Array<RCP<mat_type > >() );
1048 C_outs.back().push_back( rcp(
new mat_type(*C[1]) ) );
1049 C_outs.back().push_back( rcp(
new mat_type(*C[0]) ) );
1051 theX = tuple( theX[1], theX[0] );
1056 for (
size_type o=0; o<S_outs.size(); o++) {
1062 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1063 <<
" total tests) failed: Tolerance exceeded! Error = "
1064 << err <<
" > TOL = " << TOL <<
"."
1068 sout <<
" || <S,S> - I || after : " << err << endl;
1074 if (C_outs[o].size() > 0) {
1076 if (C_outs[o].size() > 1) {
1081 if (err > ATOL*TOL) {
1083 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1084 <<
" total tests) failed: Tolerance exceeded! Error = "
1085 << err <<
" > ATOL*TOL = " << (ATOL*TOL) <<
"."
1089 sout <<
" " << t <<
"|| S_in - X1*C1 - X2*C2 - S_out*B || : " << err << endl;
1092 if (theX.size() > 0 && theX[0] != null) {
1096 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1097 <<
" total tests) failed: Tolerance exceeded! Error = "
1098 << err <<
" > TOL = " << TOL <<
"."
1102 sout <<
" " << t <<
"|| <X[0],S> || after : " << err << endl;
1105 if (theX.size() > 1 && theX[1] != null) {
1109 <<
" *** Test (number " << (t+1) <<
" of " << numtests
1110 <<
" total tests) failed: Tolerance exceeded! Error = "
1111 << err <<
" > TOL = " << TOL <<
"."
1115 sout <<
" " << t <<
"|| <X[1],S> || after : " << err << endl;
1120 sout <<
" *** Error: OrthoManager threw exception: " << e.what() << endl;
1131 const int msgType = (numerr > 0) ?
1132 (
static_cast<int>(
Debug) |
static_cast<int>(
Errors)) :
1133 static_cast<int>(
Debug);
1137 MyOM->stream(
static_cast< MsgType >(msgType)) << sout.str() << endl;
1147 const Teuchos::RCP< const MV >& S,
1152 int numFailures = 0;
1155 const int msgType = (
static_cast<int>(
Debug) |
static_cast<int>(
Errors));
1159 RCP< mat_type > bZero (
new mat_type (1, 1));
1160 std::vector< magnitude_type > zeroNorm( 1 );
1163 OM->normalize( *zeroVec, bZero );
1166 if ( zeroNorm[0] != ZERO )
1168 MyOM->stream(
static_cast< MsgType >(msgType)) <<
" --> Normalization of zero vector FAILED!" << std::endl;
1181 const Teuchos::RCP< const MV >& S,
1184 using Teuchos::Array;
1187 using Teuchos::tuple;
1190 std::ostringstream sout;
1235 sout <<
"The test matrix S has Frobenius norm " << ATOL
1236 <<
", and the relative error tolerance is TOL = "
1237 << TOL <<
"." << endl;
1239 const int numtests = 1;
1240 for (
int t = 0; t < numtests; ++t) {
1252 RCP< mat_type > B (
new mat_type (sizeS, sizeS));
1257 Teuchos::randomSyncedMatrix(*B);
1259 const int reportedRank = OM->normalize (*S_copy, B);
1260 sout <<
"normalize() returned rank " << reportedRank << endl;
1261 if (reportedRank == 0) {
1262 sout <<
" *** Error: Cannot continue, since normalize() "
1263 "reports that S has rank 0" << endl;
1277 std::vector<int> indices (reportedRank);
1278 for (
int j = 0; j < reportedRank; ++j)
1289 RCP< mat_type > B_top (
new mat_type (Teuchos::Copy, *B, reportedRank, sizeS));
1295 sout <<
" *** Error: Tolerance exceeded: err = "
1296 << err <<
" > TOL = " << TOL << endl;
1299 sout <<
" || <S,S> - I || after : " << err << endl;
1317 if (err > ATOL*TOL) {
1318 sout <<
" *** Error: Tolerance exceeded: err = "
1319 << err <<
" > ATOL*TOL = " << (ATOL*TOL) << endl;
1322 sout <<
" " << t <<
"|| S - Q*B || : " << err << endl;
1326 sout <<
" *** Error: the OrthoManager's normalize() method "
1327 "threw an exception: " << e.what() << endl;
1334 MyOM->stream(type) << sout.str();
1335 MyOM->stream(type) << endl;
1346 const Teuchos::RCP< const MV >& S,
1347 const Teuchos::RCP< const MV >& X1,
1348 const Teuchos::RCP< const MV >& X2,
1351 using Teuchos::Array;
1352 using Teuchos::null;
1355 using Teuchos::tuple;
1359 std::ostringstream sout;
1396 sout <<
"-- The test matrix S has Frobenius norm " << ATOL
1397 <<
", and the relative error tolerance is TOL = "
1398 << TOL <<
"." << endl;
1406 const int num_X = 2;
1407 Array< RCP< const MV > > X (num_X);
1412 RCP< mat_type > B (
new mat_type (sizeS, sizeS));
1417 Array< RCP< mat_type > > C (num_X);
1418 for (
int k = 0; k < num_X; ++k)
1421 Teuchos::randomSyncedMatrix(*C[k]);
1425 const int reportedRank = OM->projectAndNormalize (*Q, C, B, X);
1428 std::vector<int> indices (reportedRank);
1429 for (
int j = 0; j < reportedRank; ++j)
1437 sout <<
"-- ||Q(1:" << reportedRank <<
")^* Q(1:" << reportedRank
1438 <<
") - I||_F = " << orthoError << endl;
1439 if (orthoError > TOL)
1441 sout <<
" *** Error: ||Q(1:" << reportedRank <<
")^* Q(1:"
1442 << reportedRank <<
") - I||_F = " << orthoError
1443 <<
" > TOL = " << TOL <<
"." << endl;
1453 MVT::MvAddMv (SCT::one(), *S, SCT::zero(), *Residual, *Residual);
1459 RCP< const mat_type > B_top (
new mat_type (Teuchos::Copy, *B, reportedRank, B->numCols()));
1464 for (
int k = 0; k < num_X; ++k)
1467 sout <<
"-- ||S - Q(:, 1:" << reportedRank <<
")*B(1:"
1468 << reportedRank <<
", :) - X1*C1 - X2*C2||_F = "
1469 << residErr << endl;
1470 if (residErr > ATOL * TOL)
1472 sout <<
" *** Error: ||S - Q(:, 1:" << reportedRank
1473 <<
")*B(1:" << reportedRank <<
", :) "
1474 <<
"- X1*C1 - X2*C2||_F = " << residErr
1475 <<
" > ATOL*TOL = " << (ATOL*TOL) <<
"." << endl;
1480 if (reportedRank == 0)
1482 sout <<
"-- Reported rank of Q is zero: skipping Q, X[k] "
1483 "orthogonality test." << endl;
1487 for (
int k = 0; k < num_X; ++k)
1491 sout <<
"-- ||<Q(1:" << reportedRank <<
"), X[" << k
1492 <<
"]>||_F = " << projErr << endl;
1493 if (projErr > ATOL*TOL)
1495 sout <<
" *** Error: ||<Q(1:" << reportedRank <<
"), X["
1496 << k <<
"]>||_F = " << projErr <<
" > ATOL*TOL = "
1497 << (ATOL*TOL) <<
"." << endl;
1503 sout <<
" *** Error: The OrthoManager subclass instance threw "
1504 "an exception: " << e.what() << endl;
1511 MyOM->stream(type) << sout.str();
1512 MyOM->stream(type) << endl;
1523 const Teuchos::RCP< const MV >& S,
1524 const Teuchos::RCP< const MV >& X1,
1525 const Teuchos::RCP< const MV >& X2,
1528 using Teuchos::Array;
1529 using Teuchos::null;
1532 using Teuchos::tuple;
1536 std::ostringstream sout;
1573 sout <<
"The test matrix S has Frobenius norm " << ATOL
1574 <<
", and the relative error tolerance is TOL = "
1575 << TOL <<
"." << endl;
1582 const int num_X = 2;
1583 Array< RCP< const MV > > X (num_X);
1590 Array< RCP< mat_type > > C (num_X);
1591 for (
int k = 0; k < num_X; ++k)
1594 Teuchos::randomSyncedMatrix(*C[k]);
1598 OM->project(*S_copy, C, X);
1605 MVT::MvAddMv (SCT::one(), *S, -SCT::one(), *S_copy, *Residual);
1607 for (
int k = 0; k < num_X; ++k)
1610 sout <<
" ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr;
1611 if (residErr > ATOL * TOL)
1613 sout <<
" *** Error: ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr
1614 <<
" > ATOL*TOL = " << (ATOL*TOL) <<
".";
1617 for (
int k = 0; k < num_X; ++k)
1623 sout <<
" *** Error: S is not orthogonal to X[" << k
1624 <<
"] by a factor of " << projErr <<
" > TOL = "
1630 sout <<
" *** Error: The OrthoManager subclass instance threw "
1631 "an exception: " << e.what() << endl;
1638 MyOM->stream(type) << sout.str();
1639 MyOM->stream(type) << endl;
1646 const Teuchos::RCP< const MV >& S,
1647 const Teuchos::RCP< const MV >& X1,
1648 const Teuchos::RCP< const MV >& X2,
1659 const Teuchos::RCP< const MV >& S,
1660 const Teuchos::RCP< const MV >& X1,
1661 const Teuchos::RCP< const MV >& X2,
1664 using Teuchos::Array;
1665 using Teuchos::null;
1668 using Teuchos::tuple;
1673 std::ostringstream sout;
1712 sout <<
"The test matrix S has Frobenius norm " << ATOL
1713 <<
", and the relative error tolerance is TOL = "
1714 << TOL <<
"." << endl;
1750 sout <<
" || <S,X1> || before : " << err << endl;
1754 sout <<
" || <S,X2> || before : " << err << endl;
1757 for (
int t = 0; t < numtests; ++t)
1759 Array< RCP< const MV > > theX;
1760 Array< RCP< mat_type > > C;
1761 if ( (t % 3) == 0 ) {
1765 else if ( (t % 3) == 1 ) {
1768 C = tuple( rcp(
new mat_type(sizeX1,sizeS)) );
1770 else if ( (t % 3) == 2 ) {
1773 C = tuple( rcp(
new mat_type(sizeX2,sizeS)) );
1777 theX = tuple(X1,X2);
1778 C = tuple( rcp(
new mat_type(sizeX1,sizeS)),
1791 Array< RCP< MV > > S_outs;
1792 Array< Array< RCP< mat_type > > > C_outs;
1798 for (
size_type i = 0; i < C.size(); ++i) {
1799 Teuchos::randomSyncedMatrix(*C[i]);
1804 OM->project(*Scopy,C,theX);
1807 S_outs.push_back( Scopy );
1808 C_outs.push_back( Array< RCP< mat_type > >(0) );
1810 C_outs.back().push_back( rcp(
new mat_type(*C[0]) ) );
1813 C_outs.back().push_back( rcp(
new mat_type(*C[1]) ) );
1817 if ( (t % 3) == 3 ) {
1821 for (
size_type i = 0; i < C.size(); ++i) {
1822 Teuchos::randomSyncedMatrix(*C[i]);
1825 theX = tuple( theX[1], theX[0] );
1829 OM->project(*Scopy,C,theX);
1832 S_outs.push_back( Scopy );
1835 C_outs.push_back( Array<RCP<mat_type > >() );
1837 C_outs.back().push_back( rcp(
new mat_type(*C[1]) ) );
1838 C_outs.back().push_back( rcp(
new mat_type(*C[0]) ) );
1840 theX = tuple( theX[1], theX[0] );
1844 for (
size_type o = 0; o < S_outs.size(); ++o) {
1848 if (C_outs[o].size() > 0) {
1850 if (C_outs[o].size() > 1) {
1855 if (err > ATOL*TOL) {
1856 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1859 sout <<
" " << t <<
"|| S_in - X1*C1 - X2*C2 - S_out || : " << err << endl;
1862 if (theX.size() > 0 && theX[0] != null) {
1865 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1868 sout <<
" " << t <<
"|| <X[0],S> || after : " << err << endl;
1871 if (theX.size() > 1 && theX[1] != null) {
1874 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1877 sout <<
" " << t <<
"|| <X[1],S> || after : " << err << endl;
1886 for (
size_type o1=0; o1<S_outs.size(); o1++) {
1887 for (
size_type o2=o1+1; o2<S_outs.size(); o2++) {
1895 sout <<
" vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1903 sout <<
" ------------------------------------------- project() threw exception" << endl;
1904 sout <<
" Error: " << e.what() << endl;
1910 if (numerr>0) type =
Errors;
1911 MyOM->stream(type) << sout.str();
1912 MyOM->stream(type) << endl;
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Class which manages the output and verbosity of the Belos solvers.
Traits class which defines basic operations on multivectors.
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 .
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ,...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
static bool HasConstantStride(const MV &mv)
Whether the given multivector mv has constant stride.
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).
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 void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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 int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
Exception thrown to signal error in an orthogonalization manager method.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Mixin for out-of-place orthogonalization.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
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.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
MultiVecTraits< Scalar, MV > MVT
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
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.
Wrapper around OrthoManager test functionality.
static int testProjectAndNormalizeOld(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > &OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::projectAndNormalize() for the specific OrthoManager instance.
static int testProjectOld(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::project() for the specific OrthoManager instance.
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 int testProjectAndNormalizeNew(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::projectAndNormalize() for the specific OrthoManager instance.
static int testProject(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
static int testNormalizeRankReveal(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > &OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::normalize() for the specific OrthoManager instance.
Teuchos::Array< Teuchos::RCP< MV > >::size_type size_type
static int testProjectNew(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::project() for the specific OrthoManager instance.
Teuchos::ScalarTraits< scalar_type > SCT
static int testNormalize(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > &OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
Test OrthoManager::normalize() for the specific OrthoManager instance.
Teuchos::ScalarTraits< magnitude_type > SMT
static magnitude_type MVDiff(const MV &X, const MV &Y)
Compute and return $|X - Y|_F$, the Frobenius (sum of squares) norm of the difference between X and Y...
Teuchos::SerialDenseMatrix< int, scalar_type > mat_type
MultiVecTraits< scalar_type, MV > MVT
static int testProjectAndNormalize(const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, const Teuchos::RCP< const MV > &S, const Teuchos::RCP< const MV > &X1, const Teuchos::RCP< const MV > &X2, const Teuchos::RCP< Belos::OutputManager< Scalar > > &MyOM)
static magnitude_type frobeniusNorm(const MV &X)
Compute and return the Frobenius norm of X.
SCT::magnitudeType magnitude_type
MsgType
Available message types recognized by the linear solvers.