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

Generated on Tue Jul 13 09:22:47 2010 for Anasazi by  doxygen 1.4.7