BelosMVOPTester.hpp

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

Generated on Wed May 12 21:30:08 2010 for Belos by  doxygen 1.4.7