|
Anasazi Version of the Day
|
|
Anasazi Version of the Day
|
00001 // @HEADER 00002 // *********************************************************************** 00003 // 00004 // Anasazi: Block Eigensolvers Package 00005 // Copyright (2004) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00025 // 00026 // *********************************************************************** 00027 // @HEADER 00028 // 00029 #ifndef ANASAZI_MVOPTESTER_HPP 00030 #define ANASAZI_MVOPTESTER_HPP 00031 00032 // Assumptions that I have made: 00033 // * I assume/verify that a multivector must have at least one vector. This seems 00034 // to be consistent with Epetra_MultiVec. 00035 // * I do not assume that an operator is deterministic; I do assume that the 00036 // operator, applied to 0, will return 0. 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 namespace Anasazi { 00056 00062 template< class ScalarType, class MV > 00063 bool TestMultiVecTraits( 00064 const Teuchos::RCP<OutputManager<ScalarType> > &om, 00065 const Teuchos::RCP<const MV> &A ) { 00066 00067 using std::endl; 00068 00069 /* MVT Contract: 00070 00071 Clone(MV,int) 00072 CloneCopy(MV) 00073 CloneCopy(MV,vector<int>) 00074 USER: will request positive number of vectors 00075 MV: will return a multivector with exactly the number of 00076 requested vectors. 00077 vectors are the same dimension as the cloned MV 00078 00079 00080 CloneView(MV,vector<int>) [const and non-const] 00081 USER: There is no assumed communication between creation and 00082 destruction of a view. I.e., after a view is created, changes to the 00083 source multivector are not reflected in the view. Likewise, until 00084 destruction of the view, changes in the view are not reflected in the 00085 source multivector. 00086 00087 GetVecLength 00088 MV: will always be positive (MV cannot have zero vectors) 00089 00090 GetNumberVecs 00091 MV: will always be positive (MV cannot have zero vectors) 00092 00093 MvAddMv 00094 USER: multivecs will be of the same dimension and same number of vecs 00095 MV: input vectors will not be modified 00096 performing C=0*A+1*B will assign B to C exactly 00097 00098 MvTimesMatAddMv 00099 USER: multivecs and serialdensematrix will be of the proper shape 00100 MV: input arguments will not be modified 00101 following arithmetic relations hold exactly: 00102 A*I = A 00103 0*B = B 00104 1*B = B 00105 00106 MvTransMv 00107 USER: SerialDenseMatrix will be large enough to hold results. 00108 MV: SerialDenseMatrix will not be resized. 00109 Inner products will satisfy |a'*b| <= |a|*|b| 00110 alpha == 0 => SerialDenseMatrix == 0 00111 00112 MvDot 00113 USER: Results vector will be large enough for results. 00114 Both multivectors will have the same number of vectors. 00115 (Epetra crashes, otherwise.) 00116 MV: Inner products will satisfy |a'*b| <= |a|*|b| 00117 Results vector will not be resized. 00118 00119 MvNorm 00120 MV: vector norm is always non-negative, and zero 00121 only for zero vectors. 00122 results vector should not be resized 00123 00124 SetBlock 00125 USER: indices will be distinct 00126 MV: assigns copies of the vectors to the specified 00127 locations, leaving the other vectors untouched. 00128 00129 MvRandom 00130 MV: Generate zero vector with "zero" probability 00131 Don't gen the same vectors twice. 00132 00133 MvInit 00134 MV: Init(alpha) sets all elements to alpha 00135 00136 MvScale (two versions) 00137 MV: scales multivector values 00138 00139 MvPrint 00140 MV: routine does not modify vectors (not tested here) 00141 *********************************************************************/ 00142 00143 typedef MultiVecTraits<ScalarType, MV> MVT; 00144 typedef Teuchos::ScalarTraits<ScalarType> SCT; 00145 typedef typename SCT::magnitudeType MagType; 00146 00147 const ScalarType one = SCT::one(); 00148 const ScalarType zero = SCT::zero(); 00149 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero(); 00150 00151 // Don't change these two without checking the initialization of ind below 00152 const int numvecs = 10; 00153 const int numvecs_2 = 5; 00154 00155 std::vector<int> ind(numvecs_2); 00156 00157 /* Initialize indices for selected copies/views 00158 The MVT specialization should not assume that 00159 these are ordered or even distinct. 00160 Also retrieve the edges. 00161 00162 However, to spice things up, grab the first vector, 00163 last vector, and choose the others randomly. 00164 */ 00165 TEUCHOS_TEST_FOR_EXCEPT(numvecs_2 != 5); 00166 ind[0] = 0; 00167 ind[1] = 5; 00168 ind[2] = 2; 00169 ind[3] = 2; 00170 ind[4] = 9; 00171 00172 /*********** GetNumberVecs() ***************************************** 00173 Verify: 00174 1) This number should be strictly positive 00175 *********************************************************************/ 00176 if ( MVT::GetNumberVecs(*A) <= 0 ) { 00177 om->stream(Warnings) 00178 << "*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl 00179 << "Returned <= 0." << endl; 00180 return false; 00181 } 00182 00183 00184 /*********** GetVecLength() ****************************************** 00185 Verify: 00186 1) This number should be strictly positive 00187 *********************************************************************/ 00188 if ( MVT::GetVecLength(*A) <= 0 ) { 00189 om->stream(Warnings) 00190 << "*** ERROR *** MultiVectorTraits::GetVecLength()" << endl 00191 << "Returned <= 0." << endl; 00192 return false; 00193 } 00194 00195 00196 /*********** Clone() and MvNorm() ************************************ 00197 Verify: 00198 1) Clone() allows us to specify the number of vectors 00199 2) Clone() returns a multivector of the same dimension 00200 3) Vector norms shouldn't be negative 00201 *********************************************************************/ 00202 { 00203 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs); 00204 std::vector<MagType> norms(numvecs); 00205 if ( MVT::GetNumberVecs(*B) != numvecs ) { 00206 om->stream(Warnings) 00207 << "*** ERROR *** MultiVecTraits::Clone()." << endl 00208 << "Did not allocate requested number of vectors." << endl; 00209 return false; 00210 } 00211 if ( MVT::GetVecLength(*B) != MVT::GetVecLength(*A) ) { 00212 om->stream(Warnings) 00213 << "*** ERROR *** MultiVecTraits::Clone()." << endl 00214 << "Did not allocate requested number of vectors." << endl; 00215 return false; 00216 } 00217 MVT::MvNorm(*B, norms); 00218 if ( (int)norms.size() != numvecs ) { 00219 om->stream(Warnings) 00220 << "*** ERROR *** MultiVecTraits::MvNorm()." << endl 00221 << "Method resized the output vector." << endl; 00222 } 00223 for (int i=0; i<numvecs; i++) { 00224 if ( norms[i] < zero_mag ) { 00225 om->stream(Warnings) 00226 << "*** ERROR *** MultiVecTraits::Clone()." << endl 00227 << "Vector had negative norm." << endl; 00228 return false; 00229 } 00230 } 00231 } 00232 00233 00234 /*********** MvRandom() and MvNorm() and MvInit() ******************** 00235 Verify: 00236 1) Set vectors to zero 00237 2) Check that norm is zero 00238 3) Perform MvRandom. 00239 4) Verify that vectors aren't zero anymore 00240 5) Perform MvRandom again. 00241 6) Verify that vector norms are different than before 00242 00243 Without knowing something about the random distribution, 00244 this is about the best that we can do, to make sure that MvRandom 00245 did at least *something*. 00246 00247 Also, make sure vector norms aren't negative. 00248 *********************************************************************/ 00249 { 00250 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs); 00251 std::vector<MagType> norms(numvecs), norms2(numvecs); 00252 00253 MVT::MvInit(*B); 00254 MVT::MvNorm(*B, norms); 00255 for (int i=0; i<numvecs; i++) { 00256 if ( norms[i] != zero_mag ) { 00257 om->stream(Warnings) 00258 << "*** ERROR *** MultiVecTraits::MvInit() " 00259 << "and MultiVecTraits::MvNorm()" << endl 00260 << "Supposedly zero vector has non-zero norm." << endl; 00261 return false; 00262 } 00263 } 00264 MVT::MvRandom(*B); 00265 MVT::MvNorm(*B, norms); 00266 MVT::MvRandom(*B); 00267 MVT::MvNorm(*B, norms2); 00268 for (int i=0; i<numvecs; i++) { 00269 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) { 00270 om->stream(Warnings) 00271 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl 00272 << "Random vector was empty (very unlikely)." << endl; 00273 return false; 00274 } 00275 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) { 00276 om->stream(Warnings) 00277 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl 00278 << "Vector had negative norm." << endl; 00279 return false; 00280 } 00281 else if ( norms[i] == norms2[i] ) { 00282 om->stream(Warnings) 00283 << "*** ERROR *** MutliVecTraits::MvRandom()." << endl 00284 << "Vectors not random enough." << endl; 00285 return false; 00286 } 00287 } 00288 } 00289 00290 00291 /*********** MvRandom() and MvNorm() and MvScale() ******************* 00292 Verify: 00293 1) Perform MvRandom. 00294 2) Verify that vectors aren't zero 00295 3) Set vectors to zero via MvScale 00296 4) Check that norm is zero 00297 *********************************************************************/ 00298 { 00299 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs); 00300 std::vector<MagType> norms(numvecs); 00301 00302 MVT::MvRandom(*B); 00303 MVT::MvScale(*B,SCT::zero()); 00304 MVT::MvNorm(*B, norms); 00305 for (int i=0; i<numvecs; i++) { 00306 if ( norms[i] != zero_mag ) { 00307 om->stream(Warnings) 00308 << "*** ERROR *** MultiVecTraits::MvScale(alpha) " 00309 << "Supposedly zero vector has non-zero norm." << endl; 00310 return false; 00311 } 00312 } 00313 00314 MVT::MvRandom(*B); 00315 std::vector<ScalarType> zeros(numvecs,SCT::zero()); 00316 MVT::MvScale(*B,zeros); 00317 MVT::MvNorm(*B, norms); 00318 for (int i=0; i<numvecs; i++) { 00319 if ( norms[i] != zero_mag ) { 00320 om->stream(Warnings) 00321 << "*** ERROR *** MultiVecTraits::MvScale(alphas) " 00322 << "Supposedly zero vector has non-zero norm." << endl; 00323 return false; 00324 } 00325 } 00326 } 00327 00328 00329 /*********** MvInit() and MvNorm() *********************************** 00330 A vector of ones of dimension n should have norm sqrt(n) 00331 1) Init vectors to all ones 00332 2) Verify that norm is sqrt(n) 00333 3) Verify that norms aren't negative 00334 00335 Note: I'm not sure that we can expect this to hold in practice. 00336 Maybe something like abs(norm-sqrt(n)) < SCT::eps() ??? 00337 The sum of 1^2==1 should be n, but what about sqrt(n)? 00338 They may be using a different square root than ScalartTraits 00339 On my iBook G4 and on jeter, this test works. 00340 Right now, this has been demoted to a warning. 00341 *********************************************************************/ 00342 { 00343 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs); 00344 std::vector<MagType> norms(numvecs); 00345 00346 MVT::MvInit(*B,one); 00347 MVT::MvNorm(*B, norms); 00348 bool BadNormWarning = false; 00349 for (int i=0; i<numvecs; i++) { 00350 if ( norms[i] < zero_mag ) { 00351 om->stream(Warnings) 00352 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl 00353 << "Vector had negative norm." << endl; 00354 return false; 00355 } 00356 else if ( norms[i] != SCT::squareroot(MVT::GetVecLength(*B)) && !BadNormWarning ) { 00357 om->stream(Warnings) 00358 << endl 00359 << "Warning testing MultiVecTraits::MvInit()." << endl 00360 << "Ones vector should have norm sqrt(dim)." << endl 00361 << "norms[i]: " << norms[i] << "\tdim: " << MVT::GetVecLength(*B) << endl << endl; 00362 BadNormWarning = true; 00363 } 00364 } 00365 } 00366 00367 00368 /*********** MvInit() and MvNorm() *********************************** 00369 A vector of zeros of dimension n should have norm 0 00370 1) Verify that norms aren't negative 00371 2) Verify that norms are zero 00372 00373 We must know this works before the next tests. 00374 *********************************************************************/ 00375 { 00376 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs); 00377 std::vector<MagType> norms(numvecs); 00378 MVT::MvInit(*B, zero_mag); 00379 MVT::MvNorm(*B, norms); 00380 for (int i=0; i<numvecs; i++) { 00381 if ( norms[i] < zero_mag ) { 00382 om->stream(Warnings) 00383 << "*** ERROR *** MultiVecTraits::MvInit()." << endl 00384 << "Vector had negative norm." << endl; 00385 return false; 00386 } 00387 else if ( norms[i] != zero_mag ) { 00388 om->stream(Warnings) 00389 << "*** ERROR *** MultiVecTraits::MvInit()." << endl 00390 << "Zero vector should have norm zero." << endl; 00391 return false; 00392 } 00393 } 00394 } 00395 00396 00397 /*********** CloneCopy(MV,vector<int>) and MvNorm ******************** 00398 1) Check quantity/length of vectors 00399 2) Check vector norms for agreement 00400 3) Zero out B and make sure that C norms are not affected 00401 *********************************************************************/ 00402 { 00403 Teuchos::RCP<MV> B, C; 00404 std::vector<MagType> norms(numvecs), norms2(ind.size()); 00405 00406 B = MVT::Clone(*A,numvecs); 00407 MVT::MvRandom(*B); 00408 MVT::MvNorm(*B, norms); 00409 C = MVT::CloneCopy(*B,ind); 00410 MVT::MvNorm(*C, norms2); 00411 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) { 00412 om->stream(Warnings) 00413 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl 00414 << "Wrong number of vectors." << endl; 00415 return false; 00416 } 00417 if ( MVT::GetVecLength(*C) != MVT::GetVecLength(*B) ) { 00418 om->stream(Warnings) 00419 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl 00420 << "Vector lengths don't match." << endl; 00421 return false; 00422 } 00423 for (int i=0; i<numvecs_2; i++) { 00424 if ( norms2[i] != norms[ind[i]] ) { 00425 om->stream(Warnings) 00426 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl 00427 << "Copied vectors do not agree:" 00428 << norms2[i] << " != " << norms[ind[i]] << endl; 00429 return false; 00430 } 00431 } 00432 MVT::MvInit(*B,zero); 00433 MVT::MvNorm(*C, norms2); 00434 for (int i=0; i<numvecs_2; i++) { 00435 if ( norms2[i] != norms2[i] ) { 00436 om->stream(Warnings) 00437 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl 00438 << "Copied vectors were not independent." << endl; 00439 return false; 00440 } 00441 } 00442 } 00443 00444 00445 /*********** CloneCopy(MV) and MvNorm ******************************** 00446 1) Check quantity 00447 2) Check value of norms 00448 3) Zero out B and make sure that C is still okay 00449 *********************************************************************/ 00450 { 00451 Teuchos::RCP<MV> B, C; 00452 std::vector<MagType> norms(numvecs), norms2(numvecs); 00453 00454 B = MVT::Clone(*A,numvecs); 00455 MVT::MvRandom(*B); 00456 MVT::MvNorm(*B, norms); 00457 C = MVT::CloneCopy(*B); 00458 MVT::MvNorm(*C, norms2); 00459 if ( MVT::GetNumberVecs(*C) != numvecs ) { 00460 om->stream(Warnings) 00461 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl 00462 << "Wrong number of vectors." << endl; 00463 return false; 00464 } 00465 for (int i=0; i<numvecs; i++) { 00466 if ( norms2[i] != norms[i] ) { 00467 om->stream(Warnings) 00468 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl 00469 << "Copied vectors do not agree." << endl; 00470 return false; 00471 } 00472 } 00473 MVT::MvInit(*B,zero); 00474 MVT::MvNorm(*C, norms); 00475 for (int i=0; i<numvecs; i++) { 00476 if ( norms2[i] != norms[i] ) { 00477 om->stream(Warnings) 00478 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl 00479 << "Copied vectors were not independent." << endl; 00480 return false; 00481 } 00482 } 00483 } 00484 00485 00486 /*********** CloneView(MV,vector<int>) and MvNorm ******************** 00487 Check that we have a view of the selected vectors 00488 1) Check quantity 00489 2) Check value of norms 00490 3) Zero out B and make sure that C is zero as well 00491 *********************************************************************/ 00492 { 00493 Teuchos::RCP<MV> B, C; 00494 std::vector<MagType> norms(numvecs), norms2(ind.size()); 00495 00496 B = MVT::Clone(*A,numvecs); 00497 MVT::MvRandom(*B); 00498 MVT::MvNorm(*B, norms); 00499 C = MVT::CloneViewNonConst(*B,ind); 00500 MVT::MvNorm(*C, norms2); 00501 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) { 00502 om->stream(Warnings) 00503 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl 00504 << "Wrong number of vectors." << endl; 00505 return false; 00506 } 00507 for (int i=0; i<numvecs_2; i++) { 00508 if ( norms2[i] != norms[ind[i]] ) { 00509 om->stream(Warnings) 00510 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl 00511 << "Viewed vectors do not agree." << endl; 00512 return false; 00513 } 00514 } 00515 } 00516 00517 00518 /*********** const CloneView(MV,vector<int>) and MvNorm() ************ 00519 Check that we have a view of the selected vectors. 00520 1) Check quantity 00521 2) Check value of norms for agreement 00522 3) Zero out B and make sure that C is zerod as well 00523 *********************************************************************/ 00524 { 00525 Teuchos::RCP<MV> B; 00526 Teuchos::RCP<const MV> constB, C; 00527 std::vector<MagType> normsB(numvecs), normsC(ind.size()); 00528 std::vector<int> allind(numvecs); 00529 for (int i=0; i<numvecs; i++) { 00530 allind[i] = i; 00531 } 00532 00533 B = MVT::Clone(*A,numvecs); 00534 MVT::MvRandom( *B ); 00535 MVT::MvNorm(*B, normsB); 00536 C = MVT::CloneView(*B,ind); 00537 MVT::MvNorm(*C, normsC); 00538 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) { 00539 om->stream(Warnings) 00540 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl 00541 << "Wrong number of vectors." << endl; 00542 return false; 00543 } 00544 for (int i=0; i<numvecs_2; i++) { 00545 if ( normsC[i] != normsB[ind[i]] ) { 00546 om->stream(Warnings) 00547 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl 00548 << "Viewed vectors do not agree." << endl; 00549 return false; 00550 } 00551 } 00552 } 00553 00554 00555 /*********** SetBlock() and MvNorm() ********************************* 00556 SetBlock() will copy the vectors from C into B 00557 1) Verify that the specified vectors were copied 00558 2) Verify that the other vectors were not modified 00559 3) Verify that C was not modified 00560 4) Change C and then check B to make sure it was not modified 00561 00562 Use a different index set than has been used so far (distinct entries). 00563 This is because duplicate entries will cause the vector to be 00564 overwritten, making it more difficult to test. 00565 *********************************************************************/ 00566 { 00567 Teuchos::RCP<MV> B, C; 00568 std::vector<MagType> normsB1(numvecs), normsB2(numvecs), 00569 normsC1(numvecs_2), normsC2(numvecs_2); 00570 00571 B = MVT::Clone(*A,numvecs); 00572 C = MVT::Clone(*A,numvecs_2); 00573 // Just do every other one, interleaving the vectors of C into B 00574 ind.resize(numvecs_2); 00575 for (int i=0; i<numvecs_2; i++) { 00576 ind[i] = 2*i; 00577 } 00578 MVT::MvRandom(*B); 00579 MVT::MvRandom(*C); 00580 00581 MVT::MvNorm(*B,normsB1); 00582 MVT::MvNorm(*C,normsC1); 00583 MVT::SetBlock(*C,ind,*B); 00584 MVT::MvNorm(*B,normsB2); 00585 MVT::MvNorm(*C,normsC2); 00586 00587 // check that C was not changed by SetBlock 00588 for (int i=0; i<numvecs_2; i++) { 00589 if ( normsC1[i] != normsC2[i] ) { 00590 om->stream(Warnings) 00591 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl 00592 << "Operation modified source vectors." << endl; 00593 return false; 00594 } 00595 } 00596 // check that the correct vectors of B were modified 00597 // and the others were not 00598 for (int i=0; i<numvecs; i++) { 00599 if (i % 2 == 0) { 00600 // should be a vector from C 00601 if ( normsB2[i] != normsC1[i/2] ) { 00602 om->stream(Warnings) 00603 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl 00604 << "Copied vectors do not agree." << endl; 00605 return false; 00606 } 00607 } 00608 else { 00609 // should be an original vector 00610 if ( normsB1[i] != normsB2[i] ) { 00611 om->stream(Warnings) 00612 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl 00613 << "Incorrect vectors were modified." << endl; 00614 return false; 00615 } 00616 } 00617 } 00618 MVT::MvInit(*C,zero); 00619 MVT::MvNorm(*B,normsB1); 00620 // verify that we copied and didn't reference 00621 for (int i=0; i<numvecs; i++) { 00622 if ( normsB1[i] != normsB2[i] ) { 00623 om->stream(Warnings) 00624 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl 00625 << "Copied vectors were not independent." << endl; 00626 return false; 00627 } 00628 } 00629 } 00630 00631 00632 /*********** MvTransMv() ********************************************* 00633 Performs C = alpha * A^H * B, where 00634 alpha is type ScalarType 00635 A,B are type MV with p and q vectors, respectively 00636 C is a SerialDenseMatrix<int,ScalarType> ALREADY sized to p by q 00637 00638 Verify: 00639 1) C is not resized by the routine 00640 3) Check that zero*(A^H B) == zero 00641 3) Check inner product inequality: 00642 [ |a1|*|b1| ... |ap|*|b1| ] 00643 [a1 ... ap]^H [b1 ... bq] <= [ ... |ai|*|bj| ... ] 00644 [ |ap|*|b1| ... |ap|*|bq| ] 00645 4) Zero B and check that C is zero 00646 5) Zero A and check that C is zero 00647 00648 Note: Should we really require that C is correctly sized already? 00649 Epetra does (and crashes if it isn't.) 00650 *********************************************************************/ 00651 { 00652 const int p = 7; 00653 const int q = 9; 00654 Teuchos::RCP<MV> B, C; 00655 std::vector<MagType> normsB(p), normsC(q); 00656 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 00657 00658 B = MVT::Clone(*A,p); 00659 C = MVT::Clone(*A,q); 00660 00661 // randomize the multivectors 00662 MVT::MvRandom(*B); 00663 MVT::MvNorm(*B,normsB); 00664 MVT::MvRandom(*C); 00665 MVT::MvNorm(*C,normsC); 00666 00667 // perform SDM = zero() * B^H * C 00668 MVT::MvTransMv( zero, *B, *C, SDM ); 00669 00670 // check the sizes: not allowed to have shrunk 00671 if ( SDM.numRows() != p || SDM.numCols() != q ) { 00672 om->stream(Warnings) 00673 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00674 << "Routine resized SerialDenseMatrix." << endl; 00675 return false; 00676 } 00677 00678 // check that zero**A^H*B == zero 00679 if ( SDM.normOne() != zero ) { 00680 om->stream(Warnings) 00681 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00682 << "Scalar argument processed incorrectly." << endl; 00683 return false; 00684 } 00685 00686 // perform SDM = one * B^H * C 00687 MVT::MvTransMv( one, *B, *C, SDM ); 00688 00689 // check the norms: a^H b = |a| |b| cos(theta) <= |a| |b| 00690 // with equality only when a and b are colinear 00691 for (int i=0; i<p; i++) { 00692 for (int j=0; j<q; j++) { 00693 if ( SCT::magnitude(SDM(i,j)) 00694 > SCT::magnitude(normsB[i]*normsC[j]) ) { 00695 om->stream(Warnings) 00696 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00697 << "Triangle inequality did not hold: " 00698 << SCT::magnitude(SDM(i,j)) 00699 << " > " 00700 << SCT::magnitude(normsB[i]*normsC[j]) 00701 << endl; 00702 return false; 00703 } 00704 } 00705 } 00706 MVT::MvInit(*C); 00707 MVT::MvRandom(*B); 00708 MVT::MvTransMv( one, *B, *C, SDM ); 00709 for (int i=0; i<p; i++) { 00710 for (int j=0; j<q; j++) { 00711 if ( SDM(i,j) != zero ) { 00712 om->stream(Warnings) 00713 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00714 << "Inner products not zero for C==0." << endl; 00715 return false; 00716 } 00717 } 00718 } 00719 MVT::MvInit(*B); 00720 MVT::MvRandom(*C); 00721 MVT::MvTransMv( one, *B, *C, SDM ); 00722 for (int i=0; i<p; i++) { 00723 for (int j=0; j<q; j++) { 00724 if ( SDM(i,j) != zero ) { 00725 om->stream(Warnings) 00726 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00727 << "Inner products not zero for B==0." << endl; 00728 return false; 00729 } 00730 } 00731 } 00732 } 00733 00734 00735 /*********** MvDot() ************************************************* 00736 Verify: 00737 1) Results vector not resized 00738 2) Inner product inequalities are satisfied 00739 3) Zero vectors give zero inner products 00740 *********************************************************************/ 00741 { 00742 const int p = 7; 00743 Teuchos::RCP<MV> B, C; 00744 std::vector<ScalarType> iprods(p); 00745 std::vector<MagType> normsB(p), normsC(p); 00746 00747 B = MVT::Clone(*A,p); 00748 C = MVT::Clone(*A,p); 00749 00750 MVT::MvRandom(*B); 00751 MVT::MvRandom(*C); 00752 MVT::MvNorm(*B,normsB); 00753 MVT::MvNorm(*C,normsC); 00754 MVT::MvDot( *B, *C, iprods ); 00755 if ( (int)iprods.size() != p ) { 00756 om->stream(Warnings) 00757 << "*** ERROR *** MultiVecTraits::MvDot." << endl 00758 << "Routine resized results vector." << endl; 00759 return false; 00760 } 00761 for (int i=0; i<p; i++) { 00762 if ( SCT::magnitude(iprods[i]) 00763 > SCT::magnitude(normsB[i]*normsC[i]) ) { 00764 om->stream(Warnings) 00765 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00766 << "Inner products not valid." << endl; 00767 return false; 00768 } 00769 } 00770 MVT::MvInit(*B); 00771 MVT::MvRandom(*C); 00772 MVT::MvDot( *B, *C, iprods ); 00773 for (int i=0; i<p; i++) { 00774 if ( iprods[i] != zero ) { 00775 om->stream(Warnings) 00776 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00777 << "Inner products not zero for B==0." << endl; 00778 return false; 00779 } 00780 } 00781 MVT::MvInit(*C); 00782 MVT::MvRandom(*B); 00783 MVT::MvDot( *B, *C, iprods ); 00784 for (int i=0; i<p; i++) { 00785 if ( iprods[i] != zero ) { 00786 om->stream(Warnings) 00787 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00788 << "Inner products not zero for C==0." << endl; 00789 return false; 00790 } 00791 } 00792 } 00793 00794 00795 /*********** MvAddMv() *********************************************** 00796 D = alpha*B + beta*C 00797 1) Use alpha==0,beta==1 and check that D == C 00798 2) Use alpha==1,beta==0 and check that D == B 00799 3) Use D==0 and D!=0 and check that result is the same 00800 4) Check that input arguments are not modified 00801 *********************************************************************/ 00802 { 00803 const int p = 7; 00804 Teuchos::RCP<MV> B, C, D; 00805 std::vector<MagType> normsB1(p), normsB2(p), 00806 normsC1(p), normsC2(p), 00807 normsD1(p), normsD2(p); 00808 ScalarType alpha = SCT::random(), 00809 beta = SCT::random(); 00810 00811 B = MVT::Clone(*A,p); 00812 C = MVT::Clone(*A,p); 00813 D = MVT::Clone(*A,p); 00814 00815 MVT::MvRandom(*B); 00816 MVT::MvRandom(*C); 00817 MVT::MvNorm(*B,normsB1); 00818 MVT::MvNorm(*C,normsC1); 00819 00820 // check that 0*B+1*C == C 00821 MVT::MvAddMv(zero,*B,one,*C,*D); 00822 MVT::MvNorm(*B,normsB2); 00823 MVT::MvNorm(*C,normsC2); 00824 MVT::MvNorm(*D,normsD1); 00825 for (int i=0; i<p; i++) { 00826 if ( normsB1[i] != normsB2[i] ) { 00827 om->stream(Warnings) 00828 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00829 << "Input arguments were modified." << endl; 00830 return false; 00831 } 00832 else if ( normsC1[i] != normsC2[i] ) { 00833 om->stream(Warnings) 00834 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00835 << "Input arguments were modified." << endl; 00836 return false; 00837 } 00838 else if ( normsC1[i] != normsD1[i] ) { 00839 om->stream(Warnings) 00840 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00841 << "Assignment did not work." << endl; 00842 return false; 00843 } 00844 } 00845 00846 // check that 1*B+0*C == B 00847 MVT::MvAddMv(one,*B,zero,*C,*D); 00848 MVT::MvNorm(*B,normsB2); 00849 MVT::MvNorm(*C,normsC2); 00850 MVT::MvNorm(*D,normsD1); 00851 for (int i=0; i<p; i++) { 00852 if ( normsB1[i] != normsB2[i] ) { 00853 om->stream(Warnings) 00854 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00855 << "Input arguments were modified." << endl; 00856 return false; 00857 } 00858 else if ( normsC1[i] != normsC2[i] ) { 00859 om->stream(Warnings) 00860 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00861 << "Input arguments were modified." << endl; 00862 return false; 00863 } 00864 else if ( normsB1[i] != normsD1[i] ) { 00865 om->stream(Warnings) 00866 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00867 << "Assignment did not work." << endl; 00868 return false; 00869 } 00870 } 00871 00872 // check that alpha*B+beta*C -> D is invariant under initial D 00873 // first, try random D 00874 MVT::MvRandom(*D); 00875 MVT::MvAddMv(alpha,*B,beta,*C,*D); 00876 MVT::MvNorm(*B,normsB2); 00877 MVT::MvNorm(*C,normsC2); 00878 MVT::MvNorm(*D,normsD1); 00879 // check that input args are not modified 00880 for (int i=0; i<p; i++) { 00881 if ( normsB1[i] != normsB2[i] ) { 00882 om->stream(Warnings) 00883 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00884 << "Input arguments were modified." << endl; 00885 return false; 00886 } 00887 else if ( normsC1[i] != normsC2[i] ) { 00888 om->stream(Warnings) 00889 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00890 << "Input arguments were modified." << endl; 00891 return false; 00892 } 00893 } 00894 // next, try zero D 00895 MVT::MvInit(*D); 00896 MVT::MvAddMv(alpha,*B,beta,*C,*D); 00897 MVT::MvNorm(*B,normsB2); 00898 MVT::MvNorm(*C,normsC2); 00899 MVT::MvNorm(*D,normsD2); 00900 // check that input args are not modified and that D is the same 00901 // as the above test 00902 for (int i=0; i<p; i++) { 00903 if ( normsB1[i] != normsB2[i] ) { 00904 om->stream(Warnings) 00905 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00906 << "Input arguments were modified." << endl; 00907 return false; 00908 } 00909 else if ( normsC1[i] != normsC2[i] ) { 00910 om->stream(Warnings) 00911 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00912 << "Input arguments were modified." << endl; 00913 return false; 00914 } 00915 else if ( normsD1[i] != normsD2[i] ) { 00916 om->stream(Warnings) 00917 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00918 << "Results varies depending on initial state of dest vectors." << endl; 00919 return false; 00920 } 00921 } 00922 } 00923 00924 /*********** MvAddMv() *********************************************** 00925 Similar to above, but where B or C are potentially the same 00926 object as D. This case is commonly used, for example, to affect 00927 A <- alpha*A 00928 via 00929 MvAddMv(alpha,A,zero,A,A) 00930 ** OR ** 00931 MvAddMv(zero,A,alpha,A,A) 00932 00933 The result is that the operation has to be "atomic". That is, 00934 B and C are no longer reliable after D is modified, so that 00935 the assignment to D must be the last thing to occur. 00936 00937 D = alpha*B + beta*C 00938 00939 1) Use alpha==0,beta==1 and check that D == C 00940 2) Use alpha==1,beta==0 and check that D == B 00941 *********************************************************************/ 00942 { 00943 const int p = 7; 00944 Teuchos::RCP<MV> B, D; 00945 Teuchos::RCP<const MV> C; 00946 std::vector<MagType> normsB(p), 00947 normsD(p); 00948 std::vector<int> lclindex(p); 00949 for (int i=0; i<p; i++) lclindex[i] = i; 00950 00951 B = MVT::Clone(*A,p); 00952 C = MVT::CloneView(*B,lclindex); 00953 D = MVT::CloneViewNonConst(*B,lclindex); 00954 00955 MVT::MvRandom(*B); 00956 MVT::MvNorm(*B,normsB); 00957 00958 // check that 0*B+1*C == C 00959 MVT::MvAddMv(zero,*B,one,*C,*D); 00960 MVT::MvNorm(*D,normsD); 00961 for (int i=0; i<p; i++) { 00962 if ( normsB[i] != normsD[i] ) { 00963 om->stream(Warnings) 00964 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl 00965 << "Assignment did not work." << endl; 00966 return false; 00967 } 00968 } 00969 00970 // check that 1*B+0*C == B 00971 MVT::MvAddMv(one,*B,zero,*C,*D); 00972 MVT::MvNorm(*D,normsD); 00973 for (int i=0; i<p; i++) { 00974 if ( normsB[i] != normsD[i] ) { 00975 om->stream(Warnings) 00976 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl 00977 << "Assignment did not work." << endl; 00978 return false; 00979 } 00980 } 00981 00982 } 00983 00984 00985 /*********** MvTimesMatAddMv() 7 by 5 ******************************** 00986 C = alpha*B*SDM + beta*C 00987 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged 00988 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero 00989 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B 00990 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged 00991 5) Test with non-square matrices 00992 6) Always check that input arguments are not modified 00993 *********************************************************************/ 00994 { 00995 const int p = 7, q = 5; 00996 Teuchos::RCP<MV> B, C; 00997 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 00998 std::vector<MagType> normsC1(q), normsC2(q), 00999 normsB1(p), normsB2(p); 01000 01001 B = MVT::Clone(*A,p); 01002 C = MVT::Clone(*A,q); 01003 01004 // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged 01005 MVT::MvRandom(*B); 01006 MVT::MvRandom(*C); 01007 MVT::MvNorm(*B,normsB1); 01008 MVT::MvNorm(*C,normsC1); 01009 SDM.random(); 01010 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C); 01011 MVT::MvNorm(*B,normsB2); 01012 MVT::MvNorm(*C,normsC2); 01013 for (int i=0; i<p; i++) { 01014 if ( normsB1[i] != normsB2[i] ) { 01015 om->stream(Warnings) 01016 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01017 << "Input vectors were modified." << endl; 01018 return false; 01019 } 01020 } 01021 for (int i=0; i<q; i++) { 01022 if ( normsC1[i] != normsC2[i] ) { 01023 om->stream(Warnings) 01024 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01025 << "Arithmetic test 1 failed." << endl; 01026 return false; 01027 } 01028 } 01029 01030 // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero 01031 MVT::MvRandom(*B); 01032 MVT::MvRandom(*C); 01033 MVT::MvNorm(*B,normsB1); 01034 MVT::MvNorm(*C,normsC1); 01035 SDM.random(); 01036 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C); 01037 MVT::MvNorm(*B,normsB2); 01038 MVT::MvNorm(*C,normsC2); 01039 for (int i=0; i<p; i++) { 01040 if ( normsB1[i] != normsB2[i] ) { 01041 om->stream(Warnings) 01042 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01043 << "Input vectors were modified." << endl; 01044 return false; 01045 } 01046 } 01047 for (int i=0; i<q; i++) { 01048 if ( normsC2[i] != zero ) { 01049 om->stream(Warnings) 01050 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01051 << "Arithmetic test 2 failed: " 01052 << normsC2[i] 01053 << " != " 01054 << zero 01055 << endl; 01056 return false; 01057 } 01058 } 01059 01060 // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B 01061 // |0| 01062 MVT::MvRandom(*B); 01063 MVT::MvRandom(*C); 01064 MVT::MvNorm(*B,normsB1); 01065 MVT::MvNorm(*C,normsC1); 01066 SDM.scale(zero); 01067 for (int i=0; i<q; i++) { 01068 SDM(i,i) = one; 01069 } 01070 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C); 01071 MVT::MvNorm(*B,normsB2); 01072 MVT::MvNorm(*C,normsC2); 01073 for (int i=0; i<p; i++) { 01074 if ( normsB1[i] != normsB2[i] ) { 01075 om->stream(Warnings) 01076 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01077 << "Input vectors were modified." << endl; 01078 return false; 01079 } 01080 } 01081 for (int i=0; i<q; i++) { 01082 if ( normsB1[i] != normsC2[i] ) { 01083 om->stream(Warnings) 01084 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01085 << "Arithmetic test 3 failed: " 01086 << normsB1[i] 01087 << " != " 01088 << normsC2[i] 01089 << endl; 01090 return false; 01091 } 01092 } 01093 01094 // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged 01095 MVT::MvRandom(*B); 01096 MVT::MvRandom(*C); 01097 MVT::MvNorm(*B,normsB1); 01098 MVT::MvNorm(*C,normsC1); 01099 SDM.scale(zero); 01100 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C); 01101 MVT::MvNorm(*B,normsB2); 01102 MVT::MvNorm(*C,normsC2); 01103 for (int i=0; i<p; i++) { 01104 if ( normsB1[i] != normsB2[i] ) { 01105 om->stream(Warnings) 01106 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01107 << "Input vectors were modified." << endl; 01108 return false; 01109 } 01110 } 01111 for (int i=0; i<q; i++) { 01112 if ( normsC1[i] != normsC2[i] ) { 01113 om->stream(Warnings) 01114 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01115 << "Arithmetic test 4 failed." << endl; 01116 return false; 01117 } 01118 } 01119 } 01120 01121 /*********** MvTimesMatAddMv() 5 by 7 ******************************** 01122 C = alpha*B*SDM + beta*C 01123 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged 01124 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero 01125 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B 01126 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged 01127 5) Test with non-square matrices 01128 6) Always check that input arguments are not modified 01129 *********************************************************************/ 01130 { 01131 const int p = 5, q = 7; 01132 Teuchos::RCP<MV> B, C; 01133 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 01134 std::vector<MagType> normsC1(q), normsC2(q), 01135 normsB1(p), normsB2(p); 01136 01137 B = MVT::Clone(*A,p); 01138 C = MVT::Clone(*A,q); 01139 01140 // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged 01141 MVT::MvRandom(*B); 01142 MVT::MvRandom(*C); 01143 MVT::MvNorm(*B,normsB1); 01144 MVT::MvNorm(*C,normsC1); 01145 SDM.random(); 01146 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C); 01147 MVT::MvNorm(*B,normsB2); 01148 MVT::MvNorm(*C,normsC2); 01149 for (int i=0; i<p; i++) { 01150 if ( normsB1[i] != normsB2[i] ) { 01151 om->stream(Warnings) 01152 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01153 << "Input vectors were modified." << endl; 01154 return false; 01155 } 01156 } 01157 for (int i=0; i<q; i++) { 01158 if ( normsC1[i] != normsC2[i] ) { 01159 om->stream(Warnings) 01160 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01161 << "Arithmetic test 5 failed." << endl; 01162 return false; 01163 } 01164 } 01165 01166 // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero 01167 MVT::MvRandom(*B); 01168 MVT::MvRandom(*C); 01169 MVT::MvNorm(*B,normsB1); 01170 MVT::MvNorm(*C,normsC1); 01171 SDM.random(); 01172 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C); 01173 MVT::MvNorm(*B,normsB2); 01174 MVT::MvNorm(*C,normsC2); 01175 for (int i=0; i<p; i++) { 01176 if ( normsB1[i] != normsB2[i] ) { 01177 om->stream(Warnings) 01178 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01179 << "Input vectors were modified." << endl; 01180 return false; 01181 } 01182 } 01183 for (int i=0; i<q; i++) { 01184 if ( normsC2[i] != zero ) { 01185 om->stream(Warnings) 01186 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01187 << "Arithmetic test 6 failed: " 01188 << normsC2[i] 01189 << " != " 01190 << zero 01191 << endl; 01192 return false; 01193 } 01194 } 01195 01196 // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B 01197 MVT::MvRandom(*B); 01198 MVT::MvRandom(*C); 01199 MVT::MvNorm(*B,normsB1); 01200 MVT::MvNorm(*C,normsC1); 01201 SDM.scale(zero); 01202 for (int i=0; i<p; i++) { 01203 SDM(i,i) = one; 01204 } 01205 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C); 01206 MVT::MvNorm(*B,normsB2); 01207 MVT::MvNorm(*C,normsC2); 01208 for (int i=0; i<p; i++) { 01209 if ( normsB1[i] != normsB2[i] ) { 01210 om->stream(Warnings) 01211 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01212 << "Input vectors were modified." << endl; 01213 return false; 01214 } 01215 } 01216 for (int i=0; i<p; i++) { 01217 if ( normsB1[i] != normsC2[i] ) { 01218 om->stream(Warnings) 01219 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01220 << "Arithmetic test 7 failed." << endl; 01221 return false; 01222 } 01223 } 01224 for (int i=p; i<q; i++) { 01225 if ( normsC2[i] != zero ) { 01226 om->stream(Warnings) 01227 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01228 << "Arithmetic test 7 failed." << endl; 01229 return false; 01230 } 01231 } 01232 01233 // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged 01234 MVT::MvRandom(*B); 01235 MVT::MvRandom(*C); 01236 MVT::MvNorm(*B,normsB1); 01237 MVT::MvNorm(*C,normsC1); 01238 SDM.scale(zero); 01239 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C); 01240 MVT::MvNorm(*B,normsB2); 01241 MVT::MvNorm(*C,normsC2); 01242 for (int i=0; i<p; i++) { 01243 if ( normsB1[i] != normsB2[i] ) { 01244 om->stream(Warnings) 01245 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01246 << "Input vectors were modified." << endl; 01247 return false; 01248 } 01249 } 01250 for (int i=0; i<q; i++) { 01251 if ( normsC1[i] != normsC2[i] ) { 01252 om->stream(Warnings) 01253 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01254 << "Arithmetic test 8 failed." << endl; 01255 return false; 01256 } 01257 } 01258 } 01259 01260 return true; 01261 01262 } 01263 01264 01265 01271 template< class ScalarType, class MV, class OP> 01272 bool TestOperatorTraits( 01273 const Teuchos::RCP<OutputManager<ScalarType> > &om, 01274 const Teuchos::RCP<const MV> &A, 01275 const Teuchos::RCP<const OP> &M) { 01276 01277 using std::endl; 01278 01279 /* OPT Contract: 01280 Apply() 01281 MV: OP*zero == zero 01282 Warn if OP is not deterministic (OP*A != OP*A) 01283 Does not modify input arguments 01284 *********************************************************************/ 01285 01286 typedef MultiVecTraits<ScalarType, MV> MVT; 01287 typedef Teuchos::ScalarTraits<ScalarType> SCT; 01288 typedef OperatorTraits<ScalarType, MV, OP> OPT; 01289 typedef typename SCT::magnitudeType MagType; 01290 01291 const int numvecs = 10; 01292 01293 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs), 01294 C = MVT::Clone(*A,numvecs); 01295 01296 std::vector<MagType> normsB1(numvecs), normsB2(numvecs), 01297 normsC1(numvecs), normsC2(numvecs); 01298 bool NonDeterministicWarning; 01299 01300 /*********** Apply() ************************************************* 01301 Verify: 01302 1) OP*B == OP*B; OP is deterministic (just warn on this) 01303 2) OP*zero == 0 01304 3) OP*B doesn't modify B 01305 4) OP*B is invariant under initial state of destination vectors 01306 *********************************************************************/ 01307 MVT::MvInit(*B); 01308 MVT::MvRandom(*C); 01309 MVT::MvNorm(*B,normsB1); 01310 OPT::Apply(*M,*B,*C); 01311 MVT::MvNorm(*B,normsB2); 01312 MVT::MvNorm(*C,normsC2); 01313 for (int i=0; i<numvecs; i++) { 01314 if (normsB2[i] != normsB1[i]) { 01315 om->stream(Warnings) 01316 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl 01317 << "Apply() modified the input vectors." << endl; 01318 return false; 01319 } 01320 if (normsC2[i] != SCT::zero()) { 01321 om->stream(Warnings) 01322 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl 01323 << "Operator applied to zero did not return zero." << endl; 01324 return false; 01325 } 01326 } 01327 01328 // If we send in a random matrix, we should not get a zero return 01329 MVT::MvRandom(*B); 01330 MVT::MvNorm(*B,normsB1); 01331 OPT::Apply(*M,*B,*C); 01332 MVT::MvNorm(*B,normsB2); 01333 MVT::MvNorm(*C,normsC2); 01334 bool ZeroWarning = false; 01335 for (int i=0; i<numvecs; i++) { 01336 if (normsB2[i] != normsB1[i]) { 01337 om->stream(Warnings) 01338 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl 01339 << "Apply() modified the input vectors." << endl; 01340 return false; 01341 } 01342 if (normsC2[i] == SCT::zero() && ZeroWarning==false ) { 01343 om->stream(Warnings) 01344 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl 01345 << "Operator applied to random vectors returned zero." << endl; 01346 ZeroWarning = true; 01347 } 01348 } 01349 01350 // Apply operator with C init'd to zero 01351 MVT::MvRandom(*B); 01352 MVT::MvNorm(*B,normsB1); 01353 MVT::MvInit(*C); 01354 OPT::Apply(*M,*B,*C); 01355 MVT::MvNorm(*B,normsB2); 01356 MVT::MvNorm(*C,normsC1); 01357 for (int i=0; i<numvecs; i++) { 01358 if (normsB2[i] != normsB1[i]) { 01359 om->stream(Warnings) 01360 << "*** ERROR *** OperatorTraits::Apply() [3]" << endl 01361 << "Apply() modified the input vectors." << endl; 01362 return false; 01363 } 01364 } 01365 01366 // Apply operator with C init'd to random 01367 // Check that result is the same as before; warn if not. 01368 // This could be a result of a bug, or a stochastic 01369 // operator. We do not want to prejudice against a 01370 // stochastic operator. 01371 MVT::MvRandom(*C); 01372 OPT::Apply(*M,*B,*C); 01373 MVT::MvNorm(*B,normsB2); 01374 MVT::MvNorm(*C,normsC2); 01375 NonDeterministicWarning = false; 01376 for (int i=0; i<numvecs; i++) { 01377 if (normsB2[i] != normsB1[i]) { 01378 om->stream(Warnings) 01379 << "*** ERROR *** OperatorTraits::Apply() [4]" << endl 01380 << "Apply() modified the input vectors." << endl; 01381 return false; 01382 } 01383 if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) { 01384 om->stream(Warnings) 01385 << endl 01386 << "*** WARNING *** OperatorTraits::Apply() [4]" << endl 01387 << "Apply() returned two different results." << endl << endl; 01388 NonDeterministicWarning = true; 01389 } 01390 } 01391 01392 return true; 01393 01394 } 01395 01396 } 01397 01398 #endif
1.7.4
