29 #ifndef ANASAZI_MVOPTESTER_HPP 30 #define ANASAZI_MVOPTESTER_HPP 53 #include "Teuchos_RCP.hpp" 54 #include "Teuchos_as.hpp" 63 template<
class ScalarType,
class MV >
66 const Teuchos::RCP<const MV> &A ) {
145 typedef Teuchos::ScalarTraits<ScalarType> SCT;
146 typedef typename SCT::magnitudeType MagType;
148 const ScalarType one = SCT::one();
149 const ScalarType zero = SCT::zero();
150 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero();
153 const int numvecs = 10;
154 const int numvecs_2 = 5;
156 std::vector<int> ind(numvecs_2);
166 TEUCHOS_TEST_FOR_EXCEPT(numvecs_2 != 5);
177 if ( MVT::GetNumberVecs(*A) <= 0 ) {
179 <<
"*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl
180 <<
"Returned <= 0." << endl;
188 if ( MVT::GetGlobalLength(*A) <= 0 ) {
190 <<
"*** ERROR *** MultiVectorTraitsExt::GetGlobalLength()" << endl
191 <<
"Returned <= 0." << endl;
202 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
203 std::vector<MagType> norms(numvecs);
204 if ( MVT::GetNumberVecs(*B) != numvecs ) {
206 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
207 <<
"Did not allocate requested number of vectors." << endl;
210 if ( MVT::GetGlobalLength(*B) != MVT::GetGlobalLength(*A) ) {
212 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
213 <<
"Did not allocate requested number of vectors." << endl;
216 MVT::MvNorm(*B, norms);
217 if ( (
int)norms.size() != numvecs ) {
219 <<
"*** ERROR *** MultiVecTraits::MvNorm()." << endl
220 <<
"Method resized the output vector." << endl;
222 for (
int i=0; i<numvecs; i++) {
223 if ( norms[i] < zero_mag ) {
225 <<
"*** ERROR *** MultiVecTraits::Clone()." << endl
226 <<
"Vector had negative norm." << endl;
249 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
250 std::vector<MagType> norms(numvecs), norms2(numvecs);
253 MVT::MvNorm(*B, norms);
254 for (
int i=0; i<numvecs; i++) {
255 if ( norms[i] != zero_mag ) {
257 <<
"*** ERROR *** MultiVecTraits::MvInit() " 258 <<
"and MultiVecTraits::MvNorm()" << endl
259 <<
"Supposedly zero vector has non-zero norm." << endl;
264 MVT::MvNorm(*B, norms);
266 MVT::MvNorm(*B, norms2);
267 for (
int i=0; i<numvecs; i++) {
268 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) {
270 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
271 <<
"Random vector was empty (very unlikely)." << endl;
274 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) {
276 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
277 <<
"Vector had negative norm." << endl;
280 else if ( norms[i] == norms2[i] ) {
282 <<
"*** ERROR *** MutliVecTraits::MvRandom()." << endl
283 <<
"Vectors not random enough." << endl;
298 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
299 std::vector<MagType> norms(numvecs);
302 MVT::MvScale(*B,SCT::zero());
303 MVT::MvNorm(*B, norms);
304 for (
int i=0; i<numvecs; i++) {
305 if ( norms[i] != zero_mag ) {
307 <<
"*** ERROR *** MultiVecTraits::MvScale(alpha) " 308 <<
"Supposedly zero vector has non-zero norm." << endl;
314 std::vector<ScalarType> zeros(numvecs,SCT::zero());
315 MVT::MvScale(*B,zeros);
316 MVT::MvNorm(*B, norms);
317 for (
int i=0; i<numvecs; i++) {
318 if ( norms[i] != zero_mag ) {
320 <<
"*** ERROR *** MultiVecTraits::MvScale(alphas) " 321 <<
"Supposedly zero vector has non-zero norm." << endl;
342 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
343 std::vector<MagType> norms(numvecs);
346 MVT::MvNorm(*B, norms);
347 bool BadNormWarning =
false;
348 for (
int i=0; i<numvecs; i++) {
349 if ( norms[i] < zero_mag ) {
351 <<
"*** ERROR *** MultiVecTraits::MvRandom()." << endl
352 <<
"Vector had negative norm." << endl;
355 else if ( norms[i] != SCT::squareroot(MVT::GetGlobalLength(*B)) && !BadNormWarning ) {
358 <<
"Warning testing MultiVecTraits::MvInit()." << endl
359 <<
"Ones vector should have norm sqrt(dim)." << endl
360 <<
"norms[i]: " << norms[i] <<
"\tdim: " << MVT::GetGlobalLength(*B) << endl << endl;
361 BadNormWarning =
true;
375 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
376 std::vector<MagType> norms(numvecs);
377 MVT::MvInit(*B, zero_mag);
378 MVT::MvNorm(*B, norms);
379 for (
int i=0; i<numvecs; i++) {
380 if ( norms[i] < zero_mag ) {
382 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
383 <<
"Vector had negative norm." << endl;
386 else if ( norms[i] != zero_mag ) {
388 <<
"*** ERROR *** MultiVecTraits::MvInit()." << endl
389 <<
"Zero vector should have norm zero." << endl;
402 Teuchos::RCP<MV> B, C;
403 std::vector<MagType> norms(numvecs), norms2(ind.size());
405 B = MVT::Clone(*A,numvecs);
407 MVT::MvNorm(*B, norms);
408 C = MVT::CloneCopy(*B,ind);
409 MVT::MvNorm(*C, norms2);
410 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
412 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
413 <<
"Wrong number of vectors." << endl;
416 if ( MVT::GetGlobalLength(*C) != MVT::GetGlobalLength(*B) ) {
418 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
419 <<
"Vector lengths don't match." << endl;
422 for (
int i=0; i<numvecs_2; i++) {
423 if ( norms2[i] != norms[ind[i]] ) {
425 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
426 <<
"Copied vectors do not agree:" 427 << norms2[i] <<
" != " << norms[ind[i]] << endl;
431 MVT::MvInit(*B,zero);
432 MVT::MvNorm(*C, norms2);
433 for (
int i=0; i<numvecs_2; i++) {
434 if ( norms2[i] != norms2[i] ) {
436 <<
"*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
437 <<
"Copied vectors were not independent." << endl;
450 Teuchos::RCP<MV> B, C;
451 std::vector<MagType> norms(numvecs), norms2(numvecs);
453 B = MVT::Clone(*A,numvecs);
455 MVT::MvNorm(*B, norms);
456 C = MVT::CloneCopy(*B);
457 MVT::MvNorm(*C, norms2);
458 if ( MVT::GetNumberVecs(*C) != numvecs ) {
460 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
461 <<
"Wrong number of vectors." << endl;
464 for (
int i=0; i<numvecs; i++) {
465 if ( norms2[i] != norms[i] ) {
467 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
468 <<
"Copied vectors do not agree." << endl;
472 MVT::MvInit(*B,zero);
473 MVT::MvNorm(*C, norms);
474 for (
int i=0; i<numvecs; i++) {
475 if ( norms2[i] != norms[i] ) {
477 <<
"*** ERROR *** MultiVecTraits::CloneCopy()." << endl
478 <<
"Copied vectors were not independent." << endl;
492 Teuchos::RCP<MV> B, C;
493 std::vector<MagType> norms(numvecs), norms2(ind.size());
495 B = MVT::Clone(*A,numvecs);
497 MVT::MvNorm(*B, norms);
498 C = MVT::CloneViewNonConst(*B,ind);
499 MVT::MvNorm(*C, norms2);
500 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
502 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
503 <<
"Wrong number of vectors." << endl;
506 for (
int i=0; i<numvecs_2; i++) {
507 if ( norms2[i] != norms[ind[i]] ) {
509 <<
"*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
510 <<
"Viewed vectors do not agree." << endl;
525 Teuchos::RCP<const MV> constB, C;
526 std::vector<MagType> normsB(numvecs), normsC(ind.size());
527 std::vector<int> allind(numvecs);
528 for (
int i=0; i<numvecs; i++) {
532 B = MVT::Clone(*A,numvecs);
534 MVT::MvNorm(*B, normsB);
535 C = MVT::CloneView(*B,ind);
536 MVT::MvNorm(*C, normsC);
537 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
539 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
540 <<
"Wrong number of vectors." << endl;
543 for (
int i=0; i<numvecs_2; i++) {
544 if ( normsC[i] != normsB[ind[i]] ) {
546 <<
"*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
547 <<
"Viewed vectors do not agree." << endl;
566 Teuchos::RCP<MV> B, C;
567 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
568 normsC1(numvecs_2), normsC2(numvecs_2);
570 B = MVT::Clone(*A,numvecs);
571 C = MVT::Clone(*A,numvecs_2);
573 ind.resize(numvecs_2);
574 for (
int i=0; i<numvecs_2; i++) {
580 MVT::MvNorm(*B,normsB1);
581 MVT::MvNorm(*C,normsC1);
582 MVT::SetBlock(*C,ind,*B);
583 MVT::MvNorm(*B,normsB2);
584 MVT::MvNorm(*C,normsC2);
587 for (
int i=0; i<numvecs_2; i++) {
588 if ( normsC1[i] != normsC2[i] ) {
590 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
591 <<
"Operation modified source vectors." << endl;
597 for (
int i=0; i<numvecs; i++) {
600 if ( normsB2[i] != normsC1[i/2] ) {
602 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
603 <<
"Copied vectors do not agree." << endl;
609 if ( normsB1[i] != normsB2[i] ) {
611 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
612 <<
"Incorrect vectors were modified." << endl;
617 MVT::MvInit(*C,zero);
618 MVT::MvNorm(*B,normsB1);
620 for (
int i=0; i<numvecs; i++) {
621 if ( normsB1[i] != normsB2[i] ) {
623 <<
"*** ERROR *** MultiVecTraits::SetBlock()." << endl
624 <<
"Copied vectors were not independent." << endl;
657 Teuchos::RCP<MV> B, C;
658 std::vector<MagType> normsB(p), normsC(q);
659 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
661 B = MVT::Clone(*A,p);
662 C = MVT::Clone(*A,q);
666 MVT::MvNorm(*B,normsB);
668 MVT::MvNorm(*C,normsC);
671 MVT::MvTransMv( zero, *B, *C, SDM );
674 if ( SDM.numRows() != p || SDM.numCols() != q ) {
676 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
677 <<
"Routine resized SerialDenseMatrix." << endl;
682 if ( SDM.normOne() != zero ) {
684 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
685 <<
"Scalar argument processed incorrectly." << endl;
690 MVT::MvTransMv( one, *B, *C, SDM );
694 for (
int i=0; i<p; i++) {
695 for (
int j=0; j<q; j++) {
696 if ( SCT::magnitude(SDM(i,j))
697 > SCT::magnitude(normsB[i]*normsC[j]) ) {
699 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
700 <<
"Triangle inequality did not hold: " 701 << SCT::magnitude(SDM(i,j))
703 << SCT::magnitude(normsB[i]*normsC[j])
711 MVT::MvTransMv( one, *B, *C, SDM );
712 for (
int i=0; i<p; i++) {
713 for (
int j=0; j<q; j++) {
714 if ( SDM(i,j) != zero ) {
716 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
717 <<
"Inner products not zero for C==0." << endl;
724 MVT::MvTransMv( one, *B, *C, SDM );
725 for (
int i=0; i<p; i++) {
726 for (
int j=0; j<q; j++) {
727 if ( SDM(i,j) != zero ) {
729 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
730 <<
"Inner products not zero for B==0." << endl;
743 Teuchos::SerialDenseMatrix<int, ScalarType> largeSDM(p+1,q+1);
744 Teuchos::SerialDenseMatrix<int, ScalarType> SDM2(Teuchos::View, largeSDM, p, q);
745 largeSDM.putScalar( one );
748 MVT::MvTransMv( one, *B, *C, SDM2 );
749 for (
int i=0; i<p; i++) {
750 for (
int j=0; j<q; j++) {
751 if ( SDM2(i,j) != zero ) {
753 <<
"*** ERROR *** MultiVecTraits::MvTransMv()." << endl
754 <<
"Inner products not zero for C==0 when using a view into Teuchos::SerialDenseMatrix<>." << endl;
770 Teuchos::RCP<MV> B, C;
771 std::vector<ScalarType> iprods(p);
772 std::vector<MagType> normsB(p), normsC(p);
774 B = MVT::Clone(*A,p);
775 C = MVT::Clone(*A,p);
779 MVT::MvNorm(*B,normsB);
780 MVT::MvNorm(*C,normsC);
781 MVT::MvDot( *B, *C, iprods );
782 if ( (
int)iprods.size() != p ) {
784 <<
"*** ERROR *** MultiVecTraits::MvDot." << endl
785 <<
"Routine resized results vector." << endl;
788 for (
int i=0; i<p; i++) {
789 if ( SCT::magnitude(iprods[i])
790 > SCT::magnitude(normsB[i]*normsC[i]) ) {
792 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
793 <<
"Inner products not valid." << endl;
799 MVT::MvDot( *B, *C, iprods );
800 for (
int i=0; i<p; i++) {
801 if ( iprods[i] != zero ) {
803 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
804 <<
"Inner products not zero for B==0." << endl;
810 MVT::MvDot( *B, *C, iprods );
811 for (
int i=0; i<p; i++) {
812 if ( iprods[i] != zero ) {
814 <<
"*** ERROR *** MultiVecTraits::MvDot()." << endl
815 <<
"Inner products not zero for C==0." << endl;
831 Teuchos::RCP<MV> B, C, D;
832 std::vector<MagType> normsB1(p), normsB2(p),
833 normsC1(p), normsC2(p),
834 normsD1(p), normsD2(p);
835 ScalarType alpha = SCT::random(),
836 beta = SCT::random();
838 B = MVT::Clone(*A,p);
839 C = MVT::Clone(*A,p);
840 D = MVT::Clone(*A,p);
844 MVT::MvNorm(*B,normsB1);
845 MVT::MvNorm(*C,normsC1);
848 MVT::MvAddMv(zero,*B,one,*C,*D);
849 MVT::MvNorm(*B,normsB2);
850 MVT::MvNorm(*C,normsC2);
851 MVT::MvNorm(*D,normsD1);
852 for (
int i=0; i<p; i++) {
853 if ( normsB1[i] != normsB2[i] ) {
855 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
856 <<
"Input arguments were modified." << endl;
859 else if ( normsC1[i] != normsC2[i] ) {
861 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
862 <<
"Input arguments were modified." << endl;
865 else if ( normsC1[i] != normsD1[i] ) {
867 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
868 <<
"Assignment did not work." << endl;
874 MVT::MvAddMv(one,*B,zero,*C,*D);
875 MVT::MvNorm(*B,normsB2);
876 MVT::MvNorm(*C,normsC2);
877 MVT::MvNorm(*D,normsD1);
878 for (
int i=0; i<p; i++) {
879 if ( normsB1[i] != normsB2[i] ) {
881 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
882 <<
"Input arguments were modified." << endl;
885 else if ( normsC1[i] != normsC2[i] ) {
887 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
888 <<
"Input arguments were modified." << endl;
891 else if ( normsB1[i] != normsD1[i] ) {
893 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
894 <<
"Assignment did not work." << endl;
902 MVT::MvAddMv(alpha,*B,beta,*C,*D);
903 MVT::MvNorm(*B,normsB2);
904 MVT::MvNorm(*C,normsC2);
905 MVT::MvNorm(*D,normsD1);
907 for (
int i=0; i<p; i++) {
908 if ( normsB1[i] != normsB2[i] ) {
910 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
911 <<
"Input arguments were modified." << endl;
914 else if ( normsC1[i] != normsC2[i] ) {
916 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
917 <<
"Input arguments were modified." << endl;
923 MVT::MvAddMv(alpha,*B,beta,*C,*D);
924 MVT::MvNorm(*B,normsB2);
925 MVT::MvNorm(*C,normsC2);
926 MVT::MvNorm(*D,normsD2);
929 for (
int i=0; i<p; i++) {
930 if ( normsB1[i] != normsB2[i] ) {
932 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
933 <<
"Input arguments were modified." << endl;
936 else if ( normsC1[i] != normsC2[i] ) {
938 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
939 <<
"Input arguments were modified." << endl;
942 else if ( normsD1[i] != normsD2[i] ) {
944 <<
"*** ERROR *** MultiVecTraits::MvAddMv()." << endl
945 <<
"Results varies depending on initial state of dest vectors." << endl;
971 Teuchos::RCP<MV> B, D;
972 Teuchos::RCP<const MV> C;
973 std::vector<MagType> normsB(p),
975 std::vector<int> lclindex(p);
976 for (
int i=0; i<p; i++) lclindex[i] = i;
978 B = MVT::Clone(*A,p);
979 C = MVT::CloneView(*B,lclindex);
980 D = MVT::CloneViewNonConst(*B,lclindex);
983 MVT::MvNorm(*B,normsB);
986 MVT::MvAddMv(zero,*B,one,*C,*D);
987 MVT::MvNorm(*D,normsD);
988 for (
int i=0; i<p; i++) {
989 if ( normsB[i] != normsD[i] ) {
991 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
992 <<
"Assignment did not work." << endl;
998 MVT::MvAddMv(one,*B,zero,*C,*D);
999 MVT::MvNorm(*D,normsD);
1000 for (
int i=0; i<p; i++) {
1001 if ( normsB[i] != normsD[i] ) {
1003 <<
"*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1004 <<
"Assignment did not work." << endl;
1022 const int p = 7, q = 5;
1023 Teuchos::RCP<MV> B, C;
1024 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1025 std::vector<MagType> normsC1(q), normsC2(q),
1026 normsB1(p), normsB2(p);
1028 B = MVT::Clone(*A,p);
1029 C = MVT::Clone(*A,q);
1034 MVT::MvNorm(*B,normsB1);
1035 MVT::MvNorm(*C,normsC1);
1037 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1038 MVT::MvNorm(*B,normsB2);
1039 MVT::MvNorm(*C,normsC2);
1040 for (
int i=0; i<p; i++) {
1041 if ( normsB1[i] != normsB2[i] ) {
1043 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1044 <<
"Input vectors were modified." << endl;
1048 for (
int i=0; i<q; i++) {
1049 if ( normsC1[i] != normsC2[i] ) {
1051 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1052 <<
"Arithmetic test 1 failed." << endl;
1060 MVT::MvNorm(*B,normsB1);
1061 MVT::MvNorm(*C,normsC1);
1063 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1064 MVT::MvNorm(*B,normsB2);
1065 MVT::MvNorm(*C,normsC2);
1066 for (
int i=0; i<p; i++) {
1067 if ( normsB1[i] != normsB2[i] ) {
1069 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1070 <<
"Input vectors were modified." << endl;
1074 for (
int i=0; i<q; i++) {
1075 if ( normsC2[i] != zero ) {
1077 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1078 <<
"Arithmetic test 2 failed: " 1091 MVT::MvNorm(*B,normsB1);
1092 MVT::MvNorm(*C,normsC1);
1094 for (
int i=0; i<q; i++) {
1097 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1098 MVT::MvNorm(*B,normsB2);
1099 MVT::MvNorm(*C,normsC2);
1100 for (
int i=0; i<p; i++) {
1101 if ( normsB1[i] != normsB2[i] ) {
1103 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1104 <<
"Input vectors were modified." << endl;
1108 for (
int i=0; i<q; i++) {
1109 if ( normsB1[i] != normsC2[i] ) {
1111 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1112 <<
"Arithmetic test 3 failed: " 1124 MVT::MvNorm(*B,normsB1);
1125 MVT::MvNorm(*C,normsC1);
1127 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1128 MVT::MvNorm(*B,normsB2);
1129 MVT::MvNorm(*C,normsC2);
1130 for (
int i=0; i<p; i++) {
1131 if ( normsB1[i] != normsB2[i] ) {
1133 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1134 <<
"Input vectors were modified." << endl;
1138 for (
int i=0; i<q; i++) {
1139 if ( normsC1[i] != normsC2[i] ) {
1141 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1142 <<
"Arithmetic test 4 failed." << endl;
1158 const int p = 5, q = 7;
1159 Teuchos::RCP<MV> B, C;
1160 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1161 std::vector<MagType> normsC1(q), normsC2(q),
1162 normsB1(p), normsB2(p);
1164 B = MVT::Clone(*A,p);
1165 C = MVT::Clone(*A,q);
1170 MVT::MvNorm(*B,normsB1);
1171 MVT::MvNorm(*C,normsC1);
1173 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1174 MVT::MvNorm(*B,normsB2);
1175 MVT::MvNorm(*C,normsC2);
1176 for (
int i=0; i<p; i++) {
1177 if ( normsB1[i] != normsB2[i] ) {
1179 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1180 <<
"Input vectors were modified." << endl;
1184 for (
int i=0; i<q; i++) {
1185 if ( normsC1[i] != normsC2[i] ) {
1187 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1188 <<
"Arithmetic test 5 failed." << endl;
1196 MVT::MvNorm(*B,normsB1);
1197 MVT::MvNorm(*C,normsC1);
1199 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1200 MVT::MvNorm(*B,normsB2);
1201 MVT::MvNorm(*C,normsC2);
1202 for (
int i=0; i<p; i++) {
1203 if ( normsB1[i] != normsB2[i] ) {
1205 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1206 <<
"Input vectors were modified." << endl;
1210 for (
int i=0; i<q; i++) {
1211 if ( normsC2[i] != zero ) {
1213 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1214 <<
"Arithmetic test 6 failed: " 1226 MVT::MvNorm(*B,normsB1);
1227 MVT::MvNorm(*C,normsC1);
1229 for (
int i=0; i<p; i++) {
1232 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1233 MVT::MvNorm(*B,normsB2);
1234 MVT::MvNorm(*C,normsC2);
1235 for (
int i=0; i<p; i++) {
1236 if ( normsB1[i] != normsB2[i] ) {
1238 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1239 <<
"Input vectors were modified." << endl;
1243 for (
int i=0; i<p; i++) {
1244 if ( normsB1[i] != normsC2[i] ) {
1246 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1247 <<
"Arithmetic test 7 failed." << endl;
1251 for (
int i=p; i<q; i++) {
1252 if ( normsC2[i] != zero ) {
1254 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1255 <<
"Arithmetic test 7 failed." << endl;
1263 MVT::MvNorm(*B,normsB1);
1264 MVT::MvNorm(*C,normsC1);
1266 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1267 MVT::MvNorm(*B,normsB2);
1268 MVT::MvNorm(*C,normsC2);
1269 for (
int i=0; i<p; i++) {
1270 if ( normsB1[i] != normsB2[i] ) {
1272 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1273 <<
"Input vectors were modified." << endl;
1277 for (
int i=0; i<q; i++) {
1278 if ( normsC1[i] != normsC2[i] ) {
1280 <<
"*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1281 <<
"Arithmetic test 8 failed." << endl;
1298 template<
class ScalarType,
class MV,
class OP>
1301 const Teuchos::RCP<const MV> &A,
1302 const Teuchos::RCP<const OP> &M) {
1314 typedef Teuchos::ScalarTraits<ScalarType> SCT;
1316 typedef typename SCT::magnitudeType MagType;
1318 const int numvecs = 10;
1320 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
1321 C = MVT::Clone(*A,numvecs);
1323 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
1324 normsC1(numvecs), normsC2(numvecs);
1325 bool NonDeterministicWarning;
1336 MVT::MvNorm(*B,normsB1);
1337 OPT::Apply(*M,*B,*C);
1338 MVT::MvNorm(*B,normsB2);
1339 MVT::MvNorm(*C,normsC2);
1340 for (
int i=0; i<numvecs; i++) {
1341 if (normsB2[i] != normsB1[i]) {
1343 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1344 <<
"Apply() modified the input vectors." << endl;
1347 if (normsC2[i] != SCT::zero()) {
1349 <<
"*** ERROR *** OperatorTraits::Apply() [1]" << endl
1350 <<
"Operator applied to zero did not return zero." << endl;
1357 MVT::MvNorm(*B,normsB1);
1358 OPT::Apply(*M,*B,*C);
1359 MVT::MvNorm(*B,normsB2);
1360 MVT::MvNorm(*C,normsC2);
1361 bool ZeroWarning =
false;
1362 for (
int i=0; i<numvecs; i++) {
1363 if (normsB2[i] != normsB1[i]) {
1365 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1366 <<
"Apply() modified the input vectors." << endl;
1369 if (normsC2[i] == SCT::zero() && ZeroWarning==
false ) {
1371 <<
"*** ERROR *** OperatorTraits::Apply() [2]" << endl
1372 <<
"Operator applied to random vectors returned zero." << endl;
1379 MVT::MvNorm(*B,normsB1);
1381 OPT::Apply(*M,*B,*C);
1382 MVT::MvNorm(*B,normsB2);
1383 MVT::MvNorm(*C,normsC1);
1384 for (
int i=0; i<numvecs; i++) {
1385 if (normsB2[i] != normsB1[i]) {
1387 <<
"*** ERROR *** OperatorTraits::Apply() [3]" << endl
1388 <<
"Apply() modified the input vectors." << endl;
1399 OPT::Apply(*M,*B,*C);
1400 MVT::MvNorm(*B,normsB2);
1401 MVT::MvNorm(*C,normsC2);
1402 NonDeterministicWarning =
false;
1403 for (
int i=0; i<numvecs; i++) {
1404 if (normsB2[i] != normsB1[i]) {
1406 <<
"*** ERROR *** OperatorTraits::Apply() [4]" << endl
1407 <<
"Apply() modified the input vectors." << endl;
1410 if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) {
1413 <<
"*** WARNING *** OperatorTraits::Apply() [4]" << endl
1414 <<
"Apply() returned two different results." << endl << endl;
1415 NonDeterministicWarning =
true;
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.