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