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