AnasaziMVOPTester.hpp

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

Generated on Wed May 12 21:24:34 2010 for Anasazi by  doxygen 1.4.7