00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef ANASAZI_MVOPTESTER_HPP
00030 #define ANASAZI_MVOPTESTER_HPP
00031
00032
00033
00034
00035
00036
00037
00046 #include "AnasaziConfigDefs.hpp"
00047 #include "AnasaziTypes.hpp"
00048
00049 #include "AnasaziMultiVecTraits.hpp"
00050 #include "AnasaziOperatorTraits.hpp"
00051 #include "AnasaziOutputManager.hpp"
00052
00053 #include "Teuchos_RCP.hpp"
00054
00055
00056 namespace Anasazi {
00057
00063 template< class ScalarType, class MV >
00064 bool TestMultiVecTraits(
00065 const Teuchos::RCP<OutputManager<ScalarType> > &om,
00066 const Teuchos::RCP<const MV> &A ) {
00067
00068 using std::endl;
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 typedef MultiVecTraits<ScalarType, MV> MVT;
00142 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00143 typedef typename SCT::magnitudeType MagType;
00144
00145 const ScalarType one = SCT::one();
00146 const ScalarType zero = SCT::zero();
00147 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero();
00148
00149
00150 const unsigned int numvecs = 10;
00151 const unsigned int numvecs_2 = 5;
00152
00153 std::vector<int> ind(numvecs_2);
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 ind[0] = 0;
00164 ind[numvecs_2-1] = numvecs-1;
00165 for (unsigned int i=1; i<numvecs_2-1; i++) {
00166 ind[i] = rand() % numvecs;
00167 }
00168
00169
00170
00171
00172
00173
00174 if ( MVT::GetNumberVecs(*A) <= 0 ) {
00175 om->stream(Warnings)
00176 << "*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl
00177 << "Returned <= 0." << endl;
00178 return false;
00179 }
00180
00181
00182
00183
00184
00185
00186 if ( MVT::GetVecLength(*A) <= 0 ) {
00187 om->stream(Warnings)
00188 << "*** ERROR *** MultiVectorTraits::GetVecLength()" << endl
00189 << "Returned <= 0." << endl;
00190 return false;
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200 {
00201 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
00202 std::vector<MagType> norms(numvecs);
00203 if ( (unsigned int)MVT::GetNumberVecs(*B) != numvecs ) {
00204 om->stream(Warnings)
00205 << "*** ERROR *** MultiVecTraits::Clone()." << endl
00206 << "Did not allocate requested number of vectors." << endl;
00207 return false;
00208 }
00209 if ( MVT::GetVecLength(*B) != MVT::GetVecLength(*A) ) {
00210 om->stream(Warnings)
00211 << "*** ERROR *** MultiVecTraits::Clone()." << endl
00212 << "Did not allocate requested number of vectors." << endl;
00213 return false;
00214 }
00215 MVT::MvNorm(*B, norms);
00216 if ( norms.size() != (unsigned int)numvecs ) {
00217 om->stream(Warnings)
00218 << "*** ERROR *** MultiVecTraits::MvNorm()." << endl
00219 << "Method resized the output vector." << endl;
00220 }
00221 for (unsigned int i=0; i<numvecs; i++) {
00222 if ( norms[i] < zero_mag ) {
00223 om->stream(Warnings)
00224 << "*** ERROR *** MultiVecTraits::Clone()." << endl
00225 << "Vector had negative norm." << endl;
00226 return false;
00227 }
00228 }
00229 }
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247 {
00248 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
00249 std::vector<MagType> norms(numvecs), norms2(numvecs);
00250
00251 MVT::MvInit(*B);
00252 MVT::MvNorm(*B, norms);
00253 for (unsigned int i=0; i<numvecs; i++) {
00254 if ( norms[i] != zero_mag ) {
00255 om->stream(Warnings)
00256 << "*** ERROR *** MultiVecTraits::MvInit() "
00257 << "and MultiVecTraits::MvNorm()" << endl
00258 << "Supposedly zero vector has non-zero norm." << endl;
00259 return false;
00260 }
00261 }
00262 MVT::MvRandom(*B);
00263 MVT::MvNorm(*B, norms);
00264 MVT::MvRandom(*B);
00265 MVT::MvNorm(*B, norms2);
00266 for (unsigned int i=0; i<numvecs; i++) {
00267 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) {
00268 om->stream(Warnings)
00269 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
00270 << "Random vector was empty (very unlikely)." << endl;
00271 return false;
00272 }
00273 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) {
00274 om->stream(Warnings)
00275 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
00276 << "Vector had negative norm." << endl;
00277 return false;
00278 }
00279 else if ( norms[i] == norms2[i] ) {
00280 om->stream(Warnings)
00281 << "*** ERROR *** MutliVecTraits::MvRandom()." << endl
00282 << "Vectors not random enough." << endl;
00283 return false;
00284 }
00285 }
00286 }
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 {
00303 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
00304 std::vector<MagType> norms(numvecs);
00305
00306 MVT::MvInit(*B,one);
00307 MVT::MvNorm(*B, norms);
00308 bool BadNormWarning = false;
00309 for (unsigned int i=0; i<numvecs; i++) {
00310 if ( norms[i] < zero_mag ) {
00311 om->stream(Warnings)
00312 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
00313 << "Vector had negative norm." << endl;
00314 return false;
00315 }
00316 else if ( norms[i] != SCT::squareroot(MVT::GetVecLength(*B)) && !BadNormWarning ) {
00317 om->stream(Warnings)
00318 << endl
00319 << "Warning testing MultiVecTraits::MvInit()." << endl
00320 << "Ones vector should have norm sqrt(dim)." << endl
00321 << "norms[i]: " << norms[i] << "\tdim: " << MVT::GetVecLength(*B) << endl << endl;
00322 BadNormWarning = true;
00323 }
00324 }
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 {
00336 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
00337 std::vector<MagType> norms(numvecs);
00338 MVT::MvInit(*B, zero_mag);
00339 MVT::MvNorm(*B, norms);
00340 for (unsigned int i=0; i<numvecs; i++) {
00341 if ( norms[i] < zero_mag ) {
00342 om->stream(Warnings)
00343 << "*** ERROR *** MultiVecTraits::MvInit()." << endl
00344 << "Vector had negative norm." << endl;
00345 return false;
00346 }
00347 else if ( norms[i] != zero_mag ) {
00348 om->stream(Warnings)
00349 << "*** ERROR *** MultiVecTraits::MvInit()." << endl
00350 << "Zero vector should have norm zero." << endl;
00351 return false;
00352 }
00353 }
00354 }
00355
00356
00357
00358
00359
00360
00361
00362 {
00363 Teuchos::RCP<MV> B, C;
00364 std::vector<MagType> norms(numvecs), norms2(ind.size());
00365
00366 B = MVT::Clone(*A,numvecs);
00367 MVT::MvRandom(*B);
00368 MVT::MvNorm(*B, norms);
00369 C = MVT::CloneCopy(*B,ind);
00370 MVT::MvNorm(*C, norms2);
00371 if ( (unsigned int)MVT::GetNumberVecs(*C) != numvecs_2 ) {
00372 om->stream(Warnings)
00373 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
00374 << "Wrong number of vectors." << endl;
00375 return false;
00376 }
00377 if ( MVT::GetVecLength(*C) != MVT::GetVecLength(*B) ) {
00378 om->stream(Warnings)
00379 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
00380 << "Vector lengths don't match." << endl;
00381 return false;
00382 }
00383 for (unsigned int i=0; i<numvecs_2; i++) {
00384 if ( norms2[i] != norms[ind[i]] ) {
00385 om->stream(Warnings)
00386 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
00387 << "Copied vectors do not agree:"
00388 << norms2[i] << " != " << norms[ind[i]] << endl;
00389 return false;
00390 }
00391 }
00392 MVT::MvInit(*B,zero);
00393 MVT::MvNorm(*C, norms2);
00394 for (unsigned int i=0; i<numvecs_2; i++) {
00395 if ( norms2[i] != norms2[i] ) {
00396 om->stream(Warnings)
00397 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
00398 << "Copied vectors were not independent." << endl;
00399 return false;
00400 }
00401 }
00402 }
00403
00404
00405
00406
00407
00408
00409
00410 {
00411 Teuchos::RCP<MV> B, C;
00412 std::vector<MagType> norms(numvecs), norms2(numvecs);
00413
00414 B = MVT::Clone(*A,numvecs);
00415 MVT::MvRandom(*B);
00416 MVT::MvNorm(*B, norms);
00417 C = MVT::CloneCopy(*B);
00418 MVT::MvNorm(*C, norms2);
00419 if ( (unsigned int)MVT::GetNumberVecs(*C) != numvecs ) {
00420 om->stream(Warnings)
00421 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00422 << "Wrong number of vectors." << endl;
00423 return false;
00424 }
00425 for (unsigned int i=0; i<numvecs; i++) {
00426 if ( norms2[i] != norms[i] ) {
00427 om->stream(Warnings)
00428 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00429 << "Copied vectors do not agree." << endl;
00430 return false;
00431 }
00432 }
00433 MVT::MvInit(*B,zero);
00434 MVT::MvNorm(*C, norms);
00435 for (unsigned int i=0; i<numvecs; i++) {
00436 if ( norms2[i] != norms[i] ) {
00437 om->stream(Warnings)
00438 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00439 << "Copied vectors were not independent." << endl;
00440 return false;
00441 }
00442 }
00443 }
00444
00445
00446
00447
00448
00449
00450
00451
00452 {
00453 Teuchos::RCP<MV> B, C;
00454 std::vector<MagType> norms(numvecs), norms2(ind.size());
00455
00456 B = MVT::Clone(*A,numvecs);
00457 MVT::MvRandom(*B);
00458 MVT::MvNorm(*B, norms);
00459 C = MVT::CloneView(*B,ind);
00460 MVT::MvNorm(*C, norms2);
00461 if ( (unsigned int)MVT::GetNumberVecs(*C) != numvecs_2 ) {
00462 om->stream(Warnings)
00463 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
00464 << "Wrong number of vectors." << endl;
00465 return false;
00466 }
00467 for (unsigned int i=0; i<numvecs_2; i++) {
00468 if ( norms2[i] != norms[ind[i]] ) {
00469 om->stream(Warnings)
00470 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
00471 << "Viewed vectors do not agree." << endl;
00472 return false;
00473 }
00474 }
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 }
00489
00490
00491
00492
00493
00494
00495
00496
00497 {
00498 Teuchos::RCP<MV> B;
00499 Teuchos::RCP<const MV> constB, C;
00500 std::vector<MagType> normsB(numvecs), normsC(ind.size());
00501 std::vector<int> allind(numvecs);
00502 for (unsigned int i=0; i<numvecs; i++) {
00503 allind[i] = i;
00504 }
00505
00506 B = MVT::Clone(*A,numvecs);
00507 MVT::MvRandom( *B );
00508
00509 constB = MVT::CloneView(*B,allind);
00510 MVT::MvNorm(*constB, normsB);
00511 C = MVT::CloneView(*constB,ind);
00512 MVT::MvNorm(*C, normsC);
00513 if ( (unsigned int)MVT::GetNumberVecs(*C) != numvecs_2 ) {
00514 om->stream(Warnings)
00515 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
00516 << "Wrong number of vectors." << endl;
00517 return false;
00518 }
00519 for (unsigned int i=0; i<numvecs_2; i++) {
00520 if ( normsC[i] != normsB[ind[i]] ) {
00521 om->stream(Warnings)
00522 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
00523 << "Viewed vectors do not agree." << endl;
00524 return false;
00525 }
00526 }
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540 }
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554 {
00555 Teuchos::RCP<MV> B, C;
00556 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
00557 normsC1(numvecs_2), normsC2(numvecs_2);
00558
00559 B = MVT::Clone(*A,numvecs);
00560 C = MVT::Clone(*A,numvecs_2);
00561
00562 ind.resize(numvecs_2);
00563 for (unsigned int i=0; i<numvecs_2; i++) {
00564 ind[i] = 2*i;
00565 }
00566 MVT::MvRandom(*B);
00567 MVT::MvRandom(*C);
00568
00569 MVT::MvNorm(*B,normsB1);
00570 MVT::MvNorm(*C,normsC1);
00571 MVT::SetBlock(*C,ind,*B);
00572 MVT::MvNorm(*B,normsB2);
00573 MVT::MvNorm(*C,normsC2);
00574
00575
00576 for (unsigned int i=0; i<numvecs_2; i++) {
00577 if ( normsC1[i] != normsC2[i] ) {
00578 om->stream(Warnings)
00579 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00580 << "Operation modified source vectors." << endl;
00581 return false;
00582 }
00583 }
00584
00585
00586 for (unsigned int i=0; i<numvecs; i++) {
00587 if (i % 2 == 0) {
00588
00589 if ( normsB2[i] != normsC1[i/2] ) {
00590 om->stream(Warnings)
00591 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00592 << "Copied vectors do not agree." << endl;
00593 return false;
00594 }
00595 }
00596 else {
00597
00598 if ( normsB1[i] != normsB2[i] ) {
00599 om->stream(Warnings)
00600 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00601 << "Incorrect vectors were modified." << endl;
00602 return false;
00603 }
00604 }
00605 }
00606 MVT::MvInit(*C,zero);
00607 MVT::MvNorm(*B,normsB1);
00608
00609 for (unsigned int i=0; i<numvecs; i++) {
00610 if ( normsB1[i] != normsB2[i] ) {
00611 om->stream(Warnings)
00612 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00613 << "Copied vectors were not independent." << endl;
00614 return false;
00615 }
00616 }
00617 }
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639 {
00640 const unsigned int p = 7;
00641 const unsigned int q = 9;
00642 Teuchos::RCP<MV> B, C;
00643 std::vector<MagType> normsB(p), normsC(q);
00644 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
00645
00646 B = MVT::Clone(*A,p);
00647 C = MVT::Clone(*A,q);
00648
00649
00650 MVT::MvRandom(*B);
00651 MVT::MvNorm(*B,normsB);
00652 MVT::MvRandom(*C);
00653 MVT::MvNorm(*C,normsC);
00654
00655
00656 MVT::MvTransMv( zero, *B, *C, SDM );
00657
00658
00659 if ( (unsigned int)SDM.numRows() != p || (unsigned int)SDM.numCols() != q ) {
00660 om->stream(Warnings)
00661 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00662 << "Routine resized SerialDenseMatrix." << endl;
00663 return false;
00664 }
00665
00666
00667 if ( SDM.normOne() != zero ) {
00668 om->stream(Warnings)
00669 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00670 << "Scalar argument processed incorrectly." << endl;
00671 return false;
00672 }
00673
00674
00675 MVT::MvTransMv( one, *B, *C, SDM );
00676
00677
00678
00679 for (unsigned int i=0; i<p; i++) {
00680 for (unsigned int j=0; j<q; j++) {
00681 if ( SCT::magnitude(SDM(i,j))
00682 > SCT::magnitude(normsB[i]*normsC[j]) ) {
00683 om->stream(Warnings)
00684 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00685 << "Triangle inequality did not hold: "
00686 << SCT::magnitude(SDM(i,j))
00687 << " > "
00688 << SCT::magnitude(normsB[i]*normsC[j])
00689 << endl;
00690 return false;
00691 }
00692 }
00693 }
00694 MVT::MvInit(*C);
00695 MVT::MvRandom(*B);
00696 MVT::MvTransMv( one, *B, *C, SDM );
00697 for (unsigned int i=0; i<p; i++) {
00698 for (unsigned int j=0; j<q; j++) {
00699 if ( SDM(i,j) != zero ) {
00700 om->stream(Warnings)
00701 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00702 << "Inner products not zero for C==0." << endl;
00703 return false;
00704 }
00705 }
00706 }
00707 MVT::MvInit(*B);
00708 MVT::MvRandom(*C);
00709 MVT::MvTransMv( one, *B, *C, SDM );
00710 for (unsigned int i=0; i<p; i++) {
00711 for (unsigned int j=0; j<q; j++) {
00712 if ( SDM(i,j) != zero ) {
00713 om->stream(Warnings)
00714 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00715 << "Inner products not zero for B==0." << endl;
00716 return false;
00717 }
00718 }
00719 }
00720 }
00721
00722
00723
00724
00725
00726
00727
00728
00729 {
00730 const unsigned int p = 7;
00731 Teuchos::RCP<MV> B, C;
00732 std::vector<ScalarType> iprods(p);
00733 std::vector<MagType> normsB(p), normsC(p);
00734
00735 B = MVT::Clone(*A,p);
00736 C = MVT::Clone(*A,p);
00737
00738 MVT::MvRandom(*B);
00739 MVT::MvRandom(*C);
00740 MVT::MvNorm(*B,normsB);
00741 MVT::MvNorm(*C,normsC);
00742 MVT::MvDot( *B, *C, iprods );
00743 if ( iprods.size() != p ) {
00744 om->stream(Warnings)
00745 << "*** ERROR *** MultiVecTraits::MvDot." << endl
00746 << "Routine resized results vector." << endl;
00747 return false;
00748 }
00749 for (unsigned int i=0; i<p; i++) {
00750 if ( SCT::magnitude(iprods[i])
00751 > SCT::magnitude(normsB[i]*normsC[i]) ) {
00752 om->stream(Warnings)
00753 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00754 << "Inner products not valid." << endl;
00755 return false;
00756 }
00757 }
00758 MVT::MvInit(*B);
00759 MVT::MvRandom(*C);
00760 MVT::MvDot( *B, *C, iprods );
00761 for (unsigned int i=0; i<p; i++) {
00762 if ( iprods[i] != zero ) {
00763 om->stream(Warnings)
00764 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00765 << "Inner products not zero for B==0." << endl;
00766 return false;
00767 }
00768 }
00769 MVT::MvInit(*C);
00770 MVT::MvRandom(*B);
00771 MVT::MvDot( *B, *C, iprods );
00772 for (unsigned int i=0; i<p; i++) {
00773 if ( iprods[i] != zero ) {
00774 om->stream(Warnings)
00775 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00776 << "Inner products not zero for C==0." << endl;
00777 return false;
00778 }
00779 }
00780 }
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790 {
00791 const unsigned int p = 7;
00792 Teuchos::RCP<MV> B, C, D;
00793 std::vector<MagType> normsB1(p), normsB2(p),
00794 normsC1(p), normsC2(p),
00795 normsD1(p), normsD2(p);
00796 ScalarType alpha = SCT::random(),
00797 beta = SCT::random();
00798
00799 B = MVT::Clone(*A,p);
00800 C = MVT::Clone(*A,p);
00801 D = MVT::Clone(*A,p);
00802
00803 MVT::MvRandom(*B);
00804 MVT::MvRandom(*C);
00805 MVT::MvNorm(*B,normsB1);
00806 MVT::MvNorm(*C,normsC1);
00807
00808
00809 MVT::MvAddMv(zero,*B,one,*C,*D);
00810 MVT::MvNorm(*B,normsB2);
00811 MVT::MvNorm(*C,normsC2);
00812 MVT::MvNorm(*D,normsD1);
00813 for (unsigned int i=0; i<p; i++) {
00814 if ( normsB1[i] != normsB2[i] ) {
00815 om->stream(Warnings)
00816 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00817 << "Input arguments were modified." << endl;
00818 return false;
00819 }
00820 else if ( normsC1[i] != normsC2[i] ) {
00821 om->stream(Warnings)
00822 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00823 << "Input arguments were modified." << endl;
00824 return false;
00825 }
00826 else if ( normsC1[i] != normsD1[i] ) {
00827 om->stream(Warnings)
00828 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00829 << "Assignment did not work." << endl;
00830 return false;
00831 }
00832 }
00833
00834
00835 MVT::MvAddMv(one,*B,zero,*C,*D);
00836 MVT::MvNorm(*B,normsB2);
00837 MVT::MvNorm(*C,normsC2);
00838 MVT::MvNorm(*D,normsD1);
00839 for (unsigned int i=0; i<p; i++) {
00840 if ( normsB1[i] != normsB2[i] ) {
00841 om->stream(Warnings)
00842 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00843 << "Input arguments were modified." << endl;
00844 return false;
00845 }
00846 else if ( normsC1[i] != normsC2[i] ) {
00847 om->stream(Warnings)
00848 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00849 << "Input arguments were modified." << endl;
00850 return false;
00851 }
00852 else if ( normsB1[i] != normsD1[i] ) {
00853 om->stream(Warnings)
00854 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00855 << "Assignment did not work." << endl;
00856 return false;
00857 }
00858 }
00859
00860
00861
00862 MVT::MvRandom(*D);
00863 MVT::MvAddMv(alpha,*B,beta,*C,*D);
00864 MVT::MvNorm(*B,normsB2);
00865 MVT::MvNorm(*C,normsC2);
00866 MVT::MvNorm(*D,normsD1);
00867
00868 for (unsigned int i=0; i<p; i++) {
00869 if ( normsB1[i] != normsB2[i] ) {
00870 om->stream(Warnings)
00871 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00872 << "Input arguments were modified." << endl;
00873 return false;
00874 }
00875 else if ( normsC1[i] != normsC2[i] ) {
00876 om->stream(Warnings)
00877 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00878 << "Input arguments were modified." << endl;
00879 return false;
00880 }
00881 }
00882
00883 MVT::MvInit(*D);
00884 MVT::MvAddMv(alpha,*B,beta,*C,*D);
00885 MVT::MvNorm(*B,normsB2);
00886 MVT::MvNorm(*C,normsC2);
00887 MVT::MvNorm(*D,normsD2);
00888
00889
00890 for (unsigned int i=0; i<p; i++) {
00891 if ( normsB1[i] != normsB2[i] ) {
00892 om->stream(Warnings)
00893 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00894 << "Input arguments were modified." << endl;
00895 return false;
00896 }
00897 else if ( normsC1[i] != normsC2[i] ) {
00898 om->stream(Warnings)
00899 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00900 << "Input arguments were modified." << endl;
00901 return false;
00902 }
00903 else if ( normsD1[i] != normsD2[i] ) {
00904 om->stream(Warnings)
00905 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00906 << "Results varies depending on initial state of dest vectors." << endl;
00907 return false;
00908 }
00909 }
00910 }
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930 {
00931 const unsigned int p = 7;
00932 Teuchos::RCP<MV> B, C, D;
00933 std::vector<MagType> normsB(p),
00934 normsD(p);
00935 std::vector<int> lclindex(p);
00936 for (unsigned int i=0; i<p; i++) lclindex[i] = i;
00937
00938 B = MVT::Clone(*A,p);
00939 C = MVT::CloneView(*B,lclindex);
00940 D = MVT::CloneView(*B,lclindex);
00941
00942 MVT::MvRandom(*B);
00943 MVT::MvNorm(*B,normsB);
00944
00945
00946 MVT::MvAddMv(zero,*B,one,*C,*D);
00947 MVT::MvNorm(*D,normsD);
00948 for (unsigned int i=0; i<p; i++) {
00949 if ( normsB[i] != normsD[i] ) {
00950 om->stream(Warnings)
00951 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
00952 << "Assignment did not work." << endl;
00953 return false;
00954 }
00955 }
00956
00957
00958 MVT::MvAddMv(one,*B,zero,*C,*D);
00959 MVT::MvNorm(*D,normsD);
00960 for (unsigned int i=0; i<p; i++) {
00961 if ( normsB[i] != normsD[i] ) {
00962 om->stream(Warnings)
00963 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
00964 << "Assignment did not work." << endl;
00965 return false;
00966 }
00967 }
00968
00969 }
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981 {
00982 const unsigned int p = 7, q = 5;
00983 Teuchos::RCP<MV> B, C;
00984 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
00985 std::vector<MagType> normsC1(q), normsC2(q),
00986 normsB1(p), normsB2(p);
00987
00988 B = MVT::Clone(*A,p);
00989 C = MVT::Clone(*A,q);
00990
00991
00992 MVT::MvRandom(*B);
00993 MVT::MvRandom(*C);
00994 MVT::MvNorm(*B,normsB1);
00995 MVT::MvNorm(*C,normsC1);
00996 SDM.random();
00997 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
00998 MVT::MvNorm(*B,normsB2);
00999 MVT::MvNorm(*C,normsC2);
01000 for (unsigned int i=0; i<p; i++) {
01001 if ( normsB1[i] != normsB2[i] ) {
01002 om->stream(Warnings)
01003 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01004 << "Input vectors were modified." << endl;
01005 return false;
01006 }
01007 }
01008 for (unsigned int i=0; i<q; i++) {
01009 if ( normsC1[i] != normsC2[i] ) {
01010 om->stream(Warnings)
01011 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01012 << "Arithmetic test 1 failed." << endl;
01013 return false;
01014 }
01015 }
01016
01017
01018 MVT::MvRandom(*B);
01019 MVT::MvRandom(*C);
01020 MVT::MvNorm(*B,normsB1);
01021 MVT::MvNorm(*C,normsC1);
01022 SDM.random();
01023 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01024 MVT::MvNorm(*B,normsB2);
01025 MVT::MvNorm(*C,normsC2);
01026 for (unsigned int i=0; i<p; i++) {
01027 if ( normsB1[i] != normsB2[i] ) {
01028 om->stream(Warnings)
01029 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01030 << "Input vectors were modified." << endl;
01031 return false;
01032 }
01033 }
01034 for (unsigned int i=0; i<q; i++) {
01035 if ( normsC2[i] != zero ) {
01036 om->stream(Warnings)
01037 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01038 << "Arithmetic test 2 failed: "
01039 << normsC2[i]
01040 << " != "
01041 << zero
01042 << endl;
01043 return false;
01044 }
01045 }
01046
01047
01048
01049 MVT::MvRandom(*B);
01050 MVT::MvRandom(*C);
01051 MVT::MvNorm(*B,normsB1);
01052 MVT::MvNorm(*C,normsC1);
01053 SDM.scale(zero);
01054 for (unsigned int i=0; i<q; i++) {
01055 SDM(i,i) = one;
01056 }
01057 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01058 MVT::MvNorm(*B,normsB2);
01059 MVT::MvNorm(*C,normsC2);
01060 for (unsigned int i=0; i<p; i++) {
01061 if ( normsB1[i] != normsB2[i] ) {
01062 om->stream(Warnings)
01063 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01064 << "Input vectors were modified." << endl;
01065 return false;
01066 }
01067 }
01068 for (unsigned int i=0; i<q; i++) {
01069 if ( normsB1[i] != normsC2[i] ) {
01070 om->stream(Warnings)
01071 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01072 << "Arithmetic test 3 failed: "
01073 << normsB1[i]
01074 << " != "
01075 << normsC2[i]
01076 << endl;
01077 return false;
01078 }
01079 }
01080
01081
01082 MVT::MvRandom(*B);
01083 MVT::MvRandom(*C);
01084 MVT::MvNorm(*B,normsB1);
01085 MVT::MvNorm(*C,normsC1);
01086 SDM.scale(zero);
01087 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01088 MVT::MvNorm(*B,normsB2);
01089 MVT::MvNorm(*C,normsC2);
01090 for (unsigned int i=0; i<p; i++) {
01091 if ( normsB1[i] != normsB2[i] ) {
01092 om->stream(Warnings)
01093 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01094 << "Input vectors were modified." << endl;
01095 return false;
01096 }
01097 }
01098 for (unsigned int i=0; i<q; i++) {
01099 if ( normsC1[i] != normsC2[i] ) {
01100 om->stream(Warnings)
01101 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01102 << "Arithmetic test 4 failed." << endl;
01103 return false;
01104 }
01105 }
01106 }
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117 {
01118 const unsigned int p = 5, q = 7;
01119 Teuchos::RCP<MV> B, C;
01120 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
01121 std::vector<MagType> normsC1(q), normsC2(q),
01122 normsB1(p), normsB2(p);
01123
01124 B = MVT::Clone(*A,p);
01125 C = MVT::Clone(*A,q);
01126
01127
01128 MVT::MvRandom(*B);
01129 MVT::MvRandom(*C);
01130 MVT::MvNorm(*B,normsB1);
01131 MVT::MvNorm(*C,normsC1);
01132 SDM.random();
01133 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
01134 MVT::MvNorm(*B,normsB2);
01135 MVT::MvNorm(*C,normsC2);
01136 for (unsigned int i=0; i<p; i++) {
01137 if ( normsB1[i] != normsB2[i] ) {
01138 om->stream(Warnings)
01139 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01140 << "Input vectors were modified." << endl;
01141 return false;
01142 }
01143 }
01144 for (unsigned int i=0; i<q; i++) {
01145 if ( normsC1[i] != normsC2[i] ) {
01146 om->stream(Warnings)
01147 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01148 << "Arithmetic test 5 failed." << endl;
01149 return false;
01150 }
01151 }
01152
01153
01154 MVT::MvRandom(*B);
01155 MVT::MvRandom(*C);
01156 MVT::MvNorm(*B,normsB1);
01157 MVT::MvNorm(*C,normsC1);
01158 SDM.random();
01159 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01160 MVT::MvNorm(*B,normsB2);
01161 MVT::MvNorm(*C,normsC2);
01162 for (unsigned int i=0; i<p; i++) {
01163 if ( normsB1[i] != normsB2[i] ) {
01164 om->stream(Warnings)
01165 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01166 << "Input vectors were modified." << endl;
01167 return false;
01168 }
01169 }
01170 for (unsigned int i=0; i<q; i++) {
01171 if ( normsC2[i] != zero ) {
01172 om->stream(Warnings)
01173 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01174 << "Arithmetic test 6 failed: "
01175 << normsC2[i]
01176 << " != "
01177 << zero
01178 << endl;
01179 return false;
01180 }
01181 }
01182
01183
01184 MVT::MvRandom(*B);
01185 MVT::MvRandom(*C);
01186 MVT::MvNorm(*B,normsB1);
01187 MVT::MvNorm(*C,normsC1);
01188 SDM.scale(zero);
01189 for (unsigned int i=0; i<p; i++) {
01190 SDM(i,i) = one;
01191 }
01192 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01193 MVT::MvNorm(*B,normsB2);
01194 MVT::MvNorm(*C,normsC2);
01195 for (unsigned int i=0; i<p; i++) {
01196 if ( normsB1[i] != normsB2[i] ) {
01197 om->stream(Warnings)
01198 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01199 << "Input vectors were modified." << endl;
01200 return false;
01201 }
01202 }
01203 for (unsigned int i=0; i<p; i++) {
01204 if ( normsB1[i] != normsC2[i] ) {
01205 om->stream(Warnings)
01206 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01207 << "Arithmetic test 7 failed." << endl;
01208 return false;
01209 }
01210 }
01211 for (unsigned int i=p; i<q; i++) {
01212 if ( normsC2[i] != zero ) {
01213 om->stream(Warnings)
01214 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01215 << "Arithmetic test 7 failed." << endl;
01216 return false;
01217 }
01218 }
01219
01220
01221 MVT::MvRandom(*B);
01222 MVT::MvRandom(*C);
01223 MVT::MvNorm(*B,normsB1);
01224 MVT::MvNorm(*C,normsC1);
01225 SDM.scale(zero);
01226 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01227 MVT::MvNorm(*B,normsB2);
01228 MVT::MvNorm(*C,normsC2);
01229 for (unsigned int i=0; i<p; i++) {
01230 if ( normsB1[i] != normsB2[i] ) {
01231 om->stream(Warnings)
01232 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01233 << "Input vectors were modified." << endl;
01234 return false;
01235 }
01236 }
01237 for (unsigned int i=0; i<q; i++) {
01238 if ( normsC1[i] != normsC2[i] ) {
01239 om->stream(Warnings)
01240 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01241 << "Arithmetic test 8 failed." << endl;
01242 return false;
01243 }
01244 }
01245 }
01246
01247 return true;
01248
01249 }
01250
01251
01252
01258 template< class ScalarType, class MV, class OP>
01259 bool TestOperatorTraits(
01260 const Teuchos::RCP<OutputManager<ScalarType> > &om,
01261 const Teuchos::RCP<const MV> &A,
01262 const Teuchos::RCP<const OP> &M) {
01263
01264 using std::endl;
01265
01266
01267
01268
01269
01270
01271
01272
01273 typedef MultiVecTraits<ScalarType, MV> MVT;
01274 typedef Teuchos::ScalarTraits<ScalarType> SCT;
01275 typedef OperatorTraits<ScalarType, MV, OP> OPT;
01276 typedef typename SCT::magnitudeType MagType;
01277
01278 const unsigned int numvecs = 10;
01279
01280 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
01281 C = MVT::Clone(*A,numvecs);
01282
01283 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
01284 normsC1(numvecs), normsC2(numvecs);
01285 bool NonDeterministicWarning;
01286
01287
01288
01289
01290
01291
01292
01293
01294 MVT::MvInit(*B);
01295 MVT::MvRandom(*C);
01296 MVT::MvNorm(*B,normsB1);
01297 OPT::Apply(*M,*B,*C);
01298 MVT::MvNorm(*B,normsB2);
01299 MVT::MvNorm(*C,normsC2);
01300 for (unsigned int i=0; i<numvecs; i++) {
01301 if (normsB2[i] != normsB1[i]) {
01302 om->stream(Warnings)
01303 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01304 << "Apply() modified the input vectors." << endl;
01305 return false;
01306 }
01307 if (normsC2[i] != SCT::zero()) {
01308 om->stream(Warnings)
01309 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01310 << "Operator applied to zero did not return zero." << endl;
01311 return false;
01312 }
01313 }
01314
01315
01316 MVT::MvRandom(*B);
01317 MVT::MvNorm(*B,normsB1);
01318 OPT::Apply(*M,*B,*C);
01319 MVT::MvNorm(*B,normsB2);
01320 MVT::MvNorm(*C,normsC2);
01321 bool ZeroWarning = false;
01322 for (unsigned int i=0; i<numvecs; i++) {
01323 if (normsB2[i] != normsB1[i]) {
01324 om->stream(Warnings)
01325 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01326 << "Apply() modified the input vectors." << endl;
01327 return false;
01328 }
01329 if (normsC2[i] == SCT::zero() && ZeroWarning==false ) {
01330 om->stream(Warnings)
01331 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01332 << "Operator applied to random vectors returned zero." << endl;
01333 ZeroWarning = true;
01334 }
01335 }
01336
01337
01338 MVT::MvRandom(*B);
01339 MVT::MvNorm(*B,normsB1);
01340 MVT::MvInit(*C);
01341 OPT::Apply(*M,*B,*C);
01342 MVT::MvNorm(*B,normsB2);
01343 MVT::MvNorm(*C,normsC1);
01344 for (unsigned int i=0; i<numvecs; i++) {
01345 if (normsB2[i] != normsB1[i]) {
01346 om->stream(Warnings)
01347 << "*** ERROR *** OperatorTraits::Apply() [3]" << endl
01348 << "Apply() modified the input vectors." << endl;
01349 return false;
01350 }
01351 }
01352
01353
01354
01355
01356
01357
01358 MVT::MvRandom(*C);
01359 OPT::Apply(*M,*B,*C);
01360 MVT::MvNorm(*B,normsB2);
01361 MVT::MvNorm(*C,normsC2);
01362 NonDeterministicWarning = false;
01363 for (unsigned int i=0; i<numvecs; i++) {
01364 if (normsB2[i] != normsB1[i]) {
01365 om->stream(Warnings)
01366 << "*** ERROR *** OperatorTraits::Apply() [4]" << endl
01367 << "Apply() modified the input vectors." << endl;
01368 return false;
01369 }
01370 if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) {
01371 om->stream(Warnings)
01372 << endl
01373 << "*** WARNING *** OperatorTraits::Apply() [4]" << endl
01374 << "Apply() returned two different results." << endl << endl;
01375 NonDeterministicWarning = true;
01376 }
01377 }
01378
01379 return true;
01380
01381 }
01382
01383 }
01384
01385 #endif