42 #ifndef BELOS_MVOPTESTER_HPP 43 #define BELOS_MVOPTESTER_HPP 62 #include "Teuchos_RCP.hpp" 63 #include "Teuchos_SetScientific.hpp" 83 template<
class ScalarType,
class MV >
86 const Teuchos::RCP<const MV> &A)
88 using Teuchos::SetScientific;
91 typedef Teuchos::ScalarTraits<ScalarType> STS;
92 typedef typename STS::magnitudeType MagType;
96 SetScientific<ScalarType> sci (om->stream (
Warnings));
101 const MagType
tol = Teuchos::as<MagType> (100) * STS::eps ();
177 const ScalarType one = STS::one();
178 const ScalarType zero = STS::zero();
179 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero();
182 const int numvecs = 10;
183 const int numvecs_2 = 5;
185 std::vector<int> ind(numvecs_2);
195 TEUCHOS_TEST_FOR_EXCEPT(numvecs_2 != 5);
206 if ( MVT::GetNumberVecs(*A) <= 0 ) {
208 <<
"*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl
209 <<
"Returned <= 0." << endl;
218 if ( MVT::GetGlobalLength(*A) <= 0 ) {
220 <<
"*** ERROR *** MultiVectorTraitsExt::GetGlobalLength()" << endl
221 <<
"Returned <= 0." << endl;
234 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
235 std::vector<MagType> norms(2*numvecs);
236 bool ResizeWarning =
false;
237 if ( MVT::GetNumberVecs(*B) != numvecs ) {
239 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
240 <<
"Did not allocate requested number of vectors." << endl;
243 if ( MVT::GetGlobalLength(*B) != MVT::GetGlobalLength(*A) ) {
245 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
246 <<
"Did not allocate requested number of vectors." << endl;
249 MVT::MvNorm(*B, norms);
250 if ( norms.size() != 2*numvecs && ResizeWarning==false ) {
252 <<
"*** WARNING *** MultiVecTraits::MvNorm()." << endl
253 <<
"Method resized the output vector." << endl;
254 ResizeWarning =
true;
256 for (
int i=0; i<numvecs; i++) {
257 if ( norms[i] < zero_mag ) {
259 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
260 <<
"Vector had negative norm." << endl;
283 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
284 std::vector<MagType> norms(numvecs), norms2(numvecs);
287 MVT::MvNorm(*B, norms);
288 for (
int i=0; i<numvecs; i++) {
289 if ( norms[i] != zero_mag ) {
291 <<
"*** ERROR *** MultiVecTraits::MvInit() " 292 <<
"and MultiVecTraits::MvNorm()" << endl
293 <<
"Supposedly zero vector has non-zero norm." << endl;
298 MVT::MvNorm(*B, norms);
300 MVT::MvNorm(*B, norms2);
301 for (
int i=0; i<numvecs; i++) {
302 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) {
304 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
305 <<
"Random vector was empty (very unlikely)." << endl;
308 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) {
310 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
311 <<
"Vector had negative norm." << endl;
314 else if ( norms[i] == norms2[i] ) {
316 <<
"*** ERROR *** MutliVecTraits::MvRandom()." << endl
317 <<
"Vectors not random enough." << endl;
332 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
333 std::vector<MagType> norms(numvecs);
336 MVT::MvScale(*B,STS::zero());
337 MVT::MvNorm(*B, norms);
338 for (
int i=0; i<numvecs; i++) {
339 if ( norms[i] != zero_mag ) {
341 <<
"*** ERROR *** MultiVecTraits::MvScale(alpha) " 342 <<
"Supposedly zero vector has non-zero norm." << endl;
348 std::vector<ScalarType> zeros(numvecs,STS::zero());
349 MVT::MvScale(*B,zeros);
350 MVT::MvNorm(*B, norms);
351 for (
int i=0; i<numvecs; i++) {
352 if ( norms[i] != zero_mag ) {
354 <<
"*** ERROR *** MultiVecTraits::MvScale(alphas) " 355 <<
"Supposedly zero vector has non-zero norm." << endl;
376 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
377 std::vector<MagType> norms(numvecs);
380 MVT::MvNorm(*B, norms);
381 bool BadNormWarning =
false;
382 for (
int i=0; i<numvecs; i++) {
383 if ( norms[i] < zero_mag ) {
385 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
386 <<
"Vector had negative norm." << endl;
389 else if ( norms[i] != STS::squareroot(MVT::GetGlobalLength(*B)) && !BadNormWarning ) {
392 <<
"Warning testing MultiVecTraits::MvInit()." << endl
393 <<
"Ones std::vector should have norm std::sqrt(dim)." << endl
394 <<
"norms[i]: " << norms[i] <<
"\tdim: " << MVT::GetGlobalLength(*B) << endl << endl;
395 BadNormWarning =
true;
409 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
410 std::vector<MagType> norms(numvecs);
411 MVT::MvInit(*B, zero_mag);
412 MVT::MvNorm(*B, norms);
413 for (
int i=0; i<numvecs; i++) {
414 if ( norms[i] < zero_mag ) {
416 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
417 <<
"Vector had negative norm." << endl;
420 else if ( norms[i] != zero_mag ) {
422 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
423 <<
"Zero std::vector should have norm zero." << endl;
436 Teuchos::RCP<MV> B, C;
437 std::vector<MagType> norms(numvecs), norms2(numvecs);
439 B = MVT::Clone(*A,numvecs);
441 MVT::MvNorm(*B, norms);
442 C = MVT::CloneCopy(*B,ind);
443 MVT::MvNorm(*C, norms2);
444 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
446 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
447 <<
"Wrong number of vectors." << endl;
450 if ( MVT::GetGlobalLength(*C) != MVT::GetGlobalLength(*B) ) {
452 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
453 <<
"Vector lengths don't match." << endl;
456 for (
int i=0; i<numvecs_2; i++) {
457 if (STS::magnitude (norms2[i] - norms[ind[i]]) >
tol) {
459 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
460 <<
"Copied vectors do not agree: " 461 << norms2[i] <<
" != " << norms[ind[i]] << endl
462 <<
"Difference " << STS::magnitude (norms2[i] - norms[ind[i]])
463 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
469 MVT::MvInit(*B,zero);
470 MVT::MvNorm(*C, norms);
471 for (
int i=0; i<numvecs_2; i++) {
473 if (STS::magnitude (norms2[i] - norms[i]) >
tol) {
475 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
476 <<
"Copied vectors were not independent." << endl
477 << norms2[i] <<
" != " << norms[i] << endl
478 <<
"Difference " << STS::magnitude (norms2[i] - norms[i])
479 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
491 Teuchos::RCP<MV> B, C;
492 std::vector<MagType> norms(numvecs), norms2(numvecs);
494 B = MVT::Clone(*A,numvecs);
496 MVT::MvNorm(*B, norms);
497 C = MVT::CloneCopy(*B);
498 MVT::MvNorm(*C, norms2);
499 if ( MVT::GetNumberVecs(*C) != numvecs ) {
501 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
502 <<
"Wrong number of vectors." << endl;
505 for (
int i=0; i<numvecs; i++) {
506 if (STS::magnitude (norms2[i] - norms[i]) >
tol) {
508 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
509 <<
"Copied vectors do not agree: " 510 << norms2[i] <<
" != " << norms[i] << endl
511 <<
"Difference " << STS::magnitude (norms2[i] - norms[i])
512 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
516 MVT::MvInit(*B,zero);
517 MVT::MvNorm(*C, norms);
518 for (
int i=0; i<numvecs; i++) {
520 if (STS::magnitude (norms2[i] - norms[i]) >
tol) {
522 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
523 <<
"Copied vectors were not independent." << endl
524 << norms2[i] <<
" != " << norms[i] << endl
525 <<
"Difference " << STS::magnitude (norms2[i] - norms[i])
526 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
540 Teuchos::RCP<MV> B, C;
541 std::vector<MagType> norms(numvecs), norms2(numvecs);
543 B = MVT::Clone(*A,numvecs);
545 MVT::MvNorm(*B, norms);
546 C = MVT::CloneViewNonConst(*B,ind);
547 MVT::MvNorm(*C, norms2);
548 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
550 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
551 <<
"Wrong number of vectors." << endl;
554 for (
int i=0; i<numvecs_2; i++) {
556 if (STS::magnitude (norms2[i] - norms[ind[i]]) >
tol) {
558 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
559 <<
"Viewed vectors do not agree." << endl;
574 Teuchos::RCP<const MV> C;
575 std::vector<MagType> normsB(numvecs), normsC(numvecs_2);
576 std::vector<int> allind(numvecs);
577 for (
int i=0; i<numvecs; i++) {
581 B = MVT::Clone(*A,numvecs);
583 MVT::MvNorm(*B, normsB);
584 C = MVT::CloneView(*B,ind);
585 MVT::MvNorm(*C, normsC);
586 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
588 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
589 <<
"Wrong number of vectors." << endl;
592 for (
int i=0; i<numvecs_2; i++) {
594 if (STS::magnitude (normsC[i] - normsB[ind[i]]) >
tol) {
596 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
597 <<
"Viewed vectors do not agree." << endl;
616 Teuchos::RCP<MV> B, C;
617 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
618 normsC1(numvecs_2), normsC2(numvecs_2);
620 B = MVT::Clone(*A,numvecs);
621 C = MVT::Clone(*A,numvecs_2);
623 ind.resize(numvecs_2);
624 for (
int i=0; i<numvecs_2; i++) {
630 MVT::MvNorm(*B,normsB1);
631 MVT::MvNorm(*C,normsC1);
632 MVT::SetBlock(*C,ind,*B);
633 MVT::MvNorm(*B,normsB2);
634 MVT::MvNorm(*C,normsC2);
637 for (
int i=0; i<numvecs_2; i++) {
639 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
641 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
642 <<
"Operation modified source vectors." << endl;
648 for (
int i=0; i<numvecs; i++) {
651 if (STS::magnitude (normsB2[i] - normsC1[i/2]) >
tol) {
653 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
654 <<
"Copied vectors do not agree: " << endl
655 << normsB2[i] <<
" != " << normsC1[i/2] << endl
656 <<
"Difference " << STS::magnitude (normsB2[i] - normsC1[i/2])
657 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
663 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
665 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
666 <<
"Incorrect vectors were modified." << endl
667 << normsB1[i] <<
" != " << normsB2[i] << endl
668 <<
"Difference " << STS::magnitude (normsB2[i] - normsB2[i])
669 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
674 MVT::MvInit(*C,zero);
675 MVT::MvNorm(*B,normsB1);
677 for (
int i=0; i<numvecs; i++) {
679 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
681 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
682 <<
"Copied vectors were not independent." << endl
683 << normsB1[i] <<
" != " << normsB2[i] << endl
684 <<
"Difference " << STS::magnitude (normsB1[i] - normsB2[i])
685 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
708 Teuchos::RCP<MV> B, C;
713 std::vector<MagType> normsB1(BSize), normsB2(BSize),
714 normsC1(CSize), normsC2(CSize);
716 B = MVT::Clone(*A,BSize);
717 C = MVT::Clone(*A,CSize);
720 for (
int i=0; i<setSize; i++) {
726 MVT::MvNorm(*B,normsB1);
727 MVT::MvNorm(*C,normsC1);
728 MVT::SetBlock(*C,ind,*B);
729 MVT::MvNorm(*B,normsB2);
730 MVT::MvNorm(*C,normsC2);
733 for (
int i=0; i<CSize; i++) {
735 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
737 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
738 <<
"Operation modified source vectors." << endl;
744 for (
int i=0; i<BSize; i++) {
747 const MagType diff = STS::magnitude (normsB2[i] - normsC1[i/2]);
750 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
751 <<
"Copied vectors do not agree: " << endl
752 << normsB2[i] <<
" != " << normsC1[i/2] << endl
753 <<
"Difference " << diff <<
" exceeds the tolerance 100*eps = " 760 const MagType diff = STS::magnitude (normsB1[i] - normsB2[i]);
764 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
765 <<
"Incorrect vectors were modified." << endl
766 << normsB1[i] <<
" != " << normsB2[i] << endl
767 <<
"Difference " << diff <<
" exceeds the tolerance 100*eps = " 773 MVT::MvInit(*C,zero);
774 MVT::MvNorm(*B,normsB1);
776 for (
int i=0; i<numvecs; i++) {
778 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
780 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
781 <<
"Copied vectors were not independent." << endl
782 << normsB1[i] <<
" != " << normsB2[i] << endl
783 <<
"Difference " << STS::magnitude (normsB1[i] - normsB2[i])
784 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
813 Teuchos::RCP<MV> B, C;
814 std::vector<MagType> normsB(p), normsC(q);
815 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
817 B = MVT::Clone(*A,p);
818 C = MVT::Clone(*A,q);
822 MVT::MvNorm(*B,normsB);
824 MVT::MvNorm(*C,normsC);
827 MVT::MvTransMv( zero, *B, *C, SDM );
830 if ( SDM.numRows() != p || SDM.numCols() != q ) {
832 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
833 <<
"Routine resized SerialDenseMatrix." << endl;
838 if ( SDM.normOne() != zero ) {
840 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
841 <<
"Scalar argument processed incorrectly." << endl;
846 MVT::MvTransMv( one, *B, *C, SDM );
850 for (
int i=0; i<p; i++) {
851 for (
int j=0; j<q; j++) {
852 if ( STS::magnitude(SDM(i,j))
853 > STS::magnitude(normsB[i]*normsC[j]) ) {
855 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
856 <<
"Triangle inequality did not hold: " 857 << STS::magnitude(SDM(i,j))
859 << STS::magnitude(normsB[i]*normsC[j])
867 MVT::MvTransMv( one, *B, *C, SDM );
868 for (
int i=0; i<p; i++) {
869 for (
int j=0; j<q; j++) {
870 if ( SDM(i,j) != zero ) {
872 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
873 <<
"Inner products not zero for C==0." << endl;
880 MVT::MvTransMv( one, *B, *C, SDM );
881 for (
int i=0; i<p; i++) {
882 for (
int j=0; j<q; j++) {
883 if ( SDM(i,j) != zero ) {
885 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
886 <<
"Inner products not zero for B==0." << endl;
903 Teuchos::RCP<MV> B, C;
904 std::vector<ScalarType> iprods(p+q);
905 std::vector<MagType> normsB(numvecs), normsC(numvecs);
907 B = MVT::Clone(*A,p);
908 C = MVT::Clone(*A,p);
912 MVT::MvNorm(*B,normsB);
913 MVT::MvNorm(*C,normsC);
914 MVT::MvDot( *B, *C, iprods );
915 if ( iprods.size() != p+q ) {
917 <<
"*** ERROR *** MultiVecTraits::MvDot." << endl
918 <<
"Routine resized results std::vector." << endl;
922 if ( STS::magnitude(iprods[i])
923 > STS::magnitude(normsB[i]*normsC[i]) ) {
925 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
926 <<
"Inner products not valid." << endl;
932 MVT::MvDot( *B, *C, iprods );
933 for (
int i=0; i<p; i++) {
934 if ( iprods[i] != zero ) {
936 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
937 <<
"Inner products not zero for B==0." << endl;
943 MVT::MvDot( *B, *C, iprods );
944 for (
int i=0; i<p; i++) {
945 if ( iprods[i] != zero ) {
947 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
948 <<
"Inner products not zero for C==0." << endl;
964 Teuchos::RCP<MV> B, C, D;
965 std::vector<MagType> normsB1(p), normsB2(p),
966 normsC1(p), normsC2(p),
967 normsD1(p), normsD2(p);
968 ScalarType alpha = STS::random(),
969 beta = STS::random();
971 B = MVT::Clone(*A,p);
972 C = MVT::Clone(*A,p);
973 D = MVT::Clone(*A,p);
977 MVT::MvNorm(*B,normsB1);
978 MVT::MvNorm(*C,normsC1);
981 MVT::MvAddMv(zero,*B,one,*C,*D);
982 MVT::MvNorm(*B,normsB2);
983 MVT::MvNorm(*C,normsC2);
984 MVT::MvNorm(*D,normsD1);
985 for (
int i=0; i<p; i++) {
986 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
988 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
989 <<
"Input arguments were modified." << endl;
992 else if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
994 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
995 <<
"Input arguments were modified." << endl;
998 else if (STS::magnitude (normsC1[i] - normsD1[i]) >
tol) {
1000 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1001 <<
"Assignment did not work." << endl;
1007 MVT::MvAddMv(one,*B,zero,*C,*D);
1008 MVT::MvNorm(*B,normsB2);
1009 MVT::MvNorm(*C,normsC2);
1010 MVT::MvNorm(*D,normsD1);
1011 for (
int i=0; i<p; i++) {
1012 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1014 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1015 <<
"Input arguments were modified." << endl;
1018 else if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1020 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1021 <<
"Input arguments were modified." << endl;
1024 else if (STS::magnitude (normsB1[i] - normsD1[i]) >
tol) {
1026 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1027 <<
"Assignment did not work." << endl;
1035 MVT::MvAddMv(alpha,*B,beta,*C,*D);
1036 MVT::MvNorm(*B,normsB2);
1037 MVT::MvNorm(*C,normsC2);
1038 MVT::MvNorm(*D,normsD1);
1040 for (
int i=0; i<p; i++) {
1041 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1043 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1044 <<
"Input arguments were modified." << endl;
1047 else if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1049 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1050 <<
"Input arguments were modified." << endl;
1056 MVT::MvAddMv(alpha,*B,beta,*C,*D);
1057 MVT::MvNorm(*B,normsB2);
1058 MVT::MvNorm(*C,normsC2);
1059 MVT::MvNorm(*D,normsD2);
1062 for (
int i=0; i<p; i++) {
1063 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1065 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1066 <<
"Input arguments were modified." << endl;
1069 else if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1071 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1072 <<
"Input arguments were modified." << endl;
1075 else if (STS::magnitude (normsD1[i] - normsD2[i]) >
tol) {
1077 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
1078 <<
"Results varies depending on initial state of dest vectors." << endl;
1104 Teuchos::RCP<MV> B, D;
1105 Teuchos::RCP<const MV> C;
1106 std::vector<MagType> normsB(p),
1108 std::vector<int> lclindex(p);
1109 for (
int i=0; i<p; i++) lclindex[i] = i;
1111 B = MVT::Clone(*A,p);
1112 C = MVT::CloneView(*B,lclindex);
1113 D = MVT::CloneViewNonConst(*B,lclindex);
1116 MVT::MvNorm(*B,normsB);
1119 MVT::MvAddMv(zero,*B,one,*C,*D);
1120 MVT::MvNorm(*D,normsD);
1121 for (
int i=0; i<p; i++) {
1122 if (STS::magnitude (normsB[i] - normsD[i]) >
tol) {
1124 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1125 <<
"Assignment did not work." << endl;
1131 MVT::MvAddMv(one,*B,zero,*C,*D);
1132 MVT::MvNorm(*D,normsD);
1133 for (
int i=0; i<p; i++) {
1134 if (STS::magnitude (normsB[i] - normsD[i]) >
tol) {
1136 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1137 <<
"Assignment did not work." << endl;
1155 const int p = 7, q = 5;
1156 Teuchos::RCP<MV> B, C;
1157 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1158 std::vector<MagType> normsC1(q), normsC2(q),
1159 normsB1(p), normsB2(p);
1161 B = MVT::Clone(*A,p);
1162 C = MVT::Clone(*A,q);
1167 MVT::MvNorm(*B,normsB1);
1168 MVT::MvNorm(*C,normsC1);
1170 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1171 MVT::MvNorm(*B,normsB2);
1172 MVT::MvNorm(*C,normsC2);
1173 for (
int i=0; i<p; i++) {
1174 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1176 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1177 <<
"Input vectors were modified." << endl;
1181 for (
int i=0; i<q; i++) {
1182 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1184 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1185 <<
"Arithmetic test 1 failed." << endl;
1193 MVT::MvNorm(*B,normsB1);
1194 MVT::MvNorm(*C,normsC1);
1196 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1197 MVT::MvNorm(*B,normsB2);
1198 MVT::MvNorm(*C,normsC2);
1199 for (
int i=0; i<p; i++) {
1200 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1202 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1203 <<
"Input vectors were modified." << endl;
1207 for (
int i=0; i<q; i++) {
1208 if ( normsC2[i] != zero ) {
1210 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1211 <<
"Arithmetic test 2 failed: " 1224 MVT::MvNorm(*B,normsB1);
1225 MVT::MvNorm(*C,normsC1);
1227 for (
int i=0; i<q; i++) {
1230 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1231 MVT::MvNorm(*B,normsB2);
1232 MVT::MvNorm(*C,normsC2);
1233 for (
int i=0; i<p; i++) {
1234 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1236 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1237 <<
"Input vectors were modified." << endl;
1241 for (
int i=0; i<q; i++) {
1242 if (STS::magnitude (normsB1[i] - normsC2[i]) >
tol) {
1244 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1245 <<
"Arithmetic test 3 failed: " 1257 MVT::MvNorm(*B,normsB1);
1258 MVT::MvNorm(*C,normsC1);
1260 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1261 MVT::MvNorm(*B,normsB2);
1262 MVT::MvNorm(*C,normsC2);
1263 for (
int i=0; i<p; i++) {
1264 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1266 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1267 <<
"Input vectors were modified." << endl;
1271 for (
int i=0; i<q; i++) {
1272 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1274 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1275 <<
"Arithmetic test 4 failed." << endl;
1291 const int p = 5, q = 7;
1292 Teuchos::RCP<MV> B, C;
1293 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1294 std::vector<MagType> normsC1(q), normsC2(q),
1295 normsB1(p), normsB2(p);
1297 B = MVT::Clone(*A,p);
1298 C = MVT::Clone(*A,q);
1303 MVT::MvNorm(*B,normsB1);
1304 MVT::MvNorm(*C,normsC1);
1306 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1307 MVT::MvNorm(*B,normsB2);
1308 MVT::MvNorm(*C,normsC2);
1309 for (
int i=0; i<p; i++) {
1310 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1312 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1313 <<
"Input vectors were modified." << endl;
1317 for (
int i=0; i<q; i++) {
1318 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1320 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1321 <<
"Arithmetic test 5 failed." << endl;
1329 MVT::MvNorm(*B,normsB1);
1330 MVT::MvNorm(*C,normsC1);
1332 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1333 MVT::MvNorm(*B,normsB2);
1334 MVT::MvNorm(*C,normsC2);
1335 for (
int i=0; i<p; i++) {
1336 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1338 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1339 <<
"Input vectors were modified." << endl;
1343 for (
int i=0; i<q; i++) {
1344 if ( normsC2[i] != zero ) {
1346 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1347 <<
"Arithmetic test 6 failed: " 1359 MVT::MvNorm(*B,normsB1);
1360 MVT::MvNorm(*C,normsC1);
1362 for (
int i=0; i<p; i++) {
1365 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1366 MVT::MvNorm(*B,normsB2);
1367 MVT::MvNorm(*C,normsC2);
1368 for (
int i=0; i<p; i++) {
1369 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1371 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1372 <<
"Input vectors were modified." << endl;
1376 for (
int i=0; i<p; i++) {
1377 if (STS::magnitude (normsB1[i] - normsC2[i]) >
tol) {
1379 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1380 <<
"Arithmetic test 7 failed." << endl;
1384 for (
int i=p; i<q; i++) {
1385 if ( normsC2[i] != zero ) {
1387 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1388 <<
"Arithmetic test 7 failed." << endl;
1396 MVT::MvNorm(*B,normsB1);
1397 MVT::MvNorm(*C,normsC1);
1399 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1400 MVT::MvNorm(*B,normsB2);
1401 MVT::MvNorm(*C,normsC2);
1402 for (
int i=0; i<p; i++) {
1403 if (STS::magnitude (normsB1[i] - normsB2[i]) >
tol) {
1405 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1406 <<
"Input vectors were modified." << endl;
1410 for (
int i=0; i<q; i++) {
1411 if (STS::magnitude (normsC1[i] - normsC2[i]) >
tol) {
1413 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1414 <<
"Arithmetic test 8 failed." << endl;
1447 template<
class ScalarType,
class MV,
class OP>
1450 const Teuchos::RCP<const MV> &A,
1451 const Teuchos::RCP<const OP> &M)
1453 using Teuchos::SetScientific;
1456 typedef Teuchos::ScalarTraits<ScalarType> STS;
1457 typedef typename STS::magnitudeType MagType;
1461 SetScientific<ScalarType> sci (om->stream (
Warnings));
1466 const MagType
tol = Teuchos::as<MagType> (100) * STS::eps ();
1476 typedef Teuchos::ScalarTraits<ScalarType> STS;
1478 typedef typename STS::magnitudeType MagType;
1480 const int numvecs = 10;
1482 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
1483 C = MVT::Clone(*A,numvecs);
1485 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
1486 normsC1(numvecs), normsC2(numvecs);
1487 bool NonDeterministicWarning;
1498 MVT::MvNorm(*B,normsB1);
1499 OPT::Apply(*M,*B,*C);
1500 MVT::MvNorm(*B,normsB2);
1501 MVT::MvNorm(*C,normsC2);
1502 for (
int i=0; i<numvecs; i++) {
1503 if (STS::magnitude (normsB2[i] - normsB1[i]) >
tol) {
1505 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1506 <<
"Apply() modified the input vectors." << endl
1507 <<
"Original: " << normsB1[i] <<
"; After: " << normsB2[i] << endl
1508 <<
"Difference " << STS::magnitude (normsB2[i] - normsB1[i])
1509 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
1512 if (normsC2[i] != STS::zero()) {
1514 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1515 <<
"Operator applied to zero did not return zero." << endl;
1522 MVT::MvNorm(*B,normsB1);
1523 OPT::Apply(*M,*B,*C);
1524 MVT::MvNorm(*B,normsB2);
1525 MVT::MvNorm(*C,normsC2);
1526 bool ZeroWarning =
false;
1527 for (
int i=0; i<numvecs; i++) {
1528 if (STS::magnitude (normsB2[i] - normsB1[i]) >
tol) {
1530 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1531 <<
"Apply() modified the input vectors." << endl
1532 <<
"Original: " << normsB1[i] <<
"; After: " << normsB2[i] << endl
1533 <<
"Difference " << STS::magnitude (normsB2[i] - normsB1[i])
1534 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
1537 if (normsC2[i] == STS::zero() && ZeroWarning==
false ) {
1539 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1540 <<
"Operator applied to random vectors returned zero." << endl;
1547 MVT::MvNorm(*B,normsB1);
1549 OPT::Apply(*M,*B,*C);
1550 MVT::MvNorm(*B,normsB2);
1551 MVT::MvNorm(*C,normsC1);
1552 for (
int i=0; i<numvecs; i++) {
1553 if (STS::magnitude (normsB2[i] - normsB1[i]) >
tol) {
1555 <<
"*** ERROR *** OperatorTraits::Apply() [3]" << endl
1556 <<
"Apply() modified the input vectors." << endl
1557 <<
"Original: " << normsB1[i] <<
"; After: " << normsB2[i] << endl
1558 <<
"Difference " << STS::magnitude (normsB2[i] - normsB1[i])
1559 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
1570 OPT::Apply(*M,*B,*C);
1571 MVT::MvNorm(*B,normsB2);
1572 MVT::MvNorm(*C,normsC2);
1573 NonDeterministicWarning =
false;
1574 for (
int i=0; i<numvecs; i++) {
1575 if (STS::magnitude (normsB2[i] - normsB1[i]) >
tol) {
1577 <<
"*** ERROR *** OperatorTraits::Apply() [4]" << endl
1578 <<
"Apply() modified the input vectors." << endl
1579 <<
"Original: " << normsB1[i] <<
"; After: " << normsB2[i] << endl
1580 <<
"Difference " << STS::magnitude (normsB2[i] - normsB1[i])
1581 <<
" exceeds the tolerance 100*eps = " <<
tol << endl;
1584 if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) {
1587 <<
"*** WARNING *** OperatorTraits::Apply() [4]" << endl
1588 <<
"Apply() returned two different results." << endl << endl;
1589 NonDeterministicWarning =
true;
Collection of types and exceptions used within the Belos solvers.
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.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
bool TestOperatorTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A, const Teuchos::RCP< const OP > &M)
Test correctness of OperatorTraits specialization and its operator implementation.
Class which defines basic traits for the operator type.
bool TestMultiVecTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A)
Test correctness of a MultiVecTraits specialization and multivector implementation.
Belos header file which uses auto-configuration information to include necessary C++ headers...