327 typedef typename Teuchos::Array<Teuchos::RCP<MV> >::size_type size_type;
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;
587 const int thisNumFailed = testNormalize(OM,S,MyOM);
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;
630 const int thisNumFailed = testNormalizeRankReveal(OM,S,MyOM);
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;
656 const int thisNumFailed = testNormalizeRankReveal(OM,S,MyOM);
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;
669 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
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;
693 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
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;
715 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
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;
759 const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
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;
782 MVDiff (
const MV& X,
const MV& Y)
790 "MVDiff: X and Y should have the same number of columns."
791 " X has " << numCols <<
" column(s) and Y has "
798 return frobeniusNorm (*Resid);
806 frobeniusNorm (
const MV& X)
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,
830 return testProjectAndNormalizeNew (OM, S, X1, X2, MyOM);
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);
954 for (size_type i=0; i<C.size(); i++) {
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);
1010 for (size_type i=0; i<C.size(); i++) {
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,
1651 return testProjectNew (OM, S, X1, X2, MyOM);
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;