Anasazi Version of the Day
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 int numvecs   = 10;
00153     const 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     TEUCHOS_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 ( 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 ( (int)norms.size() != numvecs ) {
00219         om->stream(Warnings)
00220           << "*** ERROR *** MultiVecTraits::MvNorm()." << endl
00221           << "Method resized the output vector." << endl;
00222       }
00223       for (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 (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 (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 (int i=0; i<numvecs; i++) {
00306         if ( norms[i] != zero_mag ) {
00307           om->stream(Warnings)
00308             << "*** ERROR *** MultiVecTraits::MvScale(alpha) "
00309             << "Supposedly zero vector has non-zero norm." << endl;
00310           return false;
00311         }
00312       }
00313 
00314       MVT::MvRandom(*B);
00315       std::vector<ScalarType> zeros(numvecs,SCT::zero());
00316       MVT::MvScale(*B,zeros);
00317       MVT::MvNorm(*B, norms);
00318       for (int i=0; i<numvecs; i++) {
00319         if ( norms[i] != zero_mag ) {
00320           om->stream(Warnings)
00321             << "*** ERROR *** MultiVecTraits::MvScale(alphas) "
00322             << "Supposedly zero vector has non-zero norm." << endl;
00323           return false;
00324         }
00325       }
00326     }
00327 
00328 
00329     /*********** MvInit() and MvNorm() ***********************************
00330        A 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 (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 (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 ( 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 (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 (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 ( MVT::GetNumberVecs(*C) != numvecs ) {
00460         om->stream(Warnings)
00461           << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00462           << "Wrong number of vectors." << endl;
00463         return false;
00464       }
00465       for (int i=0; i<numvecs; i++) {
00466         if ( norms2[i] != norms[i] ) {
00467           om->stream(Warnings)
00468             << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00469             << "Copied vectors do not agree." << endl;
00470           return false;
00471         }
00472       }
00473       MVT::MvInit(*B,zero);
00474       MVT::MvNorm(*C, norms); 
00475       for (int i=0; i<numvecs; i++) {
00476         if ( norms2[i] != norms[i] ) {
00477           om->stream(Warnings)
00478             << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
00479             << "Copied vectors were not independent." << endl;
00480           return false;
00481         }
00482       }
00483     }
00484 
00485 
00486     /*********** CloneView(MV,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::CloneViewNonConst(*B,ind);
00500       MVT::MvNorm(*C, norms2);
00501       if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
00502         om->stream(Warnings)
00503           << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
00504           << "Wrong number of vectors." << endl;
00505         return false;
00506       }
00507       for (int i=0; i<numvecs_2; i++) {
00508         if ( norms2[i] != norms[ind[i]] ) {
00509           om->stream(Warnings)
00510             << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
00511             << "Viewed vectors do not agree." << endl;
00512           return false;
00513         }
00514       }
00515     }
00516 
00517 
00518     /*********** const CloneView(MV,vector<int>) and MvNorm() ************
00519        Check that we have a view of the selected vectors.
00520        1) Check quantity
00521        2) Check value of norms for agreement
00522        3) Zero out B and make sure that C is zerod as well
00523     *********************************************************************/
00524     {
00525       Teuchos::RCP<MV> B;
00526       Teuchos::RCP<const MV> constB, C;
00527       std::vector<MagType> normsB(numvecs), normsC(ind.size());
00528       std::vector<int> allind(numvecs);
00529       for (int i=0; i<numvecs; i++) {
00530         allind[i] = i;
00531       }
00532 
00533       B = MVT::Clone(*A,numvecs);
00534       MVT::MvRandom( *B );
00535       MVT::MvNorm(*B, normsB);
00536       C = MVT::CloneView(*B,ind);
00537       MVT::MvNorm(*C, normsC);
00538       if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
00539         om->stream(Warnings)
00540           << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
00541           << "Wrong number of vectors." << endl;
00542         return false;
00543       }
00544       for (int i=0; i<numvecs_2; i++) {
00545         if ( normsC[i] != normsB[ind[i]] ) {
00546           om->stream(Warnings)
00547             << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
00548             << "Viewed vectors do not agree." << endl;
00549           return false;
00550         }
00551       }
00552     }
00553 
00554 
00555     /*********** SetBlock() and MvNorm() *********************************
00556        SetBlock() will copy the vectors from C into B 
00557        1) Verify that the specified vectors were copied
00558        2) Verify that the other vectors were not modified
00559        3) Verify that C was not modified
00560        4) Change C and then check B to make sure it was not modified
00561       
00562        Use a different index set than has been used so far (distinct entries).
00563        This is because duplicate entries will cause the vector to be
00564        overwritten, making it more difficult to test.
00565     *********************************************************************/
00566     {
00567       Teuchos::RCP<MV> B, C;
00568       std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
00569                            normsC1(numvecs_2), normsC2(numvecs_2);
00570 
00571       B = MVT::Clone(*A,numvecs);
00572       C = MVT::Clone(*A,numvecs_2);
00573       // Just do every other one, interleaving the vectors of C into B
00574       ind.resize(numvecs_2);
00575       for (int i=0; i<numvecs_2; i++) {
00576         ind[i] = 2*i;
00577       }
00578       MVT::MvRandom(*B);
00579       MVT::MvRandom(*C);
00580 
00581       MVT::MvNorm(*B,normsB1);
00582       MVT::MvNorm(*C,normsC1);
00583       MVT::SetBlock(*C,ind,*B);
00584       MVT::MvNorm(*B,normsB2);
00585       MVT::MvNorm(*C,normsC2);
00586 
00587       // check that C was not changed by SetBlock
00588       for (int i=0; i<numvecs_2; i++) {
00589         if ( normsC1[i] != normsC2[i] ) {
00590           om->stream(Warnings)
00591             << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00592             << "Operation modified source vectors." << endl;
00593           return false;
00594         }
00595       }
00596       // check that the correct vectors of B were modified
00597       // and the others were not
00598       for (int i=0; i<numvecs; i++) {
00599         if (i % 2 == 0) {
00600           // should be a vector from C
00601           if ( normsB2[i] != normsC1[i/2] ) {
00602             om->stream(Warnings)
00603               << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00604               << "Copied vectors do not agree." << endl;
00605             return false;
00606           }
00607         }
00608         else {
00609           // should be an original vector
00610           if ( normsB1[i] != normsB2[i] ) {
00611             om->stream(Warnings)
00612               << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00613               << "Incorrect vectors were modified." << endl;
00614             return false;
00615           }
00616         }
00617       }
00618       MVT::MvInit(*C,zero);
00619       MVT::MvNorm(*B,normsB1);
00620       // verify that we copied and didn't reference
00621       for (int i=0; i<numvecs; i++) {
00622         if ( normsB1[i] != normsB2[i] ) {
00623           om->stream(Warnings)
00624             << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
00625             << "Copied vectors were not independent." << endl;
00626           return false;
00627         }
00628       }
00629     }
00630 
00631 
00632     /*********** MvTransMv() *********************************************
00633       Performs C = alpha * A^H * B, where
00634         alpha is type ScalarType
00635         A,B are type MV with p and q vectors, respectively
00636         C is a SerialDenseMatrix<int,ScalarType> ALREADY sized to p by q
00637 
00638         Verify:
00639         1) C is not resized by the routine
00640         3) Check that zero*(A^H B) == zero
00641         3) Check inner product inequality:
00642                                         [ |a1|*|b1|    ...    |ap|*|b1| ]
00643            [a1 ... ap]^H [b1 ... bq] <= [    ...    |ai|*|bj|    ...    ]
00644                                         [ |ap|*|b1|    ...    |ap|*|bq| ]
00645         4) Zero B and check that B^H * C is zero
00646         5) Zero C and check that B^H * C is zero
00647 
00648         Note 1: Test 4 is performed with a p x q Teuchos::SDM view of 
00649                 a (p+1) x (q+1) Teuchos::SDM that is initialized to ones.
00650                 This ensures the user is correctly accessing and filling the SDM.
00651 
00652         Note 2: Should we really require that C is correctly sized already?
00653                 Epetra does (and crashes if it isn't.)
00654     *********************************************************************/
00655     {
00656       const int p = 7;
00657       const int q = 9;
00658       Teuchos::RCP<MV> B, C;
00659       std::vector<MagType> normsB(p), normsC(q);
00660       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
00661 
00662       B = MVT::Clone(*A,p);
00663       C = MVT::Clone(*A,q);
00664    
00665       // randomize the multivectors
00666       MVT::MvRandom(*B);
00667       MVT::MvNorm(*B,normsB);
00668       MVT::MvRandom(*C);
00669       MVT::MvNorm(*C,normsC);
00670    
00671       // perform SDM  = zero() * B^H * C
00672       MVT::MvTransMv( zero, *B, *C, SDM );
00673    
00674       // check the sizes: not allowed to have shrunk
00675       if ( SDM.numRows() != p || SDM.numCols() != q ) {
00676         om->stream(Warnings)
00677           << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00678           << "Routine resized SerialDenseMatrix." << endl;
00679         return false;
00680       }
00681    
00682       // check that zero**A^H*B == zero
00683       if ( SDM.normOne() != zero ) {
00684         om->stream(Warnings)
00685           << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00686           << "Scalar argument processed incorrectly." << endl;
00687         return false;
00688       }
00689    
00690       // perform SDM  = one * B^H * C
00691       MVT::MvTransMv( one, *B, *C, SDM );
00692    
00693       // check the norms: a^H b = |a| |b| cos(theta) <= |a| |b|
00694       // with equality only when a and b are colinear
00695       for (int i=0; i<p; i++) {
00696         for (int j=0; j<q; j++) {
00697           if (   SCT::magnitude(SDM(i,j)) 
00698                > SCT::magnitude(normsB[i]*normsC[j]) ) {
00699             om->stream(Warnings)
00700               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00701               << "Triangle inequality did not hold: " 
00702               << SCT::magnitude(SDM(i,j)) 
00703               << " > " 
00704               << SCT::magnitude(normsB[i]*normsC[j]) 
00705               << endl;
00706             return false;
00707           }
00708         }
00709       }
00710       MVT::MvInit(*C);
00711       MVT::MvRandom(*B);
00712       MVT::MvTransMv( one, *B, *C, SDM );
00713       for (int i=0; i<p; i++) {
00714         for (int j=0; j<q; j++) {
00715           if ( SDM(i,j) != zero ) {
00716             om->stream(Warnings)
00717               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00718               << "Inner products not zero for C==0." << endl;
00719             return false;
00720           }
00721         }
00722       }
00723       MVT::MvInit(*B);
00724       MVT::MvRandom(*C);
00725       MVT::MvTransMv( one, *B, *C, SDM );
00726       for (int i=0; i<p; i++) {
00727         for (int j=0; j<q; j++) {
00728           if ( SDM(i,j) != zero ) {
00729             om->stream(Warnings)
00730               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00731               << "Inner products not zero for B==0." << endl;
00732             return false;
00733           }
00734         }
00735       }
00736 
00737       // A larger SDM is filled with ones, initially, and a smaller
00738       // view is used for the MvTransMv method.  If the smaller SDM
00739       // is not all zeroes, then the interface is improperly writing
00740       // to the matrix object.
00741       // Note:  Since we didn't fail above, we know that the general
00742       //        inner product works, but we are checking to see if it
00743       //        works for a view too.  This is common usage in Anasazi.
00744       Teuchos::SerialDenseMatrix<int, ScalarType> largeSDM(p+1,q+1);
00745       Teuchos::SerialDenseMatrix<int, ScalarType> SDM2(Teuchos::View, largeSDM, p, q);
00746       largeSDM.putScalar( one );
00747       MVT::MvInit(*C);
00748       MVT::MvRandom(*B);
00749       MVT::MvTransMv( one, *B, *C, SDM2 );
00750       for (int i=0; i<p; i++) {
00751         for (int j=0; j<q; j++) {
00752           if ( SDM2(i,j) != zero ) {
00753             om->stream(Warnings)
00754               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00755               << "Inner products not zero for C==0 when using a view into Teuchos::SerialDenseMatrix<>." << 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 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 ( (int)iprods.size() != p ) {
00784         om->stream(Warnings)
00785           << "*** ERROR *** MultiVecTraits::MvDot." << endl
00786           << "Routine resized results vector." << endl;
00787         return false;
00788       }
00789       for (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 (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 (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 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 (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 (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 (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 (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 int p = 7;
00972       Teuchos::RCP<MV> B, D;
00973       Teuchos::RCP<const MV> C;
00974       std::vector<MagType> normsB(p),
00975                            normsD(p);
00976       std::vector<int> lclindex(p);
00977       for (int i=0; i<p; i++) lclindex[i] = i;
00978 
00979       B = MVT::Clone(*A,p);
00980       C = MVT::CloneView(*B,lclindex);
00981       D = MVT::CloneViewNonConst(*B,lclindex);
00982 
00983       MVT::MvRandom(*B);
00984       MVT::MvNorm(*B,normsB);
00985    
00986       // check that 0*B+1*C == C
00987       MVT::MvAddMv(zero,*B,one,*C,*D);
00988       MVT::MvNorm(*D,normsD);
00989       for (int i=0; i<p; i++) {
00990         if ( normsB[i] != normsD[i] ) {
00991           om->stream(Warnings)
00992             << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
00993             << "Assignment did not work." << endl;
00994           return false;
00995         }
00996       }
00997 
00998       // check that 1*B+0*C == B
00999       MVT::MvAddMv(one,*B,zero,*C,*D);
01000       MVT::MvNorm(*D,normsD);
01001       for (int i=0; i<p; i++) {
01002         if ( normsB[i] != normsD[i] ) {
01003           om->stream(Warnings)
01004             << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
01005             << "Assignment did not work." << endl;
01006           return false;
01007         }
01008       }
01009 
01010     }
01011 
01012 
01013     /*********** MvTimesMatAddMv() 7 by 5 ********************************
01014        C = alpha*B*SDM + beta*C
01015        1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
01016        2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
01017        3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
01018        4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
01019        5) Test with non-square matrices
01020        6) Always check that input arguments are not modified 
01021     *********************************************************************/
01022     {
01023       const int p = 7, q = 5;
01024       Teuchos::RCP<MV> B, C;
01025       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
01026       std::vector<MagType> normsC1(q), normsC2(q),
01027                            normsB1(p), normsB2(p);
01028       
01029       B = MVT::Clone(*A,p);
01030       C = MVT::Clone(*A,q);
01031 
01032       // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged
01033       MVT::MvRandom(*B);
01034       MVT::MvRandom(*C);
01035       MVT::MvNorm(*B,normsB1);
01036       MVT::MvNorm(*C,normsC1);
01037       SDM.random();
01038       MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
01039       MVT::MvNorm(*B,normsB2);
01040       MVT::MvNorm(*C,normsC2);
01041       for (int i=0; i<p; i++) {
01042         if ( normsB1[i] != normsB2[i] ) {
01043           om->stream(Warnings)
01044             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01045             << "Input vectors were modified." << endl;
01046           return false;
01047         }
01048       }
01049       for (int i=0; i<q; i++) {
01050         if ( normsC1[i] != normsC2[i] ) {
01051           om->stream(Warnings)
01052             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01053             << "Arithmetic test 1 failed." << endl;
01054           return false;
01055         }
01056       }
01057 
01058       // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero
01059       MVT::MvRandom(*B);
01060       MVT::MvRandom(*C);
01061       MVT::MvNorm(*B,normsB1);
01062       MVT::MvNorm(*C,normsC1);
01063       SDM.random();
01064       MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01065       MVT::MvNorm(*B,normsB2);
01066       MVT::MvNorm(*C,normsC2);
01067       for (int i=0; i<p; i++) {
01068         if ( normsB1[i] != normsB2[i] ) {
01069           om->stream(Warnings)
01070             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01071             << "Input vectors were modified." << endl;
01072           return false;
01073         }
01074       }
01075       for (int i=0; i<q; i++) {
01076         if ( normsC2[i] != zero ) {
01077           om->stream(Warnings)
01078             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01079             << "Arithmetic test 2 failed: " 
01080             << normsC2[i] 
01081             << " != " 
01082             << zero 
01083             << endl;
01084           return false;
01085         }
01086       }
01087 
01088       // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B
01089       //                        |0|
01090       MVT::MvRandom(*B);
01091       MVT::MvRandom(*C);
01092       MVT::MvNorm(*B,normsB1);
01093       MVT::MvNorm(*C,normsC1);
01094       SDM.scale(zero);
01095       for (int i=0; i<q; i++) {
01096         SDM(i,i) = one;
01097       }
01098       MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01099       MVT::MvNorm(*B,normsB2);
01100       MVT::MvNorm(*C,normsC2);
01101       for (int i=0; i<p; i++) {
01102         if ( normsB1[i] != normsB2[i] ) {
01103           om->stream(Warnings)
01104             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01105             << "Input vectors were modified." << endl;
01106           return false;
01107         }
01108       }
01109       for (int i=0; i<q; i++) {
01110         if ( normsB1[i] != normsC2[i] ) {
01111           om->stream(Warnings)
01112             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01113             << "Arithmetic test 3 failed: "
01114             << normsB1[i] 
01115             << " != "
01116             << normsC2[i]
01117             << endl;
01118           return false;
01119         }
01120       }
01121 
01122       // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged
01123       MVT::MvRandom(*B);
01124       MVT::MvRandom(*C);
01125       MVT::MvNorm(*B,normsB1);
01126       MVT::MvNorm(*C,normsC1);
01127       SDM.scale(zero);
01128       MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01129       MVT::MvNorm(*B,normsB2);
01130       MVT::MvNorm(*C,normsC2);
01131       for (int i=0; i<p; i++) {
01132         if ( normsB1[i] != normsB2[i] ) {
01133           om->stream(Warnings)
01134             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01135             << "Input vectors were modified." << endl;
01136           return false;
01137         }
01138       }
01139       for (int i=0; i<q; i++) {
01140         if ( normsC1[i] != normsC2[i] ) {
01141           om->stream(Warnings)
01142             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01143             << "Arithmetic test 4 failed." << endl;
01144           return false;
01145         }
01146       }
01147     }
01148 
01149     /*********** MvTimesMatAddMv() 5 by 7 ********************************
01150        C = alpha*B*SDM + beta*C
01151        1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
01152        2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
01153        3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
01154        4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
01155        5) Test with non-square matrices
01156        6) Always check that input arguments are not modified 
01157     *********************************************************************/
01158     {
01159       const int p = 5, q = 7;
01160       Teuchos::RCP<MV> B, C;
01161       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
01162       std::vector<MagType> normsC1(q), normsC2(q),
01163                            normsB1(p), normsB2(p);
01164       
01165       B = MVT::Clone(*A,p);
01166       C = MVT::Clone(*A,q);
01167 
01168       // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged
01169       MVT::MvRandom(*B);
01170       MVT::MvRandom(*C);
01171       MVT::MvNorm(*B,normsB1);
01172       MVT::MvNorm(*C,normsC1);
01173       SDM.random();
01174       MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
01175       MVT::MvNorm(*B,normsB2);
01176       MVT::MvNorm(*C,normsC2);
01177       for (int i=0; i<p; i++) {
01178         if ( normsB1[i] != normsB2[i] ) {
01179           om->stream(Warnings)
01180             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01181             << "Input vectors were modified." << endl;
01182           return false;
01183         }
01184       }
01185       for (int i=0; i<q; i++) {
01186         if ( normsC1[i] != normsC2[i] ) {
01187           om->stream(Warnings)
01188             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01189             << "Arithmetic test 5 failed." << endl;
01190           return false;
01191         }
01192       }
01193 
01194       // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero
01195       MVT::MvRandom(*B);
01196       MVT::MvRandom(*C);
01197       MVT::MvNorm(*B,normsB1);
01198       MVT::MvNorm(*C,normsC1);
01199       SDM.random();
01200       MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01201       MVT::MvNorm(*B,normsB2);
01202       MVT::MvNorm(*C,normsC2);
01203       for (int i=0; i<p; i++) {
01204         if ( normsB1[i] != normsB2[i] ) {
01205           om->stream(Warnings)
01206             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01207             << "Input vectors were modified." << endl;
01208           return false;
01209         }
01210       }
01211       for (int i=0; i<q; i++) {
01212         if ( normsC2[i] != zero ) {
01213           om->stream(Warnings)
01214             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01215             << "Arithmetic test 6 failed: " 
01216             << normsC2[i] 
01217             << " != " 
01218             << zero 
01219             << endl;
01220           return false;
01221         }
01222       }
01223 
01224       // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B
01225       MVT::MvRandom(*B);
01226       MVT::MvRandom(*C);
01227       MVT::MvNorm(*B,normsB1);
01228       MVT::MvNorm(*C,normsC1);
01229       SDM.scale(zero);
01230       for (int i=0; i<p; i++) {
01231         SDM(i,i) = one;
01232       }
01233       MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01234       MVT::MvNorm(*B,normsB2);
01235       MVT::MvNorm(*C,normsC2);
01236       for (int i=0; i<p; i++) {
01237         if ( normsB1[i] != normsB2[i] ) {
01238           om->stream(Warnings)
01239             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01240             << "Input vectors were modified." << endl;
01241           return false;
01242         }
01243       }
01244       for (int i=0; i<p; i++) {
01245         if ( normsB1[i] != normsC2[i] ) {
01246           om->stream(Warnings)
01247             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01248             << "Arithmetic test 7 failed." << endl;
01249           return false;
01250         }
01251       }
01252       for (int i=p; i<q; i++) {
01253         if ( normsC2[i] != zero ) {
01254           om->stream(Warnings)
01255             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01256             << "Arithmetic test 7 failed." << endl;
01257           return false;
01258         }
01259       }
01260 
01261       // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged
01262       MVT::MvRandom(*B);
01263       MVT::MvRandom(*C);
01264       MVT::MvNorm(*B,normsB1);
01265       MVT::MvNorm(*C,normsC1);
01266       SDM.scale(zero);
01267       MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01268       MVT::MvNorm(*B,normsB2);
01269       MVT::MvNorm(*C,normsC2);
01270       for (int i=0; i<p; i++) {
01271         if ( normsB1[i] != normsB2[i] ) {
01272           om->stream(Warnings)
01273             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01274             << "Input vectors were modified." << endl;
01275           return false;
01276         }
01277       }
01278       for (int i=0; i<q; i++) {
01279         if ( normsC1[i] != normsC2[i] ) {
01280           om->stream(Warnings)
01281             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01282             << "Arithmetic test 8 failed." << endl;
01283           return false;
01284         }
01285       }
01286     }
01287 
01288     return true;
01289 
01290   }
01291 
01292 
01293 
01299   template< class ScalarType, class MV, class OP>
01300   bool TestOperatorTraits( 
01301                 const Teuchos::RCP<OutputManager<ScalarType> > &om,
01302                 const Teuchos::RCP<const MV> &A,
01303                 const Teuchos::RCP<const OP> &M) {
01304 
01305     using std::endl;
01306 
01307     /* OPT Contract:
01308        Apply()
01309          MV: OP*zero == zero
01310              Warn if OP is not deterministic (OP*A != OP*A)
01311              Does not modify input arguments
01312     *********************************************************************/
01313 
01314     typedef MultiVecTraits<ScalarType, MV>     MVT;
01315     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
01316     typedef OperatorTraits<ScalarType, MV, OP> OPT;
01317     typedef typename SCT::magnitudeType        MagType;
01318 
01319     const int numvecs = 10;
01320 
01321     Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs), 
01322                              C = MVT::Clone(*A,numvecs);
01323 
01324     std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
01325                          normsC1(numvecs), normsC2(numvecs);
01326     bool NonDeterministicWarning;
01327 
01328     /*********** Apply() *************************************************
01329         Verify:
01330         1) OP*B == OP*B; OP is deterministic (just warn on this)
01331         2) OP*zero == 0
01332         3) OP*B doesn't modify B
01333         4) OP*B is invariant under initial state of destination vectors
01334     *********************************************************************/
01335     MVT::MvInit(*B);
01336     MVT::MvRandom(*C);
01337     MVT::MvNorm(*B,normsB1);
01338     OPT::Apply(*M,*B,*C);
01339     MVT::MvNorm(*B,normsB2);
01340     MVT::MvNorm(*C,normsC2);
01341     for (int i=0; i<numvecs; i++) {
01342       if (normsB2[i] != normsB1[i]) {
01343         om->stream(Warnings)
01344           << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01345           << "Apply() modified the input vectors." << endl;
01346         return false;
01347       }
01348       if (normsC2[i] != SCT::zero()) {
01349         om->stream(Warnings)
01350           << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01351           << "Operator applied to zero did not return zero." << endl;
01352         return false;
01353       }
01354     }
01355 
01356     // If we send in a random matrix, we should not get a zero return
01357     MVT::MvRandom(*B);
01358     MVT::MvNorm(*B,normsB1);
01359     OPT::Apply(*M,*B,*C);
01360     MVT::MvNorm(*B,normsB2);
01361     MVT::MvNorm(*C,normsC2);
01362     bool ZeroWarning = false;
01363     for (int i=0; i<numvecs; i++) {
01364       if (normsB2[i] != normsB1[i]) {
01365         om->stream(Warnings)
01366           << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01367           << "Apply() modified the input vectors." << endl;
01368         return false;
01369       }
01370       if (normsC2[i] == SCT::zero() && ZeroWarning==false ) {
01371         om->stream(Warnings)
01372           << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01373           << "Operator applied to random vectors returned zero." << endl;
01374         ZeroWarning = true;
01375       }
01376     }
01377 
01378     // Apply operator with C init'd to zero
01379     MVT::MvRandom(*B);
01380     MVT::MvNorm(*B,normsB1);
01381     MVT::MvInit(*C);
01382     OPT::Apply(*M,*B,*C);
01383     MVT::MvNorm(*B,normsB2);
01384     MVT::MvNorm(*C,normsC1);
01385     for (int i=0; i<numvecs; i++) {
01386       if (normsB2[i] != normsB1[i]) {
01387         om->stream(Warnings)
01388           << "*** ERROR *** OperatorTraits::Apply() [3]" << endl
01389           << "Apply() modified the input vectors." << endl;
01390         return false;
01391       }
01392     }
01393 
01394     // Apply operator with C init'd to random
01395     // Check that result is the same as before; warn if not.
01396     // This could be a result of a bug, or a stochastic
01397     //   operator. We do not want to prejudice against a 
01398     //   stochastic operator.
01399     MVT::MvRandom(*C);
01400     OPT::Apply(*M,*B,*C);
01401     MVT::MvNorm(*B,normsB2);
01402     MVT::MvNorm(*C,normsC2);
01403     NonDeterministicWarning = false;
01404     for (int i=0; i<numvecs; i++) {
01405       if (normsB2[i] != normsB1[i]) {
01406         om->stream(Warnings)
01407           << "*** ERROR *** OperatorTraits::Apply() [4]" << endl
01408           << "Apply() modified the input vectors." << endl;
01409         return false;
01410       }
01411       if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) {
01412         om->stream(Warnings)
01413           << endl
01414           << "*** WARNING *** OperatorTraits::Apply() [4]" << endl
01415           << "Apply() returned two different results." << endl << endl;
01416         NonDeterministicWarning = true;
01417       }
01418     }
01419 
01420     return true;
01421 
01422   }
01423 
01424 }
01425 
01426 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends