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