|
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 B^H * C is zero 00646 5) Zero C and check that B^H * C is zero 00647 00648 Note 1: Test 4 is performed with a p x q Teuchos::SDM view of 00649 a (p+1) x (q+1) Teuchos::SDM that is initialized to ones. 00650 This ensures the user is correctly accessing and filling the SDM. 00651 00652 Note 2: Should we really require that C is correctly sized already? 00653 Epetra does (and crashes if it isn't.) 00654 *********************************************************************/ 00655 { 00656 const int p = 7; 00657 const int q = 9; 00658 Teuchos::RCP<MV> B, C; 00659 std::vector<MagType> normsB(p), normsC(q); 00660 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 00661 00662 B = MVT::Clone(*A,p); 00663 C = MVT::Clone(*A,q); 00664 00665 // randomize the multivectors 00666 MVT::MvRandom(*B); 00667 MVT::MvNorm(*B,normsB); 00668 MVT::MvRandom(*C); 00669 MVT::MvNorm(*C,normsC); 00670 00671 // perform SDM = zero() * B^H * C 00672 MVT::MvTransMv( zero, *B, *C, SDM ); 00673 00674 // check the sizes: not allowed to have shrunk 00675 if ( SDM.numRows() != p || SDM.numCols() != q ) { 00676 om->stream(Warnings) 00677 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00678 << "Routine resized SerialDenseMatrix." << endl; 00679 return false; 00680 } 00681 00682 // check that zero**A^H*B == zero 00683 if ( SDM.normOne() != zero ) { 00684 om->stream(Warnings) 00685 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00686 << "Scalar argument processed incorrectly." << endl; 00687 return false; 00688 } 00689 00690 // perform SDM = one * B^H * C 00691 MVT::MvTransMv( one, *B, *C, SDM ); 00692 00693 // check the norms: a^H b = |a| |b| cos(theta) <= |a| |b| 00694 // with equality only when a and b are colinear 00695 for (int i=0; i<p; i++) { 00696 for (int j=0; j<q; j++) { 00697 if ( SCT::magnitude(SDM(i,j)) 00698 > SCT::magnitude(normsB[i]*normsC[j]) ) { 00699 om->stream(Warnings) 00700 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00701 << "Triangle inequality did not hold: " 00702 << SCT::magnitude(SDM(i,j)) 00703 << " > " 00704 << SCT::magnitude(normsB[i]*normsC[j]) 00705 << endl; 00706 return false; 00707 } 00708 } 00709 } 00710 MVT::MvInit(*C); 00711 MVT::MvRandom(*B); 00712 MVT::MvTransMv( one, *B, *C, SDM ); 00713 for (int i=0; i<p; i++) { 00714 for (int j=0; j<q; j++) { 00715 if ( SDM(i,j) != zero ) { 00716 om->stream(Warnings) 00717 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00718 << "Inner products not zero for C==0." << endl; 00719 return false; 00720 } 00721 } 00722 } 00723 MVT::MvInit(*B); 00724 MVT::MvRandom(*C); 00725 MVT::MvTransMv( one, *B, *C, SDM ); 00726 for (int i=0; i<p; i++) { 00727 for (int j=0; j<q; j++) { 00728 if ( SDM(i,j) != zero ) { 00729 om->stream(Warnings) 00730 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00731 << "Inner products not zero for B==0." << endl; 00732 return false; 00733 } 00734 } 00735 } 00736 00737 // A larger SDM is filled with ones, initially, and a smaller 00738 // view is used for the MvTransMv method. If the smaller SDM 00739 // is not all zeroes, then the interface is improperly writing 00740 // to the matrix object. 00741 // Note: Since we didn't fail above, we know that the general 00742 // inner product works, but we are checking to see if it 00743 // works for a view too. This is common usage in Anasazi. 00744 Teuchos::SerialDenseMatrix<int, ScalarType> largeSDM(p+1,q+1); 00745 Teuchos::SerialDenseMatrix<int, ScalarType> SDM2(Teuchos::View, largeSDM, p, q); 00746 largeSDM.putScalar( one ); 00747 MVT::MvInit(*C); 00748 MVT::MvRandom(*B); 00749 MVT::MvTransMv( one, *B, *C, SDM2 ); 00750 for (int i=0; i<p; i++) { 00751 for (int j=0; j<q; j++) { 00752 if ( SDM2(i,j) != zero ) { 00753 om->stream(Warnings) 00754 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl 00755 << "Inner products not zero for C==0 when using a view into Teuchos::SerialDenseMatrix<>." << endl; 00756 return false; 00757 } 00758 } 00759 } 00760 } 00761 00762 00763 /*********** MvDot() ************************************************* 00764 Verify: 00765 1) Results vector not resized 00766 2) Inner product inequalities are satisfied 00767 3) Zero vectors give zero inner products 00768 *********************************************************************/ 00769 { 00770 const int p = 7; 00771 Teuchos::RCP<MV> B, C; 00772 std::vector<ScalarType> iprods(p); 00773 std::vector<MagType> normsB(p), normsC(p); 00774 00775 B = MVT::Clone(*A,p); 00776 C = MVT::Clone(*A,p); 00777 00778 MVT::MvRandom(*B); 00779 MVT::MvRandom(*C); 00780 MVT::MvNorm(*B,normsB); 00781 MVT::MvNorm(*C,normsC); 00782 MVT::MvDot( *B, *C, iprods ); 00783 if ( (int)iprods.size() != p ) { 00784 om->stream(Warnings) 00785 << "*** ERROR *** MultiVecTraits::MvDot." << endl 00786 << "Routine resized results vector." << endl; 00787 return false; 00788 } 00789 for (int i=0; i<p; i++) { 00790 if ( SCT::magnitude(iprods[i]) 00791 > SCT::magnitude(normsB[i]*normsC[i]) ) { 00792 om->stream(Warnings) 00793 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00794 << "Inner products not valid." << endl; 00795 return false; 00796 } 00797 } 00798 MVT::MvInit(*B); 00799 MVT::MvRandom(*C); 00800 MVT::MvDot( *B, *C, iprods ); 00801 for (int i=0; i<p; i++) { 00802 if ( iprods[i] != zero ) { 00803 om->stream(Warnings) 00804 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00805 << "Inner products not zero for B==0." << endl; 00806 return false; 00807 } 00808 } 00809 MVT::MvInit(*C); 00810 MVT::MvRandom(*B); 00811 MVT::MvDot( *B, *C, iprods ); 00812 for (int i=0; i<p; i++) { 00813 if ( iprods[i] != zero ) { 00814 om->stream(Warnings) 00815 << "*** ERROR *** MultiVecTraits::MvDot()." << endl 00816 << "Inner products not zero for C==0." << endl; 00817 return false; 00818 } 00819 } 00820 } 00821 00822 00823 /*********** MvAddMv() *********************************************** 00824 D = alpha*B + beta*C 00825 1) Use alpha==0,beta==1 and check that D == C 00826 2) Use alpha==1,beta==0 and check that D == B 00827 3) Use D==0 and D!=0 and check that result is the same 00828 4) Check that input arguments are not modified 00829 *********************************************************************/ 00830 { 00831 const int p = 7; 00832 Teuchos::RCP<MV> B, C, D; 00833 std::vector<MagType> normsB1(p), normsB2(p), 00834 normsC1(p), normsC2(p), 00835 normsD1(p), normsD2(p); 00836 ScalarType alpha = SCT::random(), 00837 beta = SCT::random(); 00838 00839 B = MVT::Clone(*A,p); 00840 C = MVT::Clone(*A,p); 00841 D = MVT::Clone(*A,p); 00842 00843 MVT::MvRandom(*B); 00844 MVT::MvRandom(*C); 00845 MVT::MvNorm(*B,normsB1); 00846 MVT::MvNorm(*C,normsC1); 00847 00848 // check that 0*B+1*C == C 00849 MVT::MvAddMv(zero,*B,one,*C,*D); 00850 MVT::MvNorm(*B,normsB2); 00851 MVT::MvNorm(*C,normsC2); 00852 MVT::MvNorm(*D,normsD1); 00853 for (int i=0; i<p; i++) { 00854 if ( normsB1[i] != normsB2[i] ) { 00855 om->stream(Warnings) 00856 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00857 << "Input arguments were modified." << endl; 00858 return false; 00859 } 00860 else if ( normsC1[i] != normsC2[i] ) { 00861 om->stream(Warnings) 00862 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00863 << "Input arguments were modified." << endl; 00864 return false; 00865 } 00866 else if ( normsC1[i] != normsD1[i] ) { 00867 om->stream(Warnings) 00868 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00869 << "Assignment did not work." << endl; 00870 return false; 00871 } 00872 } 00873 00874 // check that 1*B+0*C == B 00875 MVT::MvAddMv(one,*B,zero,*C,*D); 00876 MVT::MvNorm(*B,normsB2); 00877 MVT::MvNorm(*C,normsC2); 00878 MVT::MvNorm(*D,normsD1); 00879 for (int i=0; i<p; i++) { 00880 if ( normsB1[i] != normsB2[i] ) { 00881 om->stream(Warnings) 00882 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00883 << "Input arguments were modified." << endl; 00884 return false; 00885 } 00886 else if ( normsC1[i] != normsC2[i] ) { 00887 om->stream(Warnings) 00888 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00889 << "Input arguments were modified." << endl; 00890 return false; 00891 } 00892 else if ( normsB1[i] != normsD1[i] ) { 00893 om->stream(Warnings) 00894 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00895 << "Assignment did not work." << endl; 00896 return false; 00897 } 00898 } 00899 00900 // check that alpha*B+beta*C -> D is invariant under initial D 00901 // first, try random D 00902 MVT::MvRandom(*D); 00903 MVT::MvAddMv(alpha,*B,beta,*C,*D); 00904 MVT::MvNorm(*B,normsB2); 00905 MVT::MvNorm(*C,normsC2); 00906 MVT::MvNorm(*D,normsD1); 00907 // check that input args are not modified 00908 for (int i=0; i<p; i++) { 00909 if ( normsB1[i] != normsB2[i] ) { 00910 om->stream(Warnings) 00911 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00912 << "Input arguments were modified." << endl; 00913 return false; 00914 } 00915 else if ( normsC1[i] != normsC2[i] ) { 00916 om->stream(Warnings) 00917 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00918 << "Input arguments were modified." << endl; 00919 return false; 00920 } 00921 } 00922 // next, try zero D 00923 MVT::MvInit(*D); 00924 MVT::MvAddMv(alpha,*B,beta,*C,*D); 00925 MVT::MvNorm(*B,normsB2); 00926 MVT::MvNorm(*C,normsC2); 00927 MVT::MvNorm(*D,normsD2); 00928 // check that input args are not modified and that D is the same 00929 // as the above test 00930 for (int i=0; i<p; i++) { 00931 if ( normsB1[i] != normsB2[i] ) { 00932 om->stream(Warnings) 00933 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00934 << "Input arguments were modified." << endl; 00935 return false; 00936 } 00937 else if ( normsC1[i] != normsC2[i] ) { 00938 om->stream(Warnings) 00939 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00940 << "Input arguments were modified." << endl; 00941 return false; 00942 } 00943 else if ( normsD1[i] != normsD2[i] ) { 00944 om->stream(Warnings) 00945 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl 00946 << "Results varies depending on initial state of dest vectors." << endl; 00947 return false; 00948 } 00949 } 00950 } 00951 00952 /*********** MvAddMv() *********************************************** 00953 Similar to above, but where B or C are potentially the same 00954 object as D. This case is commonly used, for example, to affect 00955 A <- alpha*A 00956 via 00957 MvAddMv(alpha,A,zero,A,A) 00958 ** OR ** 00959 MvAddMv(zero,A,alpha,A,A) 00960 00961 The result is that the operation has to be "atomic". That is, 00962 B and C are no longer reliable after D is modified, so that 00963 the assignment to D must be the last thing to occur. 00964 00965 D = alpha*B + beta*C 00966 00967 1) Use alpha==0,beta==1 and check that D == C 00968 2) Use alpha==1,beta==0 and check that D == B 00969 *********************************************************************/ 00970 { 00971 const int p = 7; 00972 Teuchos::RCP<MV> B, D; 00973 Teuchos::RCP<const MV> C; 00974 std::vector<MagType> normsB(p), 00975 normsD(p); 00976 std::vector<int> lclindex(p); 00977 for (int i=0; i<p; i++) lclindex[i] = i; 00978 00979 B = MVT::Clone(*A,p); 00980 C = MVT::CloneView(*B,lclindex); 00981 D = MVT::CloneViewNonConst(*B,lclindex); 00982 00983 MVT::MvRandom(*B); 00984 MVT::MvNorm(*B,normsB); 00985 00986 // check that 0*B+1*C == C 00987 MVT::MvAddMv(zero,*B,one,*C,*D); 00988 MVT::MvNorm(*D,normsD); 00989 for (int i=0; i<p; i++) { 00990 if ( normsB[i] != normsD[i] ) { 00991 om->stream(Warnings) 00992 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl 00993 << "Assignment did not work." << endl; 00994 return false; 00995 } 00996 } 00997 00998 // check that 1*B+0*C == B 00999 MVT::MvAddMv(one,*B,zero,*C,*D); 01000 MVT::MvNorm(*D,normsD); 01001 for (int i=0; i<p; i++) { 01002 if ( normsB[i] != normsD[i] ) { 01003 om->stream(Warnings) 01004 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl 01005 << "Assignment did not work." << endl; 01006 return false; 01007 } 01008 } 01009 01010 } 01011 01012 01013 /*********** MvTimesMatAddMv() 7 by 5 ******************************** 01014 C = alpha*B*SDM + beta*C 01015 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged 01016 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero 01017 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B 01018 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged 01019 5) Test with non-square matrices 01020 6) Always check that input arguments are not modified 01021 *********************************************************************/ 01022 { 01023 const int p = 7, q = 5; 01024 Teuchos::RCP<MV> B, C; 01025 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 01026 std::vector<MagType> normsC1(q), normsC2(q), 01027 normsB1(p), normsB2(p); 01028 01029 B = MVT::Clone(*A,p); 01030 C = MVT::Clone(*A,q); 01031 01032 // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged 01033 MVT::MvRandom(*B); 01034 MVT::MvRandom(*C); 01035 MVT::MvNorm(*B,normsB1); 01036 MVT::MvNorm(*C,normsC1); 01037 SDM.random(); 01038 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C); 01039 MVT::MvNorm(*B,normsB2); 01040 MVT::MvNorm(*C,normsC2); 01041 for (int i=0; i<p; i++) { 01042 if ( normsB1[i] != normsB2[i] ) { 01043 om->stream(Warnings) 01044 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01045 << "Input vectors were modified." << endl; 01046 return false; 01047 } 01048 } 01049 for (int i=0; i<q; i++) { 01050 if ( normsC1[i] != normsC2[i] ) { 01051 om->stream(Warnings) 01052 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01053 << "Arithmetic test 1 failed." << endl; 01054 return false; 01055 } 01056 } 01057 01058 // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero 01059 MVT::MvRandom(*B); 01060 MVT::MvRandom(*C); 01061 MVT::MvNorm(*B,normsB1); 01062 MVT::MvNorm(*C,normsC1); 01063 SDM.random(); 01064 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C); 01065 MVT::MvNorm(*B,normsB2); 01066 MVT::MvNorm(*C,normsC2); 01067 for (int i=0; i<p; i++) { 01068 if ( normsB1[i] != normsB2[i] ) { 01069 om->stream(Warnings) 01070 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01071 << "Input vectors were modified." << endl; 01072 return false; 01073 } 01074 } 01075 for (int i=0; i<q; i++) { 01076 if ( normsC2[i] != zero ) { 01077 om->stream(Warnings) 01078 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01079 << "Arithmetic test 2 failed: " 01080 << normsC2[i] 01081 << " != " 01082 << zero 01083 << endl; 01084 return false; 01085 } 01086 } 01087 01088 // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B 01089 // |0| 01090 MVT::MvRandom(*B); 01091 MVT::MvRandom(*C); 01092 MVT::MvNorm(*B,normsB1); 01093 MVT::MvNorm(*C,normsC1); 01094 SDM.scale(zero); 01095 for (int i=0; i<q; i++) { 01096 SDM(i,i) = one; 01097 } 01098 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C); 01099 MVT::MvNorm(*B,normsB2); 01100 MVT::MvNorm(*C,normsC2); 01101 for (int i=0; i<p; i++) { 01102 if ( normsB1[i] != normsB2[i] ) { 01103 om->stream(Warnings) 01104 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01105 << "Input vectors were modified." << endl; 01106 return false; 01107 } 01108 } 01109 for (int i=0; i<q; i++) { 01110 if ( normsB1[i] != normsC2[i] ) { 01111 om->stream(Warnings) 01112 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01113 << "Arithmetic test 3 failed: " 01114 << normsB1[i] 01115 << " != " 01116 << normsC2[i] 01117 << endl; 01118 return false; 01119 } 01120 } 01121 01122 // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged 01123 MVT::MvRandom(*B); 01124 MVT::MvRandom(*C); 01125 MVT::MvNorm(*B,normsB1); 01126 MVT::MvNorm(*C,normsC1); 01127 SDM.scale(zero); 01128 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C); 01129 MVT::MvNorm(*B,normsB2); 01130 MVT::MvNorm(*C,normsC2); 01131 for (int i=0; i<p; i++) { 01132 if ( normsB1[i] != normsB2[i] ) { 01133 om->stream(Warnings) 01134 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01135 << "Input vectors were modified." << endl; 01136 return false; 01137 } 01138 } 01139 for (int i=0; i<q; i++) { 01140 if ( normsC1[i] != normsC2[i] ) { 01141 om->stream(Warnings) 01142 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01143 << "Arithmetic test 4 failed." << endl; 01144 return false; 01145 } 01146 } 01147 } 01148 01149 /*********** MvTimesMatAddMv() 5 by 7 ******************************** 01150 C = alpha*B*SDM + beta*C 01151 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged 01152 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero 01153 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B 01154 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged 01155 5) Test with non-square matrices 01156 6) Always check that input arguments are not modified 01157 *********************************************************************/ 01158 { 01159 const int p = 5, q = 7; 01160 Teuchos::RCP<MV> B, C; 01161 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q); 01162 std::vector<MagType> normsC1(q), normsC2(q), 01163 normsB1(p), normsB2(p); 01164 01165 B = MVT::Clone(*A,p); 01166 C = MVT::Clone(*A,q); 01167 01168 // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged 01169 MVT::MvRandom(*B); 01170 MVT::MvRandom(*C); 01171 MVT::MvNorm(*B,normsB1); 01172 MVT::MvNorm(*C,normsC1); 01173 SDM.random(); 01174 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C); 01175 MVT::MvNorm(*B,normsB2); 01176 MVT::MvNorm(*C,normsC2); 01177 for (int i=0; i<p; i++) { 01178 if ( normsB1[i] != normsB2[i] ) { 01179 om->stream(Warnings) 01180 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01181 << "Input vectors were modified." << endl; 01182 return false; 01183 } 01184 } 01185 for (int i=0; i<q; i++) { 01186 if ( normsC1[i] != normsC2[i] ) { 01187 om->stream(Warnings) 01188 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01189 << "Arithmetic test 5 failed." << endl; 01190 return false; 01191 } 01192 } 01193 01194 // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero 01195 MVT::MvRandom(*B); 01196 MVT::MvRandom(*C); 01197 MVT::MvNorm(*B,normsB1); 01198 MVT::MvNorm(*C,normsC1); 01199 SDM.random(); 01200 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C); 01201 MVT::MvNorm(*B,normsB2); 01202 MVT::MvNorm(*C,normsC2); 01203 for (int i=0; i<p; i++) { 01204 if ( normsB1[i] != normsB2[i] ) { 01205 om->stream(Warnings) 01206 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01207 << "Input vectors were modified." << endl; 01208 return false; 01209 } 01210 } 01211 for (int i=0; i<q; i++) { 01212 if ( normsC2[i] != zero ) { 01213 om->stream(Warnings) 01214 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01215 << "Arithmetic test 6 failed: " 01216 << normsC2[i] 01217 << " != " 01218 << zero 01219 << endl; 01220 return false; 01221 } 01222 } 01223 01224 // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B 01225 MVT::MvRandom(*B); 01226 MVT::MvRandom(*C); 01227 MVT::MvNorm(*B,normsB1); 01228 MVT::MvNorm(*C,normsC1); 01229 SDM.scale(zero); 01230 for (int i=0; i<p; i++) { 01231 SDM(i,i) = one; 01232 } 01233 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C); 01234 MVT::MvNorm(*B,normsB2); 01235 MVT::MvNorm(*C,normsC2); 01236 for (int i=0; i<p; i++) { 01237 if ( normsB1[i] != normsB2[i] ) { 01238 om->stream(Warnings) 01239 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01240 << "Input vectors were modified." << endl; 01241 return false; 01242 } 01243 } 01244 for (int i=0; i<p; i++) { 01245 if ( normsB1[i] != normsC2[i] ) { 01246 om->stream(Warnings) 01247 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01248 << "Arithmetic test 7 failed." << endl; 01249 return false; 01250 } 01251 } 01252 for (int i=p; i<q; i++) { 01253 if ( normsC2[i] != zero ) { 01254 om->stream(Warnings) 01255 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01256 << "Arithmetic test 7 failed." << endl; 01257 return false; 01258 } 01259 } 01260 01261 // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged 01262 MVT::MvRandom(*B); 01263 MVT::MvRandom(*C); 01264 MVT::MvNorm(*B,normsB1); 01265 MVT::MvNorm(*C,normsC1); 01266 SDM.scale(zero); 01267 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C); 01268 MVT::MvNorm(*B,normsB2); 01269 MVT::MvNorm(*C,normsC2); 01270 for (int i=0; i<p; i++) { 01271 if ( normsB1[i] != normsB2[i] ) { 01272 om->stream(Warnings) 01273 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01274 << "Input vectors were modified." << endl; 01275 return false; 01276 } 01277 } 01278 for (int i=0; i<q; i++) { 01279 if ( normsC1[i] != normsC2[i] ) { 01280 om->stream(Warnings) 01281 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl 01282 << "Arithmetic test 8 failed." << endl; 01283 return false; 01284 } 01285 } 01286 } 01287 01288 return true; 01289 01290 } 01291 01292 01293 01299 template< class ScalarType, class MV, class OP> 01300 bool TestOperatorTraits( 01301 const Teuchos::RCP<OutputManager<ScalarType> > &om, 01302 const Teuchos::RCP<const MV> &A, 01303 const Teuchos::RCP<const OP> &M) { 01304 01305 using std::endl; 01306 01307 /* OPT Contract: 01308 Apply() 01309 MV: OP*zero == zero 01310 Warn if OP is not deterministic (OP*A != OP*A) 01311 Does not modify input arguments 01312 *********************************************************************/ 01313 01314 typedef MultiVecTraits<ScalarType, MV> MVT; 01315 typedef Teuchos::ScalarTraits<ScalarType> SCT; 01316 typedef OperatorTraits<ScalarType, MV, OP> OPT; 01317 typedef typename SCT::magnitudeType MagType; 01318 01319 const int numvecs = 10; 01320 01321 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs), 01322 C = MVT::Clone(*A,numvecs); 01323 01324 std::vector<MagType> normsB1(numvecs), normsB2(numvecs), 01325 normsC1(numvecs), normsC2(numvecs); 01326 bool NonDeterministicWarning; 01327 01328 /*********** Apply() ************************************************* 01329 Verify: 01330 1) OP*B == OP*B; OP is deterministic (just warn on this) 01331 2) OP*zero == 0 01332 3) OP*B doesn't modify B 01333 4) OP*B is invariant under initial state of destination vectors 01334 *********************************************************************/ 01335 MVT::MvInit(*B); 01336 MVT::MvRandom(*C); 01337 MVT::MvNorm(*B,normsB1); 01338 OPT::Apply(*M,*B,*C); 01339 MVT::MvNorm(*B,normsB2); 01340 MVT::MvNorm(*C,normsC2); 01341 for (int i=0; i<numvecs; i++) { 01342 if (normsB2[i] != normsB1[i]) { 01343 om->stream(Warnings) 01344 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl 01345 << "Apply() modified the input vectors." << endl; 01346 return false; 01347 } 01348 if (normsC2[i] != SCT::zero()) { 01349 om->stream(Warnings) 01350 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl 01351 << "Operator applied to zero did not return zero." << endl; 01352 return false; 01353 } 01354 } 01355 01356 // If we send in a random matrix, we should not get a zero return 01357 MVT::MvRandom(*B); 01358 MVT::MvNorm(*B,normsB1); 01359 OPT::Apply(*M,*B,*C); 01360 MVT::MvNorm(*B,normsB2); 01361 MVT::MvNorm(*C,normsC2); 01362 bool ZeroWarning = false; 01363 for (int i=0; i<numvecs; i++) { 01364 if (normsB2[i] != normsB1[i]) { 01365 om->stream(Warnings) 01366 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl 01367 << "Apply() modified the input vectors." << endl; 01368 return false; 01369 } 01370 if (normsC2[i] == SCT::zero() && ZeroWarning==false ) { 01371 om->stream(Warnings) 01372 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl 01373 << "Operator applied to random vectors returned zero." << endl; 01374 ZeroWarning = true; 01375 } 01376 } 01377 01378 // Apply operator with C init'd to zero 01379 MVT::MvRandom(*B); 01380 MVT::MvNorm(*B,normsB1); 01381 MVT::MvInit(*C); 01382 OPT::Apply(*M,*B,*C); 01383 MVT::MvNorm(*B,normsB2); 01384 MVT::MvNorm(*C,normsC1); 01385 for (int i=0; i<numvecs; i++) { 01386 if (normsB2[i] != normsB1[i]) { 01387 om->stream(Warnings) 01388 << "*** ERROR *** OperatorTraits::Apply() [3]" << endl 01389 << "Apply() modified the input vectors." << endl; 01390 return false; 01391 } 01392 } 01393 01394 // Apply operator with C init'd to random 01395 // Check that result is the same as before; warn if not. 01396 // This could be a result of a bug, or a stochastic 01397 // operator. We do not want to prejudice against a 01398 // stochastic operator. 01399 MVT::MvRandom(*C); 01400 OPT::Apply(*M,*B,*C); 01401 MVT::MvNorm(*B,normsB2); 01402 MVT::MvNorm(*C,normsC2); 01403 NonDeterministicWarning = false; 01404 for (int i=0; i<numvecs; i++) { 01405 if (normsB2[i] != normsB1[i]) { 01406 om->stream(Warnings) 01407 << "*** ERROR *** OperatorTraits::Apply() [4]" << endl 01408 << "Apply() modified the input vectors." << endl; 01409 return false; 01410 } 01411 if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) { 01412 om->stream(Warnings) 01413 << endl 01414 << "*** WARNING *** OperatorTraits::Apply() [4]" << endl 01415 << "Apply() returned two different results." << endl << endl; 01416 NonDeterministicWarning = true; 01417 } 01418 } 01419 01420 return true; 01421 01422 } 01423 01424 } 01425 01426 #endif
1.7.6.1
