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     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 C is zero
00646         5) Zero A and check that C is zero
00647 
00648         Note: Should we really require that C is correctly sized already?
00649         Epetra does (and crashes if it isn't.)
00650     *********************************************************************/
00651     {
00652       const int p = 7;
00653       const int q = 9;
00654       Teuchos::RCP<MV> B, C;
00655       std::vector<MagType> normsB(p), normsC(q);
00656       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
00657 
00658       B = MVT::Clone(*A,p);
00659       C = MVT::Clone(*A,q);
00660    
00661       // randomize the multivectors
00662       MVT::MvRandom(*B);
00663       MVT::MvNorm(*B,normsB);
00664       MVT::MvRandom(*C);
00665       MVT::MvNorm(*C,normsC);
00666    
00667       // perform SDM  = zero() * B^H * C
00668       MVT::MvTransMv( zero, *B, *C, SDM );
00669    
00670       // check the sizes: not allowed to have shrunk
00671       if ( SDM.numRows() != p || SDM.numCols() != q ) {
00672         om->stream(Warnings)
00673           << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00674           << "Routine resized SerialDenseMatrix." << endl;
00675         return false;
00676       }
00677    
00678       // check that zero**A^H*B == zero
00679       if ( SDM.normOne() != zero ) {
00680         om->stream(Warnings)
00681           << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00682           << "Scalar argument processed incorrectly." << endl;
00683         return false;
00684       }
00685    
00686       // perform SDM  = one * B^H * C
00687       MVT::MvTransMv( one, *B, *C, SDM );
00688    
00689       // check the norms: a^H b = |a| |b| cos(theta) <= |a| |b|
00690       // with equality only when a and b are colinear
00691       for (int i=0; i<p; i++) {
00692         for (int j=0; j<q; j++) {
00693           if (   SCT::magnitude(SDM(i,j)) 
00694                > SCT::magnitude(normsB[i]*normsC[j]) ) {
00695             om->stream(Warnings)
00696               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00697               << "Triangle inequality did not hold: " 
00698               << SCT::magnitude(SDM(i,j)) 
00699               << " > " 
00700               << SCT::magnitude(normsB[i]*normsC[j]) 
00701               << endl;
00702             return false;
00703           }
00704         }
00705       }
00706       MVT::MvInit(*C);
00707       MVT::MvRandom(*B);
00708       MVT::MvTransMv( one, *B, *C, SDM );
00709       for (int i=0; i<p; i++) {
00710         for (int j=0; j<q; j++) {
00711           if ( SDM(i,j) != zero ) {
00712             om->stream(Warnings)
00713               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00714               << "Inner products not zero for C==0." << endl;
00715             return false;
00716           }
00717         }
00718       }
00719       MVT::MvInit(*B);
00720       MVT::MvRandom(*C);
00721       MVT::MvTransMv( one, *B, *C, SDM );
00722       for (int i=0; i<p; i++) {
00723         for (int j=0; j<q; j++) {
00724           if ( SDM(i,j) != zero ) {
00725             om->stream(Warnings)
00726               << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
00727               << "Inner products not zero for B==0." << endl;
00728             return false;
00729           }
00730         }
00731       }
00732     } 
00733 
00734 
00735     /*********** MvDot() *************************************************
00736         Verify:
00737         1) Results vector not resized
00738         2) Inner product inequalities are satisfied
00739         3) Zero vectors give zero inner products
00740     *********************************************************************/
00741     {
00742       const int p = 7;
00743       Teuchos::RCP<MV> B, C;
00744       std::vector<ScalarType> iprods(p);
00745       std::vector<MagType> normsB(p), normsC(p);
00746 
00747       B = MVT::Clone(*A,p);
00748       C = MVT::Clone(*A,p);
00749 
00750       MVT::MvRandom(*B);
00751       MVT::MvRandom(*C);
00752       MVT::MvNorm(*B,normsB);
00753       MVT::MvNorm(*C,normsC);
00754       MVT::MvDot( *B, *C, iprods );
00755       if ( (int)iprods.size() != p ) {
00756         om->stream(Warnings)
00757           << "*** ERROR *** MultiVecTraits::MvDot." << endl
00758           << "Routine resized results vector." << endl;
00759         return false;
00760       }
00761       for (int i=0; i<p; i++) {
00762         if ( SCT::magnitude(iprods[i]) 
00763              > SCT::magnitude(normsB[i]*normsC[i]) ) {
00764           om->stream(Warnings)
00765             << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00766             << "Inner products not valid." << endl;
00767           return false;
00768         }
00769       }
00770       MVT::MvInit(*B);
00771       MVT::MvRandom(*C);
00772       MVT::MvDot( *B, *C, iprods );
00773       for (int i=0; i<p; i++) {
00774         if ( iprods[i] != zero ) {
00775           om->stream(Warnings)
00776             << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00777             << "Inner products not zero for B==0." << endl;
00778           return false;
00779         }
00780       }
00781       MVT::MvInit(*C);
00782       MVT::MvRandom(*B);
00783       MVT::MvDot( *B, *C, iprods );
00784       for (int i=0; i<p; i++) {
00785         if ( iprods[i] != zero ) {
00786           om->stream(Warnings)
00787             << "*** ERROR *** MultiVecTraits::MvDot()." << endl
00788             << "Inner products not zero for C==0." << endl;
00789           return false;
00790         }
00791       }
00792     }
00793 
00794 
00795     /*********** MvAddMv() ***********************************************
00796        D = alpha*B + beta*C
00797        1) Use alpha==0,beta==1 and check that D == C
00798        2) Use alpha==1,beta==0 and check that D == B
00799        3) Use D==0 and D!=0 and check that result is the same
00800        4) Check that input arguments are not modified 
00801     *********************************************************************/
00802     {
00803       const int p = 7;
00804       Teuchos::RCP<MV> B, C, D;
00805       std::vector<MagType> normsB1(p), normsB2(p),
00806                            normsC1(p), normsC2(p),
00807                            normsD1(p), normsD2(p);
00808       ScalarType alpha = SCT::random(),
00809                   beta = SCT::random();
00810 
00811       B = MVT::Clone(*A,p);
00812       C = MVT::Clone(*A,p);
00813       D = MVT::Clone(*A,p);
00814 
00815       MVT::MvRandom(*B);
00816       MVT::MvRandom(*C);
00817       MVT::MvNorm(*B,normsB1);
00818       MVT::MvNorm(*C,normsC1);
00819    
00820       // check that 0*B+1*C == C
00821       MVT::MvAddMv(zero,*B,one,*C,*D);
00822       MVT::MvNorm(*B,normsB2);
00823       MVT::MvNorm(*C,normsC2);
00824       MVT::MvNorm(*D,normsD1);
00825       for (int i=0; i<p; i++) {
00826         if ( normsB1[i] != normsB2[i] ) {
00827           om->stream(Warnings)
00828             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00829             << "Input arguments were modified." << endl;
00830           return false;
00831         }
00832         else if ( normsC1[i] != normsC2[i] ) {
00833           om->stream(Warnings)
00834             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00835             << "Input arguments were modified." << endl;
00836           return false;
00837         }
00838         else if ( normsC1[i] != normsD1[i] ) {
00839           om->stream(Warnings)
00840             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00841             << "Assignment did not work." << endl;
00842           return false;
00843         }
00844       }
00845 
00846       // check that 1*B+0*C == B
00847       MVT::MvAddMv(one,*B,zero,*C,*D);
00848       MVT::MvNorm(*B,normsB2);
00849       MVT::MvNorm(*C,normsC2);
00850       MVT::MvNorm(*D,normsD1);
00851       for (int i=0; i<p; i++) {
00852         if ( normsB1[i] != normsB2[i] ) {
00853           om->stream(Warnings)
00854             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00855             << "Input arguments were modified." << endl;
00856           return false;
00857         }
00858         else if ( normsC1[i] != normsC2[i] ) {
00859           om->stream(Warnings)
00860             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00861             << "Input arguments were modified." << endl;
00862           return false;
00863         }
00864         else if ( normsB1[i] != normsD1[i] ) {
00865           om->stream(Warnings)
00866             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00867             << "Assignment did not work." << endl;
00868           return false;
00869         }
00870       }
00871 
00872       // check that alpha*B+beta*C -> D is invariant under initial D
00873       // first, try random D
00874       MVT::MvRandom(*D);
00875       MVT::MvAddMv(alpha,*B,beta,*C,*D);
00876       MVT::MvNorm(*B,normsB2);
00877       MVT::MvNorm(*C,normsC2);
00878       MVT::MvNorm(*D,normsD1);
00879       // check that input args are not modified
00880       for (int i=0; i<p; i++) {
00881         if ( normsB1[i] != normsB2[i] ) {
00882           om->stream(Warnings)
00883             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00884             << "Input arguments were modified." << endl;
00885           return false;
00886         }
00887         else if ( normsC1[i] != normsC2[i] ) {
00888           om->stream(Warnings)
00889             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00890             << "Input arguments were modified." << endl;
00891           return false;
00892         }
00893       }
00894       // next, try zero D
00895       MVT::MvInit(*D);
00896       MVT::MvAddMv(alpha,*B,beta,*C,*D);
00897       MVT::MvNorm(*B,normsB2);
00898       MVT::MvNorm(*C,normsC2);
00899       MVT::MvNorm(*D,normsD2);
00900       // check that input args are not modified and that D is the same
00901       // as the above test
00902       for (int i=0; i<p; i++) {
00903         if ( normsB1[i] != normsB2[i] ) {
00904           om->stream(Warnings)
00905             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00906             << "Input arguments were modified." << endl;
00907           return false;
00908         }
00909         else if ( normsC1[i] != normsC2[i] ) {
00910           om->stream(Warnings)
00911             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00912             << "Input arguments were modified." << endl;
00913           return false;
00914         }
00915         else if ( normsD1[i] != normsD2[i] ) {
00916           om->stream(Warnings)
00917             << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
00918             << "Results varies depending on initial state of dest vectors." << endl;
00919           return false;
00920         }
00921       }
00922     }
00923 
00924     /*********** MvAddMv() ***********************************************
00925        Similar to above, but where B or C are potentially the same 
00926        object as D. This case is commonly used, for example, to affect
00927        A <- alpha*A
00928        via 
00929        MvAddMv(alpha,A,zero,A,A)
00930           ** OR **
00931        MvAddMv(zero,A,alpha,A,A)
00932       
00933        The result is that the operation has to be "atomic". That is, 
00934        B and C are no longer reliable after D is modified, so that 
00935        the assignment to D must be the last thing to occur.
00936 
00937        D = alpha*B + beta*C
00938 
00939        1) Use alpha==0,beta==1 and check that D == C
00940        2) Use alpha==1,beta==0 and check that D == B
00941     *********************************************************************/
00942     {
00943       const int p = 7;
00944       Teuchos::RCP<MV> B, D;
00945       Teuchos::RCP<const MV> C;
00946       std::vector<MagType> normsB(p),
00947                            normsD(p);
00948       std::vector<int> lclindex(p);
00949       for (int i=0; i<p; i++) lclindex[i] = i;
00950 
00951       B = MVT::Clone(*A,p);
00952       C = MVT::CloneView(*B,lclindex);
00953       D = MVT::CloneViewNonConst(*B,lclindex);
00954 
00955       MVT::MvRandom(*B);
00956       MVT::MvNorm(*B,normsB);
00957    
00958       // check that 0*B+1*C == C
00959       MVT::MvAddMv(zero,*B,one,*C,*D);
00960       MVT::MvNorm(*D,normsD);
00961       for (int i=0; i<p; i++) {
00962         if ( normsB[i] != normsD[i] ) {
00963           om->stream(Warnings)
00964             << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
00965             << "Assignment did not work." << endl;
00966           return false;
00967         }
00968       }
00969 
00970       // check that 1*B+0*C == B
00971       MVT::MvAddMv(one,*B,zero,*C,*D);
00972       MVT::MvNorm(*D,normsD);
00973       for (int i=0; i<p; i++) {
00974         if ( normsB[i] != normsD[i] ) {
00975           om->stream(Warnings)
00976             << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
00977             << "Assignment did not work." << endl;
00978           return false;
00979         }
00980       }
00981 
00982     }
00983 
00984 
00985     /*********** MvTimesMatAddMv() 7 by 5 ********************************
00986        C = alpha*B*SDM + beta*C
00987        1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
00988        2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
00989        3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
00990        4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
00991        5) Test with non-square matrices
00992        6) Always check that input arguments are not modified 
00993     *********************************************************************/
00994     {
00995       const int p = 7, q = 5;
00996       Teuchos::RCP<MV> B, C;
00997       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
00998       std::vector<MagType> normsC1(q), normsC2(q),
00999                            normsB1(p), normsB2(p);
01000       
01001       B = MVT::Clone(*A,p);
01002       C = MVT::Clone(*A,q);
01003 
01004       // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged
01005       MVT::MvRandom(*B);
01006       MVT::MvRandom(*C);
01007       MVT::MvNorm(*B,normsB1);
01008       MVT::MvNorm(*C,normsC1);
01009       SDM.random();
01010       MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
01011       MVT::MvNorm(*B,normsB2);
01012       MVT::MvNorm(*C,normsC2);
01013       for (int i=0; i<p; i++) {
01014         if ( normsB1[i] != normsB2[i] ) {
01015           om->stream(Warnings)
01016             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01017             << "Input vectors were modified." << endl;
01018           return false;
01019         }
01020       }
01021       for (int i=0; i<q; i++) {
01022         if ( normsC1[i] != normsC2[i] ) {
01023           om->stream(Warnings)
01024             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01025             << "Arithmetic test 1 failed." << endl;
01026           return false;
01027         }
01028       }
01029 
01030       // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero
01031       MVT::MvRandom(*B);
01032       MVT::MvRandom(*C);
01033       MVT::MvNorm(*B,normsB1);
01034       MVT::MvNorm(*C,normsC1);
01035       SDM.random();
01036       MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01037       MVT::MvNorm(*B,normsB2);
01038       MVT::MvNorm(*C,normsC2);
01039       for (int i=0; i<p; i++) {
01040         if ( normsB1[i] != normsB2[i] ) {
01041           om->stream(Warnings)
01042             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01043             << "Input vectors were modified." << endl;
01044           return false;
01045         }
01046       }
01047       for (int i=0; i<q; i++) {
01048         if ( normsC2[i] != zero ) {
01049           om->stream(Warnings)
01050             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01051             << "Arithmetic test 2 failed: " 
01052             << normsC2[i] 
01053             << " != " 
01054             << zero 
01055             << endl;
01056           return false;
01057         }
01058       }
01059 
01060       // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B
01061       //                        |0|
01062       MVT::MvRandom(*B);
01063       MVT::MvRandom(*C);
01064       MVT::MvNorm(*B,normsB1);
01065       MVT::MvNorm(*C,normsC1);
01066       SDM.scale(zero);
01067       for (int i=0; i<q; i++) {
01068         SDM(i,i) = one;
01069       }
01070       MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01071       MVT::MvNorm(*B,normsB2);
01072       MVT::MvNorm(*C,normsC2);
01073       for (int i=0; i<p; i++) {
01074         if ( normsB1[i] != normsB2[i] ) {
01075           om->stream(Warnings)
01076             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01077             << "Input vectors were modified." << endl;
01078           return false;
01079         }
01080       }
01081       for (int i=0; i<q; i++) {
01082         if ( normsB1[i] != normsC2[i] ) {
01083           om->stream(Warnings)
01084             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01085             << "Arithmetic test 3 failed: "
01086             << normsB1[i] 
01087             << " != "
01088             << normsC2[i]
01089             << endl;
01090           return false;
01091         }
01092       }
01093 
01094       // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged
01095       MVT::MvRandom(*B);
01096       MVT::MvRandom(*C);
01097       MVT::MvNorm(*B,normsB1);
01098       MVT::MvNorm(*C,normsC1);
01099       SDM.scale(zero);
01100       MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01101       MVT::MvNorm(*B,normsB2);
01102       MVT::MvNorm(*C,normsC2);
01103       for (int i=0; i<p; i++) {
01104         if ( normsB1[i] != normsB2[i] ) {
01105           om->stream(Warnings)
01106             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01107             << "Input vectors were modified." << endl;
01108           return false;
01109         }
01110       }
01111       for (int i=0; i<q; i++) {
01112         if ( normsC1[i] != normsC2[i] ) {
01113           om->stream(Warnings)
01114             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01115             << "Arithmetic test 4 failed." << endl;
01116           return false;
01117         }
01118       }
01119     }
01120 
01121     /*********** MvTimesMatAddMv() 5 by 7 ********************************
01122        C = alpha*B*SDM + beta*C
01123        1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
01124        2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
01125        3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
01126        4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
01127        5) Test with non-square matrices
01128        6) Always check that input arguments are not modified 
01129     *********************************************************************/
01130     {
01131       const int p = 5, q = 7;
01132       Teuchos::RCP<MV> B, C;
01133       Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
01134       std::vector<MagType> normsC1(q), normsC2(q),
01135                            normsB1(p), normsB2(p);
01136       
01137       B = MVT::Clone(*A,p);
01138       C = MVT::Clone(*A,q);
01139 
01140       // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged
01141       MVT::MvRandom(*B);
01142       MVT::MvRandom(*C);
01143       MVT::MvNorm(*B,normsB1);
01144       MVT::MvNorm(*C,normsC1);
01145       SDM.random();
01146       MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
01147       MVT::MvNorm(*B,normsB2);
01148       MVT::MvNorm(*C,normsC2);
01149       for (int i=0; i<p; i++) {
01150         if ( normsB1[i] != normsB2[i] ) {
01151           om->stream(Warnings)
01152             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01153             << "Input vectors were modified." << endl;
01154           return false;
01155         }
01156       }
01157       for (int i=0; i<q; i++) {
01158         if ( normsC1[i] != normsC2[i] ) {
01159           om->stream(Warnings)
01160             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01161             << "Arithmetic test 5 failed." << endl;
01162           return false;
01163         }
01164       }
01165 
01166       // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero
01167       MVT::MvRandom(*B);
01168       MVT::MvRandom(*C);
01169       MVT::MvNorm(*B,normsB1);
01170       MVT::MvNorm(*C,normsC1);
01171       SDM.random();
01172       MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
01173       MVT::MvNorm(*B,normsB2);
01174       MVT::MvNorm(*C,normsC2);
01175       for (int i=0; i<p; i++) {
01176         if ( normsB1[i] != normsB2[i] ) {
01177           om->stream(Warnings)
01178             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01179             << "Input vectors were modified." << endl;
01180           return false;
01181         }
01182       }
01183       for (int i=0; i<q; i++) {
01184         if ( normsC2[i] != zero ) {
01185           om->stream(Warnings)
01186             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01187             << "Arithmetic test 6 failed: " 
01188             << normsC2[i] 
01189             << " != " 
01190             << zero 
01191             << endl;
01192           return false;
01193         }
01194       }
01195 
01196       // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B
01197       MVT::MvRandom(*B);
01198       MVT::MvRandom(*C);
01199       MVT::MvNorm(*B,normsB1);
01200       MVT::MvNorm(*C,normsC1);
01201       SDM.scale(zero);
01202       for (int i=0; i<p; i++) {
01203         SDM(i,i) = one;
01204       }
01205       MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
01206       MVT::MvNorm(*B,normsB2);
01207       MVT::MvNorm(*C,normsC2);
01208       for (int i=0; i<p; i++) {
01209         if ( normsB1[i] != normsB2[i] ) {
01210           om->stream(Warnings)
01211             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01212             << "Input vectors were modified." << endl;
01213           return false;
01214         }
01215       }
01216       for (int i=0; i<p; i++) {
01217         if ( normsB1[i] != normsC2[i] ) {
01218           om->stream(Warnings)
01219             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01220             << "Arithmetic test 7 failed." << endl;
01221           return false;
01222         }
01223       }
01224       for (int i=p; i<q; i++) {
01225         if ( normsC2[i] != zero ) {
01226           om->stream(Warnings)
01227             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01228             << "Arithmetic test 7 failed." << endl;
01229           return false;
01230         }
01231       }
01232 
01233       // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged
01234       MVT::MvRandom(*B);
01235       MVT::MvRandom(*C);
01236       MVT::MvNorm(*B,normsB1);
01237       MVT::MvNorm(*C,normsC1);
01238       SDM.scale(zero);
01239       MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
01240       MVT::MvNorm(*B,normsB2);
01241       MVT::MvNorm(*C,normsC2);
01242       for (int i=0; i<p; i++) {
01243         if ( normsB1[i] != normsB2[i] ) {
01244           om->stream(Warnings)
01245             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01246             << "Input vectors were modified." << endl;
01247           return false;
01248         }
01249       }
01250       for (int i=0; i<q; i++) {
01251         if ( normsC1[i] != normsC2[i] ) {
01252           om->stream(Warnings)
01253             << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
01254             << "Arithmetic test 8 failed." << endl;
01255           return false;
01256         }
01257       }
01258     }
01259 
01260     return true;
01261 
01262   }
01263 
01264 
01265 
01271   template< class ScalarType, class MV, class OP>
01272   bool TestOperatorTraits( 
01273                 const Teuchos::RCP<OutputManager<ScalarType> > &om,
01274                 const Teuchos::RCP<const MV> &A,
01275                 const Teuchos::RCP<const OP> &M) {
01276 
01277     using std::endl;
01278 
01279     /* OPT Contract:
01280        Apply()
01281          MV: OP*zero == zero
01282              Warn if OP is not deterministic (OP*A != OP*A)
01283              Does not modify input arguments
01284     *********************************************************************/
01285 
01286     typedef MultiVecTraits<ScalarType, MV>     MVT;
01287     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
01288     typedef OperatorTraits<ScalarType, MV, OP> OPT;
01289     typedef typename SCT::magnitudeType        MagType;
01290 
01291     const int numvecs = 10;
01292 
01293     Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs), 
01294                              C = MVT::Clone(*A,numvecs);
01295 
01296     std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
01297                          normsC1(numvecs), normsC2(numvecs);
01298     bool NonDeterministicWarning;
01299 
01300     /*********** Apply() *************************************************
01301         Verify:
01302         1) OP*B == OP*B; OP is deterministic (just warn on this)
01303         2) OP*zero == 0
01304         3) OP*B doesn't modify B
01305         4) OP*B is invariant under initial state of destination vectors
01306     *********************************************************************/
01307     MVT::MvInit(*B);
01308     MVT::MvRandom(*C);
01309     MVT::MvNorm(*B,normsB1);
01310     OPT::Apply(*M,*B,*C);
01311     MVT::MvNorm(*B,normsB2);
01312     MVT::MvNorm(*C,normsC2);
01313     for (int i=0; i<numvecs; i++) {
01314       if (normsB2[i] != normsB1[i]) {
01315         om->stream(Warnings)
01316           << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01317           << "Apply() modified the input vectors." << endl;
01318         return false;
01319       }
01320       if (normsC2[i] != SCT::zero()) {
01321         om->stream(Warnings)
01322           << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
01323           << "Operator applied to zero did not return zero." << endl;
01324         return false;
01325       }
01326     }
01327 
01328     // If we send in a random matrix, we should not get a zero return
01329     MVT::MvRandom(*B);
01330     MVT::MvNorm(*B,normsB1);
01331     OPT::Apply(*M,*B,*C);
01332     MVT::MvNorm(*B,normsB2);
01333     MVT::MvNorm(*C,normsC2);
01334     bool ZeroWarning = false;
01335     for (int i=0; i<numvecs; i++) {
01336       if (normsB2[i] != normsB1[i]) {
01337         om->stream(Warnings)
01338           << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01339           << "Apply() modified the input vectors." << endl;
01340         return false;
01341       }
01342       if (normsC2[i] == SCT::zero() && ZeroWarning==false ) {
01343         om->stream(Warnings)
01344           << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
01345           << "Operator applied to random vectors returned zero." << endl;
01346         ZeroWarning = true;
01347       }
01348     }
01349 
01350     // Apply operator with C init'd to zero
01351     MVT::MvRandom(*B);
01352     MVT::MvNorm(*B,normsB1);
01353     MVT::MvInit(*C);
01354     OPT::Apply(*M,*B,*C);
01355     MVT::MvNorm(*B,normsB2);
01356     MVT::MvNorm(*C,normsC1);
01357     for (int i=0; i<numvecs; i++) {
01358       if (normsB2[i] != normsB1[i]) {
01359         om->stream(Warnings)
01360           << "*** ERROR *** OperatorTraits::Apply() [3]" << endl
01361           << "Apply() modified the input vectors." << endl;
01362         return false;
01363       }
01364     }
01365 
01366     // Apply operator with C init'd to random
01367     // Check that result is the same as before; warn if not.
01368     // This could be a result of a bug, or a stochastic
01369     //   operator. We do not want to prejudice against a 
01370     //   stochastic operator.
01371     MVT::MvRandom(*C);
01372     OPT::Apply(*M,*B,*C);
01373     MVT::MvNorm(*B,normsB2);
01374     MVT::MvNorm(*C,normsC2);
01375     NonDeterministicWarning = false;
01376     for (int i=0; i<numvecs; i++) {
01377       if (normsB2[i] != normsB1[i]) {
01378         om->stream(Warnings)
01379           << "*** ERROR *** OperatorTraits::Apply() [4]" << endl
01380           << "Apply() modified the input vectors." << endl;
01381         return false;
01382       }
01383       if (normsC1[i] != normsC2[i] && !NonDeterministicWarning) {
01384         om->stream(Warnings)
01385           << endl
01386           << "*** WARNING *** OperatorTraits::Apply() [4]" << endl
01387           << "Apply() returned two different results." << endl << endl;
01388         NonDeterministicWarning = true;
01389       }
01390     }
01391 
01392     return true;
01393 
01394   }
01395 
01396 }
01397 
01398 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 09:56:58 2011 for Anasazi by  doxygen 1.6.3