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