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

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7