Anasazi Version of the Day
AnasaziMVOPTester.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends