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