29 #ifndef ANASAZI_MVOPTESTER_HPP 30 #define ANASAZI_MVOPTESTER_HPP 53 #include "Teuchos_MatrixMarket_SetScientific.hpp" 54 #include "Teuchos_RCP.hpp" 55 #include "Teuchos_as.hpp" 64 template<
class ScalarType,
class MV >
67 const Teuchos::RCP<const MV> &A ) {
70 using Teuchos::MatrixMarket::details::SetScientific;
147 typedef Teuchos::ScalarTraits<ScalarType> SCT;
148 typedef typename SCT::magnitudeType MagType;
150 const ScalarType one = SCT::one();
151 const ScalarType zero = SCT::zero();
152 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero();
153 const MagType tol = SCT::eps()*100;
156 const int numvecs = 10;
157 const int numvecs_2 = 5;
159 std::vector<int> ind(numvecs_2);
169 TEUCHOS_TEST_FOR_EXCEPT(numvecs_2 != 5);
180 if ( MVT::GetNumberVecs(*A) <= 0 ) {
182 <<
"*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl
183 <<
"Returned <= 0." << endl;
191 if ( MVT::GetGlobalLength(*A) <= 0 ) {
193 <<
"*** ERROR *** MultiVectorTraitsExt::GetGlobalLength()" << endl
194 <<
"Returned <= 0." << endl;
205 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
206 std::vector<MagType> norms(2*numvecs);
207 bool ResizeWarning =
false;
208 if ( MVT::GetNumberVecs(*B) != numvecs ) {
210 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
211 <<
"Did not allocate requested number of vectors." << endl;
214 if ( MVT::GetGlobalLength(*B) != MVT::GetGlobalLength(*A) ) {
216 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
217 <<
"Did not allocate requested number of vectors." << endl;
220 MVT::MvNorm(*B, norms);
221 if ( (
int)norms.size() != 2*numvecs && (ResizeWarning ==
false) ) {
223 <<
"*** WARNING *** MultiVecTraits::MvNorm()." << endl
224 <<
"Method resized the output vector." << endl;
225 ResizeWarning =
true;
227 for (
int i=0; i<numvecs; i++) {
228 if ( norms[i] < zero_mag ) {
230 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
231 <<
"Vector had negative norm." << endl;
254 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
255 std::vector<MagType> norms(numvecs), norms2(numvecs);
258 MVT::MvNorm(*B, norms);
259 for (
int i=0; i<numvecs; i++) {
260 if ( norms[i] != zero_mag ) {
262 <<
"*** ERROR *** MultiVecTraits::MvInit() " 263 <<
"and MultiVecTraits::MvNorm()" << endl
264 <<
"Supposedly zero vector has non-zero norm." << endl;
269 MVT::MvNorm(*B, norms);
271 MVT::MvNorm(*B, norms2);
272 for (
int i=0; i<numvecs; i++) {
273 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) {
275 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
276 <<
"Random vector was empty (very unlikely)." << endl;
279 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) {
281 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
282 <<
"Vector had negative norm." << endl;
285 else if ( norms[i] == norms2[i] ) {
287 <<
"*** ERROR *** MutliVecTraits::MvRandom()." << endl
288 <<
"Vectors not random enough." << endl;
303 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
304 std::vector<MagType> norms(numvecs);
307 MVT::MvScale(*B,SCT::zero());
308 MVT::MvNorm(*B, norms);
309 for (
int i=0; i<numvecs; i++) {
310 if ( norms[i] != zero_mag ) {
312 <<
"*** ERROR *** MultiVecTraits::MvScale(alpha) " 313 <<
"Supposedly zero vector has non-zero norm." << endl;
319 std::vector<ScalarType> zeros(numvecs,SCT::zero());
320 MVT::MvScale(*B,zeros);
321 MVT::MvNorm(*B, norms);
322 for (
int i=0; i<numvecs; i++) {
323 if ( norms[i] != zero_mag ) {
325 <<
"*** ERROR *** MultiVecTraits::MvScale(alphas) " 326 <<
"Supposedly zero vector has non-zero norm." << endl;
347 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
348 std::vector<MagType> norms(numvecs);
351 MVT::MvNorm(*B, norms);
352 bool BadNormWarning =
false;
353 for (
int i=0; i<numvecs; i++) {
354 if ( norms[i] < zero_mag ) {
356 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
357 <<
"Vector had negative norm." << endl;
360 else if ( norms[i] != SCT::squareroot(MVT::GetGlobalLength(*B)) && !BadNormWarning ) {
363 <<
"Warning testing MultiVecTraits::MvInit()." << endl
364 <<
"Ones vector should have norm sqrt(dim)." << endl
365 <<
"norms[i]: " << norms[i] <<
"\tdim: " << MVT::GetGlobalLength(*B) << endl << endl;
366 BadNormWarning =
true;
380 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
381 std::vector<MagType> norms(numvecs);
382 MVT::MvInit(*B, zero_mag);
383 MVT::MvNorm(*B, norms);
384 for (
int i=0; i<numvecs; i++) {
385 if ( norms[i] < zero_mag ) {
387 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
388 <<
"Vector had negative norm." << endl;
391 else if ( norms[i] != zero_mag ) {
393 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
394 <<
"Zero vector should have norm zero." << endl;
407 Teuchos::RCP<MV> B, C;
408 std::vector<MagType> norms(numvecs), norms2(ind.size());
410 B = MVT::Clone(*A,numvecs);
412 MVT::MvNorm(*B, norms);
413 C = MVT::CloneCopy(*B,ind);
414 MVT::MvNorm(*C, norms2);
415 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
417 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
418 <<
"Wrong number of vectors." << endl;
421 if ( MVT::GetGlobalLength(*C) != MVT::GetGlobalLength(*B) ) {
423 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
424 <<
"Vector lengths don't match." << endl;
427 for (
int i=0; i<numvecs_2; i++) {
428 if ( SCT::magnitude( norms2[i] - norms[ind[i]] ) > tol ) {
430 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
431 <<
"Copied vectors do not agree:" 432 << norms2[i] <<
" != " << norms[ind[i]] << endl
433 <<
"Difference " << SCT::magnitude (norms2[i] - norms[ind[i]])
434 <<
" exceeds the tolerance 100*eps = " << tol << endl;
439 MVT::MvInit(*B,zero);
440 MVT::MvNorm(*C, norms2);
441 for (
int i=0; i<numvecs_2; i++) {
442 if ( SCT::magnitude( norms2[i] - norms2[i] ) > tol ) {
444 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
445 <<
"Copied vectors were not independent." << endl
446 <<
"Difference " << SCT::magnitude (norms2[i] - norms[i])
447 <<
" exceeds the tolerance 100*eps = " << tol << endl;
460 Teuchos::RCP<MV> B, C;
461 std::vector<MagType> norms(numvecs), norms2(numvecs);
463 B = MVT::Clone(*A,numvecs);
465 MVT::MvNorm(*B, norms);
466 C = MVT::CloneCopy(*B);
467 MVT::MvNorm(*C, norms2);
468 if ( MVT::GetNumberVecs(*C) != numvecs ) {
470 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
471 <<
"Wrong number of vectors." << endl;
474 for (
int i=0; i<numvecs; i++) {
475 if ( SCT::magnitude( norms2[i] - norms[i] ) > tol ) {
477 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
478 <<
"Copied vectors do not agree." << endl
479 <<
"Difference " << SCT::magnitude (norms2[i] - norms[i])
480 <<
" exceeds the tolerance 100*eps = " << tol << endl;
484 MVT::MvInit(*B,zero);
485 MVT::MvNorm(*C, norms);
486 for (
int i=0; i<numvecs; i++) {
487 if ( SCT::magnitude( norms2[i] - norms[i] ) > tol ) {
489 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
490 <<
"Copied vectors were not independent." << endl
491 <<
"Difference " << SCT::magnitude (norms2[i] - norms[i])
492 <<
" exceeds the tolerance 100*eps = " << tol << endl;
506 Teuchos::RCP<MV> B, C;
507 std::vector<MagType> norms(numvecs), norms2(ind.size());
509 B = MVT::Clone(*A,numvecs);
511 MVT::MvNorm(*B, norms);
512 C = MVT::CloneViewNonConst(*B,ind);
513 MVT::MvNorm(*C, norms2);
514 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
516 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
517 <<
"Wrong number of vectors." << endl;
520 for (
int i=0; i<numvecs_2; i++) {
521 if ( SCT::magnitude( norms2[i] - norms[ind[i]] ) > tol ) {
523 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
524 <<
"Viewed vectors do not agree." << endl;
539 Teuchos::RCP<const MV> constB, C;
540 std::vector<MagType> normsB(numvecs), normsC(ind.size());
541 std::vector<int> allind(numvecs);
542 for (
int i=0; i<numvecs; i++) {
546 B = MVT::Clone(*A,numvecs);
548 MVT::MvNorm(*B, normsB);
549 C = MVT::CloneView(*B,ind);
550 MVT::MvNorm(*C, normsC);
551 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
553 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
554 <<
"Wrong number of vectors." << endl;
557 for (
int i=0; i<numvecs_2; i++) {
558 if ( SCT::magnitude( normsC[i] - normsB[ind[i]] ) > tol ) {
560 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
561 <<
"Viewed vectors do not agree." << endl;
580 Teuchos::RCP<MV> B, C;
581 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
582 normsC1(numvecs_2), normsC2(numvecs_2);
584 B = MVT::Clone(*A,numvecs);
585 C = MVT::Clone(*A,numvecs_2);
587 ind.resize(numvecs_2);
588 for (
int i=0; i<numvecs_2; i++) {
594 MVT::MvNorm(*B,normsB1);
595 MVT::MvNorm(*C,normsC1);
596 MVT::SetBlock(*C,ind,*B);
597 MVT::MvNorm(*B,normsB2);
598 MVT::MvNorm(*C,normsC2);
601 for (
int i=0; i<numvecs_2; i++) {
602 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
604 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
605 <<
"Operation modified source vectors." << endl;
611 for (
int i=0; i<numvecs; i++) {
614 if ( SCT::magnitude(normsB2[i]-normsC1[i/2]) > tol ) {
616 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
617 <<
"Copied vectors do not agree." << endl
618 <<
"Difference " << SCT::magnitude (normsB2[i] - normsC1[i/2])
619 <<
" exceeds the tolerance 100*eps = " << tol << endl;
625 if ( SCT::magnitude(normsB1[i]-normsB2[i]) > tol ) {
627 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
628 <<
"Incorrect vectors were modified." << endl;
633 MVT::MvInit(*C,zero);
634 MVT::MvNorm(*B,normsB1);
636 for (
int i=0; i<numvecs; i++) {
637 if ( SCT::magnitude(normsB1[i]-normsB2[i]) > tol ) {
639 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
640 <<
"Copied vectors were not independent." << endl;
673 Teuchos::RCP<MV> B, C;
674 std::vector<MagType> normsB(p), normsC(q);
675 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
677 B = MVT::Clone(*A,p);
678 C = MVT::Clone(*A,q);
682 MVT::MvNorm(*B,normsB);
684 MVT::MvNorm(*C,normsC);
687 MVT::MvTransMv( zero, *B, *C, SDM );
690 if ( SDM.numRows() != p || SDM.numCols() != q ) {
692 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
693 <<
"Routine resized SerialDenseMatrix." << endl;
698 if ( SDM.normOne() != zero ) {
700 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
701 <<
"Scalar argument processed incorrectly." << endl;
706 MVT::MvTransMv( one, *B, *C, SDM );
710 for (
int i=0; i<p; i++) {
711 for (
int j=0; j<q; j++) {
712 if ( SCT::magnitude(SDM(i,j))
713 > SCT::magnitude(normsB[i]*normsC[j]) ) {
715 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
716 <<
"Triangle inequality did not hold: " 717 << SCT::magnitude(SDM(i,j))
719 << SCT::magnitude(normsB[i]*normsC[j])
727 MVT::MvTransMv( one, *B, *C, SDM );
728 for (
int i=0; i<p; i++) {
729 for (
int j=0; j<q; j++) {
730 if ( SDM(i,j) != zero ) {
732 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
733 <<
"Inner products not zero for C==0." << endl;
740 MVT::MvTransMv( one, *B, *C, SDM );
741 for (
int i=0; i<p; i++) {
742 for (
int j=0; j<q; j++) {
743 if ( SDM(i,j) != zero ) {
745 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
746 <<
"Inner products not zero for B==0." << endl;
759 Teuchos::SerialDenseMatrix<int, ScalarType> largeSDM(p+1,q+1);
760 Teuchos::SerialDenseMatrix<int, ScalarType> SDM2(Teuchos::View, largeSDM, p, q);
761 largeSDM.putScalar( one );
764 MVT::MvTransMv( one, *B, *C, SDM2 );
765 for (
int i=0; i<p; i++) {
766 for (
int j=0; j<q; j++) {
767 if ( SDM2(i,j) != zero ) {
769 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
770 <<
"Inner products not zero for C==0 when using a view into Teuchos::SerialDenseMatrix<>." << endl;
787 Teuchos::RCP<MV> B, C;
788 std::vector<ScalarType> iprods(q);
789 std::vector<MagType> normsB(p), normsC(p);
791 B = MVT::Clone(*A,p);
792 C = MVT::Clone(*A,p);
796 MVT::MvNorm(*B,normsB);
797 MVT::MvNorm(*C,normsC);
798 MVT::MvDot( *B, *C, iprods );
799 if ( (
int)iprods.size() != q ) {
801 <<
"*** ERROR *** MultiVecTraits::MvDot." << endl
802 <<
"Routine resized results vector." << endl;
805 for (
int i=0; i<p; i++) {
806 if ( SCT::magnitude(iprods[i])
807 > SCT::magnitude(normsB[i]*normsC[i]) ) {
809 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
810 <<
"Inner products not valid." << endl;
816 MVT::MvDot( *B, *C, iprods );
817 for (
int i=0; i<p; i++) {
818 if ( iprods[i] != zero ) {
820 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
821 <<
"Inner products not zero for B==0." << endl;
827 MVT::MvDot( *B, *C, iprods );
828 for (
int i=0; i<p; i++) {
829 if ( iprods[i] != zero ) {
831 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
832 <<
"Inner products not zero for C==0." << endl;
848 Teuchos::RCP<MV> B, C, D;
849 std::vector<MagType> normsB1(p), normsB2(p),
850 normsC1(p), normsC2(p),
851 normsD1(p), normsD2(p);
852 ScalarType alpha = SCT::random(),
853 beta = SCT::random();
855 B = MVT::Clone(*A,p);
856 C = MVT::Clone(*A,p);
857 D = MVT::Clone(*A,p);
861 MVT::MvNorm(*B,normsB1);
862 MVT::MvNorm(*C,normsC1);
865 MVT::MvAddMv(zero,*B,one,*C,*D);
866 MVT::MvNorm(*B,normsB2);
867 MVT::MvNorm(*C,normsC2);
868 MVT::MvNorm(*D,normsD1);
869 for (
int i=0; i<p; i++) {
870 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
872 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
873 <<
"Input arguments were modified." << endl;
876 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
878 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
879 <<
"Input arguments were modified." << endl;
882 else if ( SCT::magnitude(normsC1[i]-normsD1[i]) > tol ) {
884 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
885 <<
"Assignment did not work." << endl;
891 MVT::MvAddMv(one,*B,zero,*C,*D);
892 MVT::MvNorm(*B,normsB2);
893 MVT::MvNorm(*C,normsC2);
894 MVT::MvNorm(*D,normsD1);
895 for (
int i=0; i<p; i++) {
896 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
898 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
899 <<
"Input arguments were modified." << endl;
902 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
904 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
905 <<
"Input arguments were modified." << endl;
908 else if ( SCT::magnitude( normsB1[i] - normsD1[i] ) > tol ) {
910 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
911 <<
"Assignment did not work." << endl;
919 MVT::MvAddMv(alpha,*B,beta,*C,*D);
920 MVT::MvNorm(*B,normsB2);
921 MVT::MvNorm(*C,normsC2);
922 MVT::MvNorm(*D,normsD1);
924 for (
int i=0; i<p; i++) {
925 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
927 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
928 <<
"Input arguments were modified." << endl;
931 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
933 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
934 <<
"Input arguments were modified." << endl;
940 MVT::MvAddMv(alpha,*B,beta,*C,*D);
941 MVT::MvNorm(*B,normsB2);
942 MVT::MvNorm(*C,normsC2);
943 MVT::MvNorm(*D,normsD2);
946 for (
int i=0; i<p; i++) {
947 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
949 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
950 <<
"Input arguments were modified." << endl;
953 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
955 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
956 <<
"Input arguments were modified." << endl;
959 else if ( SCT::magnitude( normsD1[i] - normsD2[i] ) > tol ) {
961 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
962 <<
"Results varies depending on initial state of dest vectors." << endl;
988 Teuchos::RCP<MV> B, D;
989 Teuchos::RCP<const MV> C;
990 std::vector<MagType> normsB(p),
992 std::vector<int> lclindex(p);
993 for (
int i=0; i<p; i++) lclindex[i] = i;
995 B = MVT::Clone(*A,p);
996 C = MVT::CloneView(*B,lclindex);
997 D = MVT::CloneViewNonConst(*B,lclindex);
1000 MVT::MvNorm(*B,normsB);
1003 MVT::MvAddMv(zero,*B,one,*C,*D);
1004 MVT::MvNorm(*D,normsD);
1005 for (
int i=0; i<p; i++) {
1006 if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
1008 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1009 <<
"Assignment did not work." << endl;
1015 MVT::MvAddMv(one,*B,zero,*C,*D);
1016 MVT::MvNorm(*D,normsD);
1017 for (
int i=0; i<p; i++) {
1018 if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
1020 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1021 <<
"Assignment did not work." << endl;
1039 const int p = 7, q = 5;
1040 Teuchos::RCP<MV> B, C;
1041 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1042 std::vector<MagType> normsC1(q), normsC2(q),
1043 normsB1(p), normsB2(p);
1045 B = MVT::Clone(*A,p);
1046 C = MVT::Clone(*A,q);
1051 MVT::MvNorm(*B,normsB1);
1052 MVT::MvNorm(*C,normsC1);
1054 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1055 MVT::MvNorm(*B,normsB2);
1056 MVT::MvNorm(*C,normsC2);
1057 for (
int i=0; i<p; i++) {
1058 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1060 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1061 <<
"Input vectors were modified." << endl;
1065 for (
int i=0; i<q; i++) {
1066 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1068 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1069 <<
"Arithmetic test 1 failed." << endl;
1077 MVT::MvNorm(*B,normsB1);
1078 MVT::MvNorm(*C,normsC1);
1080 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1081 MVT::MvNorm(*B,normsB2);
1082 MVT::MvNorm(*C,normsC2);
1083 for (
int i=0; i<p; i++) {
1084 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1086 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1087 <<
"Input vectors were modified." << endl;
1091 for (
int i=0; i<q; i++) {
1092 if ( normsC2[i] != zero ) {
1094 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1095 <<
"Arithmetic test 2 failed: " 1108 MVT::MvNorm(*B,normsB1);
1109 MVT::MvNorm(*C,normsC1);
1111 for (
int i=0; i<q; i++) {
1114 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1115 MVT::MvNorm(*B,normsB2);
1116 MVT::MvNorm(*C,normsC2);
1117 for (
int i=0; i<p; i++) {
1118 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1120 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1121 <<
"Input vectors were modified." << endl;
1125 for (
int i=0; i<q; i++) {
1126 if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1128 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1129 <<
"Arithmetic test 3 failed: " 1141 MVT::MvNorm(*B,normsB1);
1142 MVT::MvNorm(*C,normsC1);
1144 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1145 MVT::MvNorm(*B,normsB2);
1146 MVT::MvNorm(*C,normsC2);
1147 for (
int i=0; i<p; i++) {
1148 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1150 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1151 <<
"Input vectors were modified." << endl;
1155 for (
int i=0; i<q; i++) {
1156 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1158 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1159 <<
"Arithmetic test 4 failed." << endl;
1175 const int p = 5, q = 7;
1176 Teuchos::RCP<MV> B, C;
1177 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1178 std::vector<MagType> normsC1(q), normsC2(q),
1179 normsB1(p), normsB2(p);
1181 B = MVT::Clone(*A,p);
1182 C = MVT::Clone(*A,q);
1187 MVT::MvNorm(*B,normsB1);
1188 MVT::MvNorm(*C,normsC1);
1190 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1191 MVT::MvNorm(*B,normsB2);
1192 MVT::MvNorm(*C,normsC2);
1193 for (
int i=0; i<p; i++) {
1194 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1196 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1197 <<
"Input vectors were modified." << endl;
1201 for (
int i=0; i<q; i++) {
1202 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1204 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1205 <<
"Arithmetic test 5 failed." << endl;
1213 MVT::MvNorm(*B,normsB1);
1214 MVT::MvNorm(*C,normsC1);
1216 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1217 MVT::MvNorm(*B,normsB2);
1218 MVT::MvNorm(*C,normsC2);
1219 for (
int i=0; i<p; i++) {
1220 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1222 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1223 <<
"Input vectors were modified." << endl;
1227 for (
int i=0; i<q; i++) {
1228 if ( normsC2[i] != zero ) {
1230 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1231 <<
"Arithmetic test 6 failed: " 1243 MVT::MvNorm(*B,normsB1);
1244 MVT::MvNorm(*C,normsC1);
1246 for (
int i=0; i<p; i++) {
1249 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1250 MVT::MvNorm(*B,normsB2);
1251 MVT::MvNorm(*C,normsC2);
1252 for (
int i=0; i<p; i++) {
1253 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1255 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1256 <<
"Input vectors were modified." << endl;
1260 for (
int i=0; i<p; i++) {
1261 if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1263 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1264 <<
"Arithmetic test 7 failed." << endl;
1268 for (
int i=p; i<q; i++) {
1269 if ( normsC2[i] != zero ) {
1271 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1272 <<
"Arithmetic test 7 failed." << endl;
1280 MVT::MvNorm(*B,normsB1);
1281 MVT::MvNorm(*C,normsC1);
1283 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1284 MVT::MvNorm(*B,normsB2);
1285 MVT::MvNorm(*C,normsC2);
1286 for (
int i=0; i<p; i++) {
1287 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1289 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1290 <<
"Input vectors were modified." << endl;
1294 for (
int i=0; i<q; i++) {
1295 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1297 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1298 <<
"Arithmetic test 8 failed." << endl;
1315 template<
class ScalarType,
class MV,
class OP>
1318 const Teuchos::RCP<const MV> &A,
1319 const Teuchos::RCP<const OP> &M) {
1331 typedef Teuchos::ScalarTraits<ScalarType> SCT;
1333 typedef typename SCT::magnitudeType MagType;
1335 const MagType tol = SCT::eps()*100;
1337 const int numvecs = 10;
1339 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
1340 C = MVT::Clone(*A,numvecs);
1342 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
1343 normsC1(numvecs), normsC2(numvecs);
1354 MVT::MvNorm(*B,normsB1);
1355 OPT::Apply(*M,*B,*C);
1356 MVT::MvNorm(*B,normsB2);
1357 MVT::MvNorm(*C,normsC2);
1358 for (
int i=0; i<numvecs; i++) {
1359 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1361 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1362 <<
"Apply() modified the input vectors." << endl;
1365 if (normsC2[i] != SCT::zero()) {
1367 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1368 <<
"Operator applied to zero did not return zero." << endl;
1375 MVT::MvNorm(*B,normsB1);
1376 OPT::Apply(*M,*B,*C);
1377 MVT::MvNorm(*B,normsB2);
1378 MVT::MvNorm(*C,normsC2);
1379 bool ZeroWarning =
false;
1380 for (
int i=0; i<numvecs; i++) {
1381 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1383 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1384 <<
"Apply() modified the input vectors." << endl;
1387 if (normsC2[i] == SCT::zero() && ZeroWarning==
false ) {
1389 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1390 <<
"Operator applied to random vectors returned zero." << endl;
1397 MVT::MvNorm(*B,normsB1);
1399 OPT::Apply(*M,*B,*C);
1400 MVT::MvNorm(*B,normsB2);
1401 MVT::MvNorm(*C,normsC1);
1402 for (
int i=0; i<numvecs; i++) {
1403 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1405 <<
"*** ERROR *** OperatorTraits::Apply() [3]" << endl
1406 <<
"Apply() modified the input vectors." << endl;
1417 OPT::Apply(*M,*B,*C);
1418 MVT::MvNorm(*B,normsB2);
1419 MVT::MvNorm(*C,normsC2);
1420 for (
int i=0; i<numvecs; i++) {
1421 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1423 <<
"*** ERROR *** OperatorTraits::Apply() [4]" << endl
1424 <<
"Apply() modified the input vectors." << endl;
1427 if ( SCT::magnitude( normsC1[i] - normsC2[i]) > tol ) {
1430 <<
"*** WARNING *** OperatorTraits::Apply() [4]" << endl
1431 <<
"Apply() returned two different results." << endl << endl;
bool TestOperatorTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A, const Teuchos::RCP< const OP > &M)
This function tests the correctness of an operator implementation with respect to an OperatorTraits s...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Abstract class definition for Anasazi Output Managers.
Output managers remove the need for the eigensolver to know any information about the required output...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
bool TestMultiVecTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A)
This is a function to test the correctness of a MultiVecTraits specialization and multivector impleme...
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Types and exceptions used within Anasazi solvers and interfaces.