BelosMVOPTester.hpp

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

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7