Belos Version of the Day
BelosOrthoManagerTest.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 
00045 
00046 #include <BelosConfigDefs.hpp>
00047 #include <BelosMultiVecTraits.hpp>
00048 #include <BelosOutputManager.hpp>
00049 #include <BelosOrthoManagerFactory.hpp>
00050 #include <Teuchos_StandardCatchMacros.hpp>
00051 #include <Teuchos_TimeMonitor.hpp>
00052 #include <iostream>
00053 #include <stdexcept>
00054 
00055 using std::endl;
00056 
00057 namespace Belos {
00058   namespace Test {
00059 
00064     template<class Scalar, class MV>
00065     class OrthoManagerBenchmarker {
00066     private:
00067       typedef Scalar scalar_type;
00068       typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00069       typedef MultiVecTraits<Scalar, MV> MVT;
00070       typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00071 
00072     public:
00078       static void 
00079       baseline (const Teuchos::RCP<const MV>& X,
00080     const int numCols,
00081     const int numBlocks,
00082     const int numTrials)
00083       {
00084   using Teuchos::Array;
00085   using Teuchos::RCP;
00086   using Teuchos::rcp;
00087   using Teuchos::Time;
00088   using Teuchos::TimeMonitor;
00089 
00090   // Make some blocks to "orthogonalize."  Fill with random
00091   // data.  We only need X so that we can make clones (it knows
00092   // its data distribution).
00093   Array<RCP<MV> > V (numBlocks);
00094   for (int k = 0; k < numBlocks; ++k)
00095     {
00096       V[k] = MVT::Clone (*X, numCols);
00097       MVT::MvRandom (*V[k]);
00098     }
00099 
00100   // Make timers with informative labels
00101   RCP<Time> timer = TimeMonitor::getNewCounter ("Baseline for OrthoManager benchmark");
00102 
00103   // Baseline benchmark just copies data.  It's sort of a lower
00104   // bound proxy for the volume of data movement done by a real
00105   // OrthoManager.
00106   {
00107     TimeMonitor monitor (*timer);
00108     for (int trial = 0; trial < numTrials; ++trial)
00109       {
00110         for (int k = 0; k < numBlocks; ++k)
00111     {
00112       for (int j = 0; j < k; ++j)
00113         MVT::Assign (*V[j], *V[k]);
00114       MVT::Assign (*X, *V[k]);
00115     }
00116       }
00117   }
00118       }
00119 
00151       static void
00152       benchmark (const Teuchos::RCP<OrthoManager<Scalar, MV> >& orthoMan,
00153      const std::string& orthoManName,
00154      const std::string& normalization,
00155      const Teuchos::RCP<const MV>& X,
00156      const int numCols,
00157      const int numBlocks,
00158      const int numTrials,
00159      const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00160      std::ostream& resultStream,
00161      const bool displayResultsCompactly=false)
00162       {
00163   using Teuchos::Array;
00164   using Teuchos::ArrayView;
00165   using Teuchos::RCP;
00166   using Teuchos::rcp;
00167   using Teuchos::Time;
00168   using Teuchos::TimeMonitor;
00169   using std::endl;
00170         
00171   TEST_FOR_EXCEPTION(orthoMan.is_null(), std::invalid_argument,
00172          "orthoMan is null");
00173   TEST_FOR_EXCEPTION(X.is_null(), std::invalid_argument,
00174          "X is null");
00175   TEST_FOR_EXCEPTION(numCols < 1, std::invalid_argument, 
00176          "numCols = " << numCols << " < 1");
00177   TEST_FOR_EXCEPTION(numBlocks < 1, std::invalid_argument, 
00178          "numBlocks = " << numBlocks << " < 1");
00179   TEST_FOR_EXCEPTION(numTrials < 1, std::invalid_argument, 
00180          "numTrials = " << numTrials << " < 1");
00181   // Debug output stream
00182   std::ostream& debugOut = outMan->stream(Debug);
00183 
00184   // If you like, you can add the "baseline" as an approximate
00185   // lower bound for orthogonalization performance.  It may be
00186   // useful as a sanity check to make sure that your
00187   // orthogonalizations are really computing something, though
00188   // testing accuracy can help with that too.
00189   //
00190   //baseline (X, numCols, numBlocks, numTrials);
00191 
00192   // Make space to put the projection and normalization
00193   // coefficients.
00194   Array<RCP<mat_type> > C (numBlocks);
00195   for (int k = 0; k < numBlocks; ++k)
00196     C[k] = rcp (new mat_type (numCols, numCols));
00197   RCP<mat_type> B (new mat_type (numCols, numCols));
00198 
00199   // Make some blocks to orthogonalize.  Fill with random data.
00200   // We won't be orthogonalizing X, or even modifying X.  We
00201   // only need X so that we can make clones (since X knows its
00202   // data distribution).
00203   Array<RCP<MV> > V (numBlocks);
00204   for (int k = 0; k < numBlocks; ++k)
00205     {
00206       V[k] = MVT::Clone (*X, numCols);
00207       MVT::MvRandom (*V[k]);
00208     }
00209 
00210   // Make timers with informative labels.  We time an additional
00211   // first run to measure the startup costs, if any, of the
00212   // OrthoManager instance.
00213   RCP<Time> firstRunTimer;
00214   {
00215     std::ostringstream os;
00216     os << "OrthoManager: " << orthoManName << " first run";
00217     firstRunTimer = TimeMonitor::getNewCounter (os.str());
00218   }
00219   RCP<Time> timer;
00220   {
00221     std::ostringstream os;
00222     os << "OrthoManager: " << orthoManName << " total over " 
00223        << numTrials << " trials (excluding first run above)";
00224     timer = TimeMonitor::getNewCounter (os.str());
00225   }
00226   // The first run lets us measure the startup costs, if any, of
00227   // the OrthoManager instance, without these costs influencing
00228   // the following timing runs.
00229   {
00230     TimeMonitor monitor (*firstRunTimer);
00231     {
00232       (void) orthoMan->normalize (*V[0], B);
00233       for (int k = 1; k < numBlocks; ++k)
00234         {
00235     // k is the number of elements in the ArrayView.  We
00236     // have to assign first to an ArrayView-of-RCP-of-MV,
00237     // rather than to an ArrayView-of-RCP-of-const-MV, since
00238     // the latter requires a reinterpret cast.  Don't you
00239     // love C++ type inference?
00240     ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
00241     ArrayView<RCP<const MV> > V_0k = 
00242       Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
00243     (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
00244         }
00245     }
00246     // "Test" that the trial run actually orthogonalized
00247     // correctly.  Results are printed to the OutputManager's
00248     // Belos::Debug output stream, so depending on the
00249     // OutputManager's chosen verbosity level, you may or may
00250     // not see the results of the test.
00251     //
00252     // NOTE (mfh 22 Jan 2011) For now, these results have to be
00253     // inspected visually.  We should add a simple automatic
00254     // test.
00255     debugOut << "Orthogonality of V[0:" << (numBlocks-1) 
00256        << "]:" << endl;
00257     for (int k = 0; k < numBlocks; ++k)
00258       {
00259         // Orthogonality of each block
00260         debugOut << "For block V[" << k << "]:" << endl;
00261         debugOut << "  ||<V[" << k << "], V[" << k << "]> - I|| = " 
00262            << orthoMan->orthonormError(*V[k]) << endl;
00263         // Relative orthogonality with the previous blocks
00264         for (int j = 0; j < k; ++j)
00265     debugOut << "  ||< V[" << j << "], V[" << k << "] >|| = " 
00266        << orthoMan->orthogError(*V[j], *V[k]) << endl;
00267       }
00268   }
00269   
00270   // Run the benchmark for numTrials trials.  Time all trials as
00271   // a single run.
00272   {
00273     TimeMonitor monitor (*timer);
00274     
00275     for (int trial = 0; trial < numTrials; ++trial)
00276       {
00277         (void) orthoMan->normalize (*V[0], B);
00278         for (int k = 1; k < numBlocks; ++k)
00279     {
00280       ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
00281       ArrayView<RCP<const MV> > V_0k = 
00282         Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
00283       (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
00284     }
00285       }
00286   }
00287 
00288   // Report timing results.
00289   if (displayResultsCompactly)
00290     {
00291       // The "compact" format is suitable for automatic parsing,
00292       // using a CSV (Comma-Delimited Values) parser.  The first
00293       // "comment" line may be parsed to extract column
00294       // ("field") labels; the second line contains the actual
00295       // data, in ASCII comma-delimited format.
00296       using std::endl;
00297       resultStream << "#orthoManName"
00298        << ",normalization"
00299        << ",numRows"
00300        << ",numCols"
00301        << ",numBlocks"
00302        << ",firstRunTimeInSeconds"
00303        << ",timeInSeconds"
00304        << ",numTrials" 
00305        << endl;
00306       resultStream << orthoManName 
00307        << "," << (orthoManName=="Simple" ? normalization : "N/A")
00308        << "," << MVT::GetVecLength(*X)
00309        << "," << numCols
00310        << "," << numBlocks
00311        << "," << firstRunTimer->totalElapsedTime()
00312        << "," << timer->totalElapsedTime()
00313        << "," << numTrials
00314        << endl;
00315     }
00316   else
00317     TimeMonitor::summarize (resultStream);
00318       }
00319     };
00320 
00324     template< class Scalar, class MV >
00325     class OrthoManagerTester {
00326     private:
00327       typedef typename Teuchos::Array<Teuchos::RCP<MV> >::size_type size_type;
00328 
00329     public:
00330       typedef Scalar scalar_type;
00331       typedef Teuchos::ScalarTraits<scalar_type> SCT;
00332       typedef typename SCT::magnitudeType magnitude_type;
00333       typedef Teuchos::ScalarTraits<magnitude_type> SMT;
00334       typedef MultiVecTraits<scalar_type, MV> MVT;
00335       typedef Teuchos::SerialDenseMatrix<int, scalar_type> mat_type;
00336 
00353       static int
00354       runTests (const Teuchos::RCP<OrthoManager<Scalar, MV> >& OM,
00355     const bool isRankRevealing,
00356     const Teuchos::RCP<MV>& S,
00357     const int sizeX1,
00358     const int sizeX2,
00359     const Teuchos::RCP<OutputManager<Scalar> >& MyOM)
00360       {
00361   using Teuchos::Array;
00362   using Teuchos::null;
00363   using Teuchos::RCP;
00364   using Teuchos::rcp;
00365   using Teuchos::rcp_dynamic_cast;
00366   using Teuchos::tuple;
00367 
00368   // Number of tests that have failed thus far.
00369   int numFailed = 0; 
00370 
00371   // Relative tolerance against which all tests are performed.
00372   const magnitude_type TOL = 1.0e-12;
00373   // Absolute tolerance constant
00374   //const magnitude_type ATOL = 10;
00375 
00376   const scalar_type ZERO = SCT::zero();
00377   const scalar_type ONE = SCT::one();
00378 
00379   // Debug output stream
00380   std::ostream& debugOut = MyOM->stream(Debug);
00381 
00382   // Number of columns in the input "prototype" multivector S.
00383   const int sizeS = MVT::GetNumberVecs (*S);
00384 
00385   // Create multivectors X1 and X2, using the same map as multivector
00386   // S.  Then, test orthogonalizing X2 against X1.  After doing so, X1
00387   // and X2 should each be M-orthonormal, and should be mutually
00388   // M-orthogonal.
00389   debugOut << "Generating X1,X2 for testing... ";
00390   RCP< MV > X1 = MVT::Clone (*S, sizeX1);
00391   RCP< MV > X2 = MVT::Clone (*S, sizeX2);
00392   debugOut << "done." << endl;
00393   {
00394     magnitude_type err;
00395 
00396     //
00397     // Fill X1 with random values, and test the normalization error.
00398     //
00399     debugOut << "Filling X1 with random values... ";
00400     MVT::MvRandom(*X1);
00401     debugOut << "done." << endl
00402        << "Calling normalize() on X1... ";
00403     // The Anasazi and Belos OrthoManager interfaces differ.
00404     // For example, Anasazi's normalize() method accepts either
00405     // one or two arguments, whereas Belos' normalize() requires
00406     // two arguments.
00407     const int initialX1Rank = OM->normalize(*X1, Teuchos::null);
00408     TEST_FOR_EXCEPTION(initialX1Rank != sizeX1, 
00409            std::runtime_error, 
00410            "normalize(X1) returned rank "
00411            << initialX1Rank << " from " << sizeX1
00412            << " vectors. Cannot continue.");
00413     debugOut << "done." << endl 
00414        << "Calling orthonormError() on X1... ";
00415     err = OM->orthonormError(*X1);
00416     TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
00417            "After normalize(X1), orthonormError(X1) = " 
00418            << err << " > TOL = " << TOL);
00419     debugOut << "done: ||<X1,X1> - I|| = " << err << endl;
00420 
00421     //
00422     // Fill X2 with random values, project against X1 and normalize,
00423     // and test the orthogonalization error.
00424     //
00425     debugOut << "Filling X2 with random values... ";
00426     MVT::MvRandom(*X2);
00427     debugOut << "done." << endl
00428        << "Calling projectAndNormalize(X2, C, B, tuple(X1))... "
00429        << std::flush;
00430     // The projectAndNormalize() interface also differs between 
00431     // Anasazi and Belos.  Anasazi's projectAndNormalize() puts 
00432     // the multivector and the array of multivectors first, and
00433     // the (array of) SerialDenseMatrix arguments (which are 
00434     // optional) afterwards.  Belos puts the (array of) 
00435     // SerialDenseMatrix arguments in the middle, and they are 
00436     // not optional.
00437     int initialX2Rank;
00438     {
00439       Array<RCP<mat_type> > C (1);
00440       RCP<mat_type> B = Teuchos::null;
00441       initialX2Rank = 
00442         OM->projectAndNormalize (*X2, C, B, tuple<RCP<const MV> >(X1));
00443     }
00444     TEST_FOR_EXCEPTION(initialX2Rank != sizeX2, 
00445            std::runtime_error, 
00446            "projectAndNormalize(X2,X1) returned rank " 
00447            << initialX2Rank << " from " << sizeX2 
00448            << " vectors. Cannot continue.");
00449     debugOut << "done." << endl
00450        << "Calling orthonormError() on X2... ";
00451     err = OM->orthonormError (*X2);
00452     TEST_FOR_EXCEPTION(err > TOL,
00453            std::runtime_error,
00454            "projectAndNormalize(X2,X1) did not meet tolerance: "
00455            "orthonormError(X2) = " << err << " > TOL = " << TOL);
00456     debugOut << "done: || <X2,X2> - I || = " << err << endl
00457        << "Calling orthogError(X2, X1)... ";
00458     err = OM->orthogError (*X2, *X1);
00459     TEST_FOR_EXCEPTION(err > TOL,
00460            std::runtime_error,
00461            "projectAndNormalize(X2,X1) did not meet tolerance: "
00462            "orthogError(X2,X1) = " << err << " > TOL = " << TOL);
00463     debugOut << "done: || <X2,X1> || = " << err << endl;
00464   }
00465 
00466 
00467   //
00468   // If OM is an OutOfPlaceNormalizerMixin, exercise the
00469   // out-of-place normalization routines.
00470   //
00471   typedef Belos::OutOfPlaceNormalizerMixin<Scalar, MV> mixin_type;
00472   RCP<mixin_type> tsqr = rcp_dynamic_cast<mixin_type>(OM);
00473   if (! tsqr.is_null())
00474     {
00475       magnitude_type err;
00476       debugOut << endl 
00477          << "=== OutOfPlaceNormalizerMixin tests ===" 
00478          << endl << endl;
00479       //
00480       // Fill X1_in with random values, and test the normalization
00481       // error with normalizeOutOfPlace().
00482       //
00483       // Don't overwrite X1, else you'll mess up the tests that
00484       // follow!
00485       //
00486       RCP<MV> X1_in = MVT::CloneCopy (*X1);
00487       debugOut << "Filling X1_in with random values... ";
00488       MVT::MvRandom(*X1_in);
00489       debugOut << "done." << endl;
00490       debugOut << "Filling X1_out with different random values...";
00491       RCP<MV> X1_out = MVT::Clone(*X1_in, MVT::GetNumberVecs(*X1_in));
00492       MVT::MvRandom(*X1_out);
00493       debugOut << "done." << endl
00494          << "Calling normalizeOutOfPlace(*X1_in, *X1_out, null)... ";
00495       const int initialX1Rank = 
00496         tsqr->normalizeOutOfPlace(*X1_in, *X1_out, Teuchos::null);
00497       TEST_FOR_EXCEPTION(initialX1Rank != sizeX1, std::runtime_error,
00498              "normalizeOutOfPlace(*X1_in, *X1_out, null) "
00499              "returned rank " << initialX1Rank << " from " 
00500              << sizeX1 << " vectors. Cannot continue.");
00501       debugOut << "done." << endl 
00502          << "Calling orthonormError() on X1_out... ";
00503       err = OM->orthonormError(*X1_out);
00504       TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
00505              "After calling normalizeOutOfPlace(*X1_in, "
00506              "*X1_out, null), orthonormError(X1) = " 
00507              << err << " > TOL = " << TOL);
00508       debugOut << "done: ||<X1_out,X1_out> - I|| = " << err << endl;
00509 
00510       //
00511       // Fill X2_in with random values, project against X1_out
00512       // and normalize via projectAndNormalizeOutOfPlace(), and
00513       // test the orthogonalization error.
00514       //
00515       // Don't overwrite X2, else you'll mess up the tests that
00516       // follow!
00517       //
00518       RCP<MV> X2_in = MVT::CloneCopy (*X2);
00519       debugOut << "Filling X2_in with random values... ";
00520       MVT::MvRandom(*X2_in);
00521       debugOut << "done." << endl
00522          << "Filling X2_out with different random values...";
00523       RCP<MV> X2_out = MVT::Clone(*X2_in, MVT::GetNumberVecs(*X2_in));
00524       MVT::MvRandom(*X2_out);
00525       debugOut << "done." << endl
00526          << "Calling projectAndNormalizeOutOfPlace(X2_in, X2_out, "
00527          << "C, B, X1_out)...";
00528       int initialX2Rank;
00529       {
00530         Array<RCP<mat_type> > C (1);
00531         RCP<mat_type> B = Teuchos::null;
00532         initialX2Rank = 
00533     tsqr->projectAndNormalizeOutOfPlace (*X2_in, *X2_out, C, B, 
00534                  tuple<RCP<const MV> >(X1_out));
00535       }
00536       TEST_FOR_EXCEPTION(initialX2Rank != sizeX2, 
00537              std::runtime_error, 
00538              "projectAndNormalizeOutOfPlace(*X2_in, "
00539              "*X2_out, C, B, tuple(X1_out)) returned rank "
00540              << initialX2Rank << " from " << sizeX2 
00541              << " vectors. Cannot continue.");
00542       debugOut << "done." << endl
00543          << "Calling orthonormError() on X2_out... ";
00544       err = OM->orthonormError (*X2_out);
00545       TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
00546              "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
00547              "C, B, tuple(X1_out)) did not meet tolerance: "
00548              "orthonormError(X2_out) = " 
00549              << err << " > TOL = " << TOL);
00550       debugOut << "done: || <X2_out,X2_out> - I || = " << err << endl
00551          << "Calling orthogError(X2_out, X1_out)... ";
00552       err = OM->orthogError (*X2_out, *X1_out);
00553       TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
00554              "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
00555              "C, B, tuple(X1_out)) did not meet tolerance: "
00556              "orthogError(X2_out, X1_out) = " 
00557              << err << " > TOL = " << TOL);
00558       debugOut << "done: || <X2_out,X1_out> || = " << err << endl;
00559       debugOut << endl 
00560          << "=== Done with OutOfPlaceNormalizerMixin tests ===" 
00561          << endl << endl;
00562     }
00563 
00564   {
00565     //
00566     // Test project() on a random multivector S, by projecting S
00567     // against various combinations of X1 and X2.
00568     //
00569     MVT::MvRandom(*S);
00570 
00571     debugOut << "Testing project() by projecting a random multivector S "
00572       "against various combinations of X1 and X2 " << endl;
00573     const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
00574     numFailed += thisNumFailed;
00575     if (thisNumFailed > 0)
00576       debugOut << "  *** " << thisNumFailed 
00577          << (thisNumFailed > 1 ? " tests" : " test") 
00578          << " failed." << endl;
00579   }
00580 
00581   if (isRankRevealing)
00582     {
00583       // run a X1,Y2 range multivector against P_{X1,X1} P_{Y2,Y2}
00584       // note, this is allowed under the restrictions on project(), 
00585       // because <X1,Y2> = 0
00586       // also, <Y2,Y2> = I, but <X1,X1> != I, so biOrtho must be set to false
00587       // it should require randomization, as 
00588       // P_{X1,X1} P_{Y2,Y2} (X1*C1 + Y2*C2) = P_{X1,X1} X1*C1 = 0
00589       mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
00590       C1.random();
00591       C2.random();
00592       // S := X1*C1
00593       MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
00594       // S := S + X2*C2
00595       MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
00596 
00597       debugOut << "Testing project() by projecting [X1 X2]-range multivector "
00598         "against P_X1 P_X2 " << endl;
00599       const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
00600       numFailed += thisNumFailed;
00601       if (thisNumFailed > 0)
00602         debugOut << "  *** " << thisNumFailed 
00603            << (thisNumFailed > 1 ? " tests" : " test") 
00604            << " failed." << endl;
00605     }
00606 
00607   // This test is only distinct from the rank-1 multivector test
00608   // (below) if S has at least 3 columns.
00609   if (isRankRevealing && sizeS > 2) 
00610     {
00611       MVT::MvRandom(*S);
00612       RCP<MV> mid = MVT::Clone(*S,1);
00613       mat_type c(sizeS,1);
00614       MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
00615       std::vector<int> ind(1); 
00616       ind[0] = sizeS-1;
00617       MVT::SetBlock(*mid,ind,*S);
00618 
00619       debugOut << "Testing normalize() on a rank-deficient multivector " << endl;
00620       const int thisNumFailed = testNormalize(OM,S,MyOM);
00621       numFailed += thisNumFailed;
00622       if (thisNumFailed > 0)
00623         debugOut << "  *** " << thisNumFailed 
00624            << (thisNumFailed > 1 ? " tests" : " test") 
00625            << " failed." << endl;
00626     }
00627 
00628   // This test will only exercise rank deficiency if S has at least 2
00629   // columns.
00630   if (isRankRevealing && sizeS > 1) 
00631     {
00632       // rank-1
00633       RCP<MV> one = MVT::Clone(*S,1);
00634       MVT::MvRandom(*one);
00635       // put multiple of column 0 in columns 0:sizeS-1
00636       for (int i=0; i<sizeS; i++) 
00637         {
00638     std::vector<int> ind(1); 
00639     ind[0] = i;
00640     RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
00641     MVT::MvAddMv(SCT::random(),*one,ZERO,*one,*Si);
00642         }
00643       debugOut << "Testing normalize() on a rank-1 multivector " << endl;
00644       const int thisNumFailed = testNormalize(OM,S,MyOM);
00645       numFailed += thisNumFailed;
00646       if (thisNumFailed > 0)
00647         debugOut << "  *** " << thisNumFailed 
00648            << (thisNumFailed > 1 ? " tests" : " test") 
00649            << " failed." << endl;
00650     }
00651 
00652   {
00653     std::vector<int> ind(1); 
00654     MVT::MvRandom(*S);
00655 
00656     debugOut << "Testing projectAndNormalize() on a random multivector " << endl;
00657     const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
00658     numFailed += thisNumFailed;
00659     if (thisNumFailed > 0)
00660       debugOut << "  *** " << thisNumFailed 
00661          << (thisNumFailed > 1 ? " tests" : " test") 
00662          << " failed." << endl;
00663   }
00664 
00665   if (isRankRevealing)
00666     {
00667       // run a X1,X2 range multivector against P_X1 P_X2
00668       // this is allowed as <X1,X2> == 0
00669       // it should require randomization, as 
00670       // P_X1 P_X2 (X1*C1 + X2*C2) = P_X1 X1*C1 = 0
00671       // and 
00672       // P_X2 P_X1 (X2*C2 + X1*C1) = P_X2 X2*C2 = 0
00673       mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
00674       C1.random();
00675       C2.random();
00676       MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
00677       MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
00678 
00679       debugOut << "Testing projectAndNormalize() by projecting [X1 X2]-range "
00680         "multivector against P_X1 P_X2 " << endl;
00681       const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
00682       numFailed += thisNumFailed;
00683       if (thisNumFailed > 0)
00684         debugOut << "  *** " << thisNumFailed 
00685            << (thisNumFailed > 1 ? " tests" : " test") 
00686            << " failed." << endl;
00687     }
00688 
00689   // This test is only distinct from the rank-1 multivector test
00690   // (below) if S has at least 3 columns.
00691   if (isRankRevealing && sizeS > 2) 
00692     {
00693       MVT::MvRandom(*S);
00694       RCP<MV> mid = MVT::Clone(*S,1);
00695       mat_type c(sizeS,1);
00696       MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
00697       std::vector<int> ind(1); 
00698       ind[0] = sizeS-1;
00699       MVT::SetBlock(*mid,ind,*S);
00700 
00701       debugOut << "Testing projectAndNormalize() on a rank-deficient "
00702         "multivector " << endl;
00703       const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
00704       numFailed += thisNumFailed;
00705       if (thisNumFailed > 0)
00706         debugOut << "  *** " << thisNumFailed 
00707            << (thisNumFailed > 1 ? " tests" : " test") 
00708            << " failed." << endl;
00709     }
00710 
00711   // This test will only exercise rank deficiency if S has at least 2
00712   // columns.
00713   if (isRankRevealing && sizeS > 1) 
00714     {
00715       // rank-1
00716       RCP<MV> one = MVT::Clone(*S,1);
00717       MVT::MvRandom(*one);
00718       // Put a multiple of column 0 in columns 0:sizeS-1.
00719       for (int i=0; i<sizeS; i++) 
00720         {
00721     std::vector<int> ind(1); 
00722     ind[0] = i;
00723     RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
00724     MVT::MvAddMv(SCT::random(),*one,ZERO,*one,*Si);
00725         }
00726       debugOut << "Testing projectAndNormalize() on a rank-1 multivector " << endl;
00727       bool constantStride = true;
00728       if (! MVT::HasConstantStride(*S))
00729         {
00730     debugOut << "-- S does not have constant stride" << endl;
00731     constantStride = false;
00732         }
00733       if (! MVT::HasConstantStride(*X1))
00734         {
00735     debugOut << "-- X1 does not have constant stride" << endl;
00736     constantStride = false;
00737         }
00738       if (! MVT::HasConstantStride(*X2))
00739         {
00740     debugOut << "-- X2 does not have constant stride" << endl;
00741     constantStride = false;
00742         }
00743       if (! constantStride)
00744         debugOut << "-- Skipping this test, since TSQR does not work on "
00745     "multivectors with nonconstant stride" << endl;
00746       else
00747         {
00748     const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
00749     numFailed += thisNumFailed;
00750     if (thisNumFailed > 0)
00751       debugOut << "  *** " << thisNumFailed 
00752          << (thisNumFailed > 1 ? " tests" : " test") 
00753          << " failed." << endl;
00754         }
00755     }
00756 
00757   if (numFailed != 0)
00758     MyOM->stream(Errors) << numFailed << " total test failures." << endl;
00759 
00760   return numFailed;
00761       }
00762 
00763     private:
00764 
00769       static magnitude_type
00770       MVDiff (const MV& X, const MV& Y)
00771       {
00772   using Teuchos::RCP;
00773 
00774   const scalar_type ONE = SCT::one();
00775   const int numCols = MVT::GetNumberVecs(X);
00776   TEST_FOR_EXCEPTION( (MVT::GetNumberVecs(Y) != numCols),
00777           std::logic_error,
00778           "MVDiff: X and Y should have the same number of columns."
00779           "  X has " << numCols << " column(s) and Y has " 
00780           << MVT::GetNumberVecs(Y) << " columns." );
00781   // Resid := X
00782   RCP< MV > Resid = MVT::CloneCopy(X);
00783   // Resid := Resid - Y
00784   MVT::MvAddMv (-ONE, Y, ONE, *Resid, *Resid);
00785 
00786   return frobeniusNorm (*Resid);
00787       }
00788 
00789 
00793       static magnitude_type
00794       frobeniusNorm (const MV& X)
00795       {
00796   const scalar_type ONE = SCT::one();
00797   const int numCols = MVT::GetNumberVecs(X);
00798   mat_type C (numCols, numCols);
00799 
00800   // $C := X^* X$
00801   MVT::MvTransMv (ONE, X, X, C);
00802 
00803   magnitude_type err (0);
00804   for (int i = 0; i < numCols; ++i)
00805     err += SCT::magnitude (C(i,i));
00806 
00807   return SCT::magnitude (SCT::squareroot (err));
00808       }
00809 
00810 
00811       static int
00812       testProjectAndNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, 
00813              const Teuchos::RCP< const MV >& S, 
00814              const Teuchos::RCP< const MV >& X1, 
00815              const Teuchos::RCP< const MV >& X2,
00816              const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
00817       {
00818   return testProjectAndNormalizeNew (OM, S, X1, X2, MyOM);
00819       }
00820 
00825       static int 
00826       testProjectAndNormalizeOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
00827           const Teuchos::RCP< const MV >& S, 
00828           const Teuchos::RCP< const MV >& X1, 
00829           const Teuchos::RCP< const MV >& X2,
00830           const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
00831       {
00832   using Teuchos::Array;
00833   using Teuchos::null;
00834   using Teuchos::RCP;
00835   using Teuchos::rcp;
00836   using Teuchos::tuple;
00837 
00838   const scalar_type ONE = SCT::one();
00839   const magnitude_type ZERO = SCT::magnitude(SCT::zero());
00840 
00841   // Relative tolerance against which all tests are performed.
00842   const magnitude_type TOL = 1.0e-12;
00843   // Absolute tolerance constant
00844   const magnitude_type ATOL = 10;
00845 
00846   const int sizeS = MVT::GetNumberVecs(*S);
00847   const int sizeX1 = MVT::GetNumberVecs(*X1);
00848   const int sizeX2 = MVT::GetNumberVecs(*X2);
00849   int numerr = 0;
00850   std::ostringstream sout;
00851 
00852   //
00853   // output tests:
00854   //   <S_out,S_out> = I
00855   //   <S_out,X1> = 0
00856   //   <S_out,X2> = 0
00857   //   S_in = S_out B + X1 C1 + X2 C2
00858   // 
00859   // we will loop over an integer specifying the test combinations
00860   // the bit pattern for the different tests is listed in parenthesis
00861   //
00862   // for the projectors, test the following combinations:
00863   // none              (00)
00864   // P_X1              (01)
00865   // P_X2              (10)
00866   // P_X1 P_X2         (11)
00867   // P_X2 P_X1         (11)
00868   // the latter two should be tested to give the same answer
00869   //
00870   // for each of these, we should test with C1, C2 and B
00871   //
00872   // if hasM:
00873   // with and without MX1   (1--) 
00874   // with and without MX2  (1---) 
00875   // with and without MS  (1----) 
00876   //
00877   // as hasM controls the upper level bits, we need only run test cases 0-3 if hasM==false
00878   // otherwise, we run test cases 0-31
00879   //
00880 
00881   int numtests = 4;
00882 
00883   // test ortho error before orthonormalizing
00884   if (X1 != null) {
00885     magnitude_type err = OM->orthogError(*S,*X1);
00886     sout << "   || <S,X1> || before     : " << err << endl;
00887   }
00888   if (X2 != null) {
00889     magnitude_type err = OM->orthogError(*S,*X2);
00890     sout << "   || <S,X2> || before     : " << err << endl;
00891   }
00892 
00893   for (int t=0; t<numtests; t++) {
00894 
00895     Array< RCP< const MV > > theX;
00896     RCP<mat_type > B = rcp( new mat_type(sizeS,sizeS) );
00897     Array<RCP<mat_type > > C;
00898     if ( (t && 3) == 0 ) {
00899       // neither <X1,Y1> nor <X2,Y2>
00900       // C, theX and theY are already empty
00901     }
00902     else if ( (t && 3) == 1 ) {
00903       // X1
00904       theX = tuple(X1);
00905       C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
00906     }
00907     else if ( (t && 3) == 2 ) {
00908       // X2
00909       theX = tuple(X2);
00910       C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
00911     }
00912     else {
00913       // X1 and X2, and the reverse.
00914       theX = tuple(X1,X2);
00915       C = tuple( rcp(new mat_type(sizeX1,sizeS)), 
00916            rcp(new mat_type(sizeX2,sizeS)) );
00917     }
00918 
00919     // We wrap up all the OrthoManager calls in a try-catch
00920     // block, in order to check whether any of the methods throw
00921     // an exception.  For the tests we perform, every thrown
00922     // exception is a failure.
00923     try {
00924       // call routine
00925       // if (t && 3) == 3, {
00926       //    call with reversed input: X2 X1
00927       // }
00928       // test all outputs for correctness
00929       // test all outputs for equivalence
00930 
00931       // here is where the outputs go
00932       Array<RCP<MV> > S_outs;
00933       Array<Array<RCP<mat_type > > > C_outs;
00934       Array<RCP<mat_type > > B_outs;
00935       RCP<MV> Scopy;
00936       Array<int> ret_out;
00937 
00938       // copies of S,MS
00939       Scopy = MVT::CloneCopy(*S);
00940       // randomize this data, it should be overwritten
00941       B->random();
00942       for (size_type i=0; i<C.size(); i++) {
00943         C[i]->random();
00944       }
00945       // Run test.  Since S was specified by the caller and
00946       // Scopy is a copy of S, we don't know what rank to expect
00947       // here -- though we do require that S have rank at least
00948       // one.
00949       //
00950       // Note that Anasazi and Belos differ, among other places, 
00951       // in the order of arguments to projectAndNormalize().
00952       int ret = OM->projectAndNormalize(*Scopy,C,B,theX);
00953       sout << "projectAndNormalize() returned rank " << ret << endl;
00954       if (ret == 0) {
00955         sout << "  *** Error: returned rank is zero, cannot continue tests" << endl;
00956         numerr++;
00957         break;
00958       }
00959       ret_out.push_back(ret);
00960       // projectAndNormalize() is only required to return a 
00961       // basis of rank "ret"
00962       // this is what we will test:
00963       //   the first "ret" columns in Scopy
00964       //   the first "ret" rows in B
00965       // save just the parts that we want
00966       // we allocate S and MS for each test, so we can save these as views
00967       // however, save copies of the C and B
00968       if (ret < sizeS) {
00969         std::vector<int> ind(ret);
00970         for (int i=0; i<ret; i++) {
00971     ind[i] = i;
00972         }
00973         S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
00974         B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
00975       }
00976       else {
00977         S_outs.push_back( Scopy );
00978         B_outs.push_back( rcp( new mat_type(*B) ) );
00979       }
00980       C_outs.push_back( Array<RCP<mat_type > >(0) );
00981       if (C.size() > 0) {
00982         C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
00983       }
00984       if (C.size() > 1) {
00985         C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
00986       }
00987 
00988       // do we run the reversed input?
00989       if ( (t && 3) == 3 ) {
00990         // copies of S,MS
00991         Scopy = MVT::CloneCopy(*S);
00992 
00993         // Fill the B and C[i] matrices with random data.  The
00994         // data will be overwritten by projectAndNormalize().
00995         // Filling these matrices here is only to catch some
00996         // bugs in projectAndNormalize().
00997         B->random();
00998         for (size_type i=0; i<C.size(); i++) {
00999     C[i]->random();
01000         }
01001         // flip the inputs
01002         theX = tuple( theX[1], theX[0] );
01003         // Run test.
01004         // Note that Anasazi and Belos differ, among other places, 
01005         // in the order of arguments to projectAndNormalize().
01006         ret = OM->projectAndNormalize(*Scopy,C,B,theX);
01007         sout << "projectAndNormalize() returned rank " << ret << endl;
01008         if (ret == 0) {
01009     sout << "  *** Error: returned rank is zero, cannot continue tests" << endl;
01010     numerr++;
01011     break;
01012         }
01013         ret_out.push_back(ret);
01014         // projectAndNormalize() is only required to return a 
01015         // basis of rank "ret"
01016         // this is what we will test:
01017         //   the first "ret" columns in Scopy
01018         //   the first "ret" rows in B
01019         // save just the parts that we want
01020         // we allocate S and MS for each test, so we can save these as views
01021         // however, save copies of the C and B
01022         if (ret < sizeS) {
01023     std::vector<int> ind(ret);
01024     for (int i=0; i<ret; i++) {
01025       ind[i] = i;
01026     }
01027     S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
01028     B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
01029         }
01030         else {
01031     S_outs.push_back( Scopy );
01032     B_outs.push_back( rcp( new mat_type(*B) ) );
01033         }
01034         C_outs.push_back( Array<RCP<mat_type > >() );
01035         // reverse the Cs to compensate for the reverse projectors
01036         C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
01037         C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
01038         // flip the inputs back
01039         theX = tuple( theX[1], theX[0] );
01040       }
01041 
01042 
01043       // test all outputs for correctness
01044       for (size_type o=0; o<S_outs.size(); o++) {
01045         // S^T M S == I
01046         {
01047     magnitude_type err = OM->orthonormError(*S_outs[o]);
01048     if (err > TOL) {
01049       sout << endl
01050            << "  *** Test (number " << (t+1) << " of " << numtests 
01051            << " total tests) failed: Tolerance exceeded!  Error = "
01052            << err << " > TOL = " << TOL << "." 
01053            << endl << endl;
01054       numerr++;
01055     }
01056     sout << "   || <S,S> - I || after  : " << err << endl;
01057         }
01058         // S_in = X1*C1 + C2*C2 + S_out*B
01059         {
01060     RCP<MV> tmp = MVT::Clone(*S,sizeS);
01061     MVT::MvTimesMatAddMv(ONE,*S_outs[o],*B_outs[o],ZERO,*tmp);
01062     if (C_outs[o].size() > 0) {
01063       MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
01064       if (C_outs[o].size() > 1) {
01065         MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
01066       }
01067     }
01068     magnitude_type err = MVDiff(*tmp,*S);
01069     if (err > ATOL*TOL) {
01070       sout << endl
01071            << "  *** Test (number " << (t+1) << " of " << numtests 
01072            << " total tests) failed: Tolerance exceeded!  Error = "
01073            << err << " > ATOL*TOL = " << (ATOL*TOL) << "."
01074            << endl << endl;
01075       numerr++;
01076     }
01077     sout << "  " << t << "|| S_in - X1*C1 - X2*C2 - S_out*B || : " << err << endl;
01078         }
01079         // <X1,S> == 0
01080         if (theX.size() > 0 && theX[0] != null) {
01081     magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
01082     if (err > TOL) {
01083       sout << endl
01084            << "  *** Test (number " << (t+1) << " of " << numtests 
01085            << " total tests) failed: Tolerance exceeded!  Error = "
01086            << err << " > TOL = " << TOL << "."
01087            << endl << endl;
01088       numerr++;
01089     }
01090     sout << "  " << t << "|| <X[0],S> || after      : " << err << endl;
01091         }
01092         // <X2,S> == 0
01093         if (theX.size() > 1 && theX[1] != null) {
01094     magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
01095     if (err > TOL) {
01096       sout << endl
01097            << "  *** Test (number " << (t+1) << " of " << numtests 
01098            << " total tests) failed: Tolerance exceeded!  Error = "
01099            << err << " > TOL = " << TOL << "."
01100            << endl << endl;
01101       numerr++;
01102     }
01103     sout << "  " << t << "|| <X[1],S> || after      : " << err << endl;
01104         }
01105       }
01106     }
01107     catch (Belos::OrthoError& e) {
01108       sout << "  *** Error: OrthoManager threw exception: " << e.what() << endl;
01109       numerr++;
01110     }
01111 
01112   } // test for
01113 
01114   // NOTE (mfh 05 Nov 2010) Since Belos::MsgType is an enum,
01115   // doing bitwise logical computations on Belos::MsgType values
01116   // (such as "Debug | Errors") and passing the result into
01117   // MyOM->stream() confuses the compiler.  As a result, we have
01118   // to do some type casts to make it work.
01119   const int msgType = (numerr > 0) ? 
01120     (static_cast<int>(Debug) | static_cast<int>(Errors)) :
01121     static_cast<int>(Debug);
01122 
01123   // We report debug-level messages always.  We also report
01124   // errors if at least one test failed.
01125   MyOM->stream(static_cast< MsgType >(msgType)) << sout.str() << endl;
01126   return numerr;
01127       }
01128 
01133       static int 
01134       testNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM, 
01135          const Teuchos::RCP< const MV >& S,
01136          const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
01137       {
01138   using Teuchos::Array;
01139   using Teuchos::RCP;
01140   using Teuchos::rcp;
01141   using Teuchos::tuple;
01142 
01143   const scalar_type ONE = SCT::one();
01144   std::ostringstream sout;
01145   // Total number of failed tests in this call of this routine.
01146   int numerr = 0;
01147 
01148   // Relative tolerance against which all tests are performed.
01149   // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
01150   // The following bounds hold for all $m \times n$ matrices $A$:
01151   // \[
01152   // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
01153   // \]
01154   // where $r$ is the (column) rank of $A$.  We bound this above
01155   // by the number of columns in $A$.  
01156   //
01157   // An accurate normalization in the Euclidean norm of a matrix
01158   // $A$ with at least as many rows m as columns n, should
01159   // produce orthogonality $\|Q^* Q - I\|_2$ less than a factor
01160   // of machine precision times a low-order polynomial in m and
01161   // n, and residual $\|A - Q B\|_2$ (where $A = Q B$ is the
01162   // computed normalization) less than that bound times the norm
01163   // of $A$.
01164   //
01165   // Since we are measuring both of these quantitites in the
01166   // Frobenius norm instead, we should scale this bound by
01167   // $\sqrt{n}$.
01168 
01169   const int numRows = MVT::GetVecLength(*S);
01170   const int numCols = MVT::GetNumberVecs(*S);
01171   const int sizeS = MVT::GetNumberVecs(*S);
01172 
01173   // A good heuristic is to scale the bound by the square root
01174   // of the number of floating-point operations.  One could
01175   // perhaps support this theoretically, since we are using
01176   // uniform random test problems.
01177   const magnitude_type fudgeFactor = 
01178     SMT::squareroot(magnitude_type(numRows) * 
01179         magnitude_type(numCols) * 
01180         magnitude_type(numCols));
01181   const magnitude_type TOL = SMT::eps() * fudgeFactor *
01182     SMT::squareroot(magnitude_type(numCols));
01183 
01184   // Absolute tolerance scaling: the Frobenius norm of the test
01185   // matrix S.  TOL*ATOL is the absolute tolerance for the
01186   // residual $\|A - Q*B\|_F$.
01187   const magnitude_type ATOL = frobeniusNorm (*S);
01188 
01189   sout << "The test matrix S has Frobenius norm " << ATOL 
01190        << ", and the relative error tolerance is TOL = " 
01191        << TOL << "." << endl;
01192 
01193   const int numtests = 1;
01194   for (int t = 0; t < numtests; ++t) {
01195 
01196     try {
01197       // call routine
01198       // test all outputs for correctness
01199 
01200       // S_copy gets a copy of S; we normalize in place, so we
01201       // need a copy to check whether the normalization
01202       // succeeded.
01203       RCP< MV > S_copy = MVT::CloneCopy (*S);
01204 
01205       // Matrix of coefficients from the normalization.
01206       RCP< mat_type > B (new mat_type (sizeS, sizeS));
01207       // The contents of B will be overwritten, but fill with
01208       // random data just to make sure that the normalization
01209       // operated on all the elements of B on which it should
01210       // operate.
01211       B->random();
01212 
01213       const int reportedRank = OM->normalize (*S_copy, B);
01214       sout << "normalize() returned rank " << reportedRank << endl;
01215       if (reportedRank == 0) {
01216         sout << "  *** Error: Cannot continue, since normalize() "
01217     "reports that S has rank 0" << endl;
01218         numerr++;
01219         break;
01220       }
01221       //
01222       // We don't know in this routine whether the input
01223       // multivector S has full rank; it is only required to
01224       // have nonzero rank.  Thus, we extract the first
01225       // reportedRank columns of S_copy and the first
01226       // reportedRank rows of B, and perform tests on them.
01227       //
01228 
01229       // Construct S_view, a view of the first reportedRank
01230       // columns of S_copy.
01231       std::vector<int> indices (reportedRank);
01232       for (int j = 0; j < reportedRank; ++j)
01233         indices[j] = j;
01234       RCP< MV > S_view = MVT::CloneViewNonConst (*S_copy, indices);
01235       // Construct B_top, a copy of the first reportedRank rows
01236       // of B.
01237       //
01238       // NOTE: We create this as a copy and not a view, because
01239       // otherwise it would not be safe with respect to RCPs.
01240       // This is because mat_type uses raw pointers
01241       // inside, so that a view would become invalid when B
01242       // would fall out of scope.
01243       RCP< mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, sizeS));
01244 
01245       // Check ||<S_view,S_view> - I||
01246       {
01247         const magnitude_type err = OM->orthonormError(*S_view);
01248         if (err > TOL) {
01249     sout << "  *** Error: Tolerance exceeded: err = " 
01250          << err << " > TOL = " << TOL << endl;
01251     numerr++;
01252         }
01253         sout << "   || <S,S> - I || after  : " << err << endl;
01254       }
01255       // Check the residual ||Residual|| = ||S_view * B_top -
01256       // S_orig||, where S_orig is a view of the first
01257       // reportedRank columns of S.
01258       {
01259         // Residual is allocated with reportedRank columns.  It
01260         // will contain the result of testing the residual error
01261         // of the normalization (i.e., $\|S - S_in*B\|$).  It
01262         // should have the dimensions of S.  Its initial value
01263         // is a copy of the first reportedRank columns of S.
01264         RCP< MV > Residual = MVT::CloneCopy (*S);
01265 
01266         // Residual := Residual - S_view * B_view
01267         MVT::MvTimesMatAddMv (-ONE, *S_view, *B_top, ONE, *Residual);
01268 
01269         // Compute ||Residual||
01270         const magnitude_type err = frobeniusNorm (*Residual);
01271         if (err > ATOL*TOL) {
01272     sout << "  *** Error: Tolerance exceeded: err = " 
01273          << err << " > ATOL*TOL = " << (ATOL*TOL) << endl;
01274     numerr++;
01275         }
01276         sout << "  " << t << "|| S - Q*B || : " << err << endl;
01277       }
01278     }
01279     catch (Belos::OrthoError& e) {
01280       sout << "  *** Error: the OrthoManager's normalize() method "
01281         "threw an exception: " << e.what() << endl;
01282       numerr++;
01283     }
01284 
01285   } // test for
01286 
01287   const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
01288   MyOM->stream(type) << sout.str();
01289   MyOM->stream(type) << endl;
01290 
01291   return numerr;
01292       }
01293 
01298       static int
01299       testProjectAndNormalizeNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, 
01300           const Teuchos::RCP< const MV >& S, 
01301           const Teuchos::RCP< const MV >& X1, 
01302           const Teuchos::RCP< const MV >& X2,
01303           const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
01304       {
01305   using Teuchos::Array;
01306   using Teuchos::null;
01307   using Teuchos::RCP;
01308   using Teuchos::rcp;
01309   using Teuchos::tuple;
01310 
01311   // We collect all the output in this string wrapper, and print
01312   // it at the end.
01313   std::ostringstream sout;
01314   // Total number of failed tests in this call of this routine.
01315   int numerr = 0;
01316 
01317   const int numRows = MVT::GetVecLength(*S);
01318   const int numCols = MVT::GetNumberVecs(*S);
01319   const int sizeS = MVT::GetNumberVecs(*S);
01320 
01321   // Relative tolerance against which all tests are performed.
01322   // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
01323   // The following bounds hold for all $m \times n$ matrices $A$:
01324   // \[
01325   // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
01326   // \]
01327   // where $r$ is the (column) rank of $A$.  We bound this above
01328   // by the number of columns in $A$.  
01329   //
01330   // Since we are measuring both of these quantitites in the
01331   // Frobenius norm instead, we scale all error tests by 
01332   // $\sqrt{n}$.
01333   //
01334   // A good heuristic is to scale the bound by the square root
01335   // of the number of floating-point operations.  One could
01336   // perhaps support this theoretically, since we are using
01337   // uniform random test problems.
01338   const magnitude_type fudgeFactor = 
01339     SMT::squareroot(magnitude_type(numRows) * 
01340         magnitude_type(numCols) * 
01341         magnitude_type(numCols));
01342   const magnitude_type TOL = SMT::eps() * fudgeFactor *
01343     SMT::squareroot(magnitude_type(numCols));
01344 
01345   // Absolute tolerance scaling: the Frobenius norm of the test
01346   // matrix S.  TOL*ATOL is the absolute tolerance for the
01347   // residual $\|A - Q*B\|_F$.
01348   const magnitude_type ATOL = frobeniusNorm (*S);
01349 
01350   sout << "-- The test matrix S has Frobenius norm " << ATOL 
01351        << ", and the relative error tolerance is TOL = " 
01352        << TOL << "." << endl;
01353 
01354   // Q will contain the result of projectAndNormalize() on S.
01355   RCP< MV > Q = MVT::CloneCopy(*S);
01356   // We use this for collecting the residual error components
01357   RCP< MV > Residual = MVT::CloneCopy(*S);
01358   // Number of elements in the X array of blocks against which
01359   // to project S.
01360   const int num_X = 2;
01361   Array< RCP< const MV > > X (num_X);
01362   X[0] = MVT::CloneCopy(*X1);
01363   X[1] = MVT::CloneCopy(*X2);
01364 
01365   // Coefficients for the normalization
01366   RCP< mat_type > B (new mat_type (sizeS, sizeS));
01367     
01368   // Array of coefficients matrices from the projection.  
01369   // For our first test, we allocate each of these matrices
01370   // with the proper dimensions.
01371   Array< RCP< mat_type > > C (num_X);
01372   for (int k = 0; k < num_X; ++k)
01373     {
01374       C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
01375       C[k]->random(); // will be overwritten
01376     }
01377   try {
01378     // Q*B := (I - X X^*) S
01379     const int reportedRank = OM->projectAndNormalize (*Q, C, B, X);
01380 
01381     // Pick out the first reportedRank columns of Q.
01382     std::vector<int> indices (reportedRank);
01383     for (int j = 0; j < reportedRank; ++j)
01384       indices[j] = j;
01385     RCP< const MV > Q_left = MVT::CloneView (*Q, indices);
01386     
01387     // Test whether the first reportedRank columns of Q are
01388     // orthogonal.
01389     {
01390       const magnitude_type orthoError = OM->orthonormError (*Q_left);
01391       sout << "-- ||Q(1:" << reportedRank << ")^* Q(1:" << reportedRank 
01392      << ") - I||_F = " << orthoError << endl;
01393       if (orthoError > TOL)
01394         {
01395     sout << "   *** Error: ||Q(1:" << reportedRank << ")^* Q(1:" 
01396          << reportedRank << ") - I||_F = " << orthoError 
01397          << " > TOL = " << TOL << "." << endl;
01398     numerr++;
01399         }
01400     }
01401 
01402     // Compute the residual: if successful, S = Q*B + 
01403     // X (X^* S =: C) in exact arithmetic.  So, the residual is
01404     // S - Q*B - X1 C1 - X2 C2.
01405     //
01406     // Residual := S
01407     MVT::MvAddMv (SCT::one(), *S, SCT::zero(), *Residual, *Residual);
01408     {
01409       // Pick out the first reportedRank rows of B.  Make a deep
01410       // copy, since mat_type is not safe with respect
01411       // to RCP-based memory management (it uses raw pointers
01412       // inside).
01413       RCP< const mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, B->numCols()));
01414       // Residual := Residual - Q(:, 1:reportedRank) * B(1:reportedRank, :)
01415       MVT::MvTimesMatAddMv (-SCT::one(), *Q_left, *B_top, SCT::one(), *Residual);
01416     }
01417     // Residual := Residual - X[k]*C[k]
01418     for (int k = 0; k < num_X; ++k)
01419       MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
01420     const magnitude_type residErr = frobeniusNorm (*Residual);
01421     sout << "-- ||S - Q(:, 1:" << reportedRank << ")*B(1:" 
01422          << reportedRank << ", :) - X1*C1 - X2*C2||_F = " 
01423          << residErr << endl;
01424     if (residErr > ATOL * TOL)
01425       {
01426         sout << "   *** Error: ||S - Q(:, 1:" << reportedRank 
01427        << ")*B(1:" << reportedRank << ", :) "
01428        << "- X1*C1 - X2*C2||_F = " << residErr 
01429        << " > ATOL*TOL = " << (ATOL*TOL) << "." << endl;
01430         numerr++;
01431       }
01432     // Verify that Q(1:reportedRank) is orthogonal to X[k], for
01433     // all k.  This test only makes sense if reportedRank > 0.
01434     if (reportedRank == 0)
01435       {
01436         sout << "-- Reported rank of Q is zero: skipping Q, X[k] "
01437     "orthogonality test." << endl;
01438       }
01439     else
01440       {
01441         for (int k = 0; k < num_X; ++k)
01442     {
01443       // Q should be orthogonal to X[k], for all k.
01444       const magnitude_type projErr = OM->orthogError(*X[k], *Q_left);
01445       sout << "-- ||<Q(1:" << reportedRank << "), X[" << k 
01446            << "]>||_F = " << projErr << endl;
01447       if (projErr > ATOL*TOL)
01448         {
01449           sout << "   *** Error: ||<Q(1:" << reportedRank << "), X["
01450          << k << "]>||_F = " << projErr << " > ATOL*TOL = " 
01451          << (ATOL*TOL) << "." << endl;
01452           numerr++;
01453         }
01454     }
01455       }
01456   } catch (Belos::OrthoError& e) {
01457     sout << "  *** Error: The OrthoManager subclass instance threw "
01458       "an exception: " << e.what() << endl;
01459     numerr++;
01460   }
01461 
01462   // Print out the collected diagnostic messages, which possibly
01463   // include error messages.
01464   const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
01465   MyOM->stream(type) << sout.str();
01466   MyOM->stream(type) << endl;
01467 
01468   return numerr;
01469       }
01470 
01471 
01475       static int
01476       testProjectNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, 
01477           const Teuchos::RCP< const MV >& S, 
01478           const Teuchos::RCP< const MV >& X1, 
01479           const Teuchos::RCP< const MV >& X2,
01480           const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
01481       {
01482   using Teuchos::Array;
01483   using Teuchos::null;
01484   using Teuchos::RCP;
01485   using Teuchos::rcp;
01486   using Teuchos::tuple;
01487 
01488   // We collect all the output in this string wrapper, and print
01489   // it at the end.
01490   std::ostringstream sout;
01491   // Total number of failed tests in this call of this routine.
01492   int numerr = 0;
01493 
01494   const int numRows = MVT::GetVecLength(*S);
01495   const int numCols = MVT::GetNumberVecs(*S);
01496   const int sizeS = MVT::GetNumberVecs(*S);
01497 
01498   // Relative tolerance against which all tests are performed.
01499   // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
01500   // The following bounds hold for all $m \times n$ matrices $A$:
01501   // \[
01502   // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
01503   // \]
01504   // where $r$ is the (column) rank of $A$.  We bound this above
01505   // by the number of columns in $A$.  
01506   //
01507   // Since we are measuring both of these quantitites in the
01508   // Frobenius norm instead, we scale all error tests by 
01509   // $\sqrt{n}$.
01510   //
01511   // A good heuristic is to scale the bound by the square root
01512   // of the number of floating-point operations.  One could
01513   // perhaps support this theoretically, since we are using
01514   // uniform random test problems.
01515   const magnitude_type fudgeFactor = 
01516     SMT::squareroot(magnitude_type(numRows) * 
01517         magnitude_type(numCols) * 
01518         magnitude_type(numCols));
01519   const magnitude_type TOL = SMT::eps() * fudgeFactor *
01520     SMT::squareroot(magnitude_type(numCols));
01521 
01522   // Absolute tolerance scaling: the Frobenius norm of the test
01523   // matrix S.  TOL*ATOL is the absolute tolerance for the
01524   // residual $\|A - Q*B\|_F$.
01525   const magnitude_type ATOL = frobeniusNorm (*S);
01526 
01527   sout << "The test matrix S has Frobenius norm " << ATOL 
01528        << ", and the relative error tolerance is TOL = " 
01529        << TOL << "." << endl;
01530 
01531   // Make some copies of S, X1, and X2.  The OrthoManager's
01532   // project() method shouldn't modify X1 or X2, but this is a a
01533   // test and we don't know that it doesn't!
01534   RCP< MV > S_copy = MVT::CloneCopy(*S);
01535   RCP< MV > Residual = MVT::CloneCopy(*S);
01536   const int num_X = 2;
01537   Array< RCP< const MV > > X (num_X);
01538   X[0] = MVT::CloneCopy(*X1);
01539   X[1] = MVT::CloneCopy(*X2);
01540     
01541   // Array of coefficients matrices from the projection.  
01542   // For our first test, we allocate each of these matrices
01543   // with the proper dimensions.
01544   Array< RCP< mat_type > > C (num_X);
01545   for (int k = 0; k < num_X; ++k)
01546     {
01547       C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
01548       C[k]->random(); // will be overwritten
01549     }
01550   try {
01551     // Compute the projection: S_copy := (I - X X^*) S
01552     OM->project(*S_copy, C, X);
01553 
01554     // Compute the residual: if successful, S = S_copy + X (X^*
01555     // S =: C) in exact arithmetic.  So, the residual is
01556     // S - S_copy - X1 C1 - X2 C2.
01557     //
01558     // Residual := S - S_copy
01559     MVT::MvAddMv (SCT::one(), *S, -SCT::one(), *S_copy, *Residual);
01560     // Residual := Residual - X[k]*C[k]
01561     for (int k = 0; k < num_X; ++k)
01562       MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
01563     magnitude_type residErr = frobeniusNorm (*Residual);
01564     sout << "  ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr;
01565     if (residErr > ATOL * TOL)
01566       {
01567         sout << "  *** Error: ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr 
01568        << " > ATOL*TOL = " << (ATOL*TOL) << ".";
01569         numerr++;
01570       }
01571     for (int k = 0; k < num_X; ++k)
01572       {
01573         // S_copy should be orthogonal to X[k] now.
01574         const magnitude_type projErr = OM->orthogError(*X[k], *S_copy);
01575         if (projErr > TOL)
01576     {
01577       sout << "  *** Error: S is not orthogonal to X[" << k 
01578            << "] by a factor of " << projErr << " > TOL = " 
01579            << TOL << ".";
01580       numerr++;
01581     }
01582       }
01583   } catch (Belos::OrthoError& e) {
01584     sout << "  *** Error: The OrthoManager subclass instance threw "
01585       "an exception: " << e.what() << endl;
01586     numerr++;
01587   }
01588 
01589   // Print out the collected diagnostic messages, which possibly
01590   // include error messages.
01591   const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
01592   MyOM->stream(type) << sout.str();
01593   MyOM->stream(type) << endl;
01594 
01595   return numerr;
01596       }
01597 
01598       static int 
01599       testProject (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, 
01600        const Teuchos::RCP< const MV >& S, 
01601        const Teuchos::RCP< const MV >& X1, 
01602        const Teuchos::RCP< const MV >& X2,
01603        const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
01604       {
01605   return testProjectNew (OM, S, X1, X2, MyOM);
01606       }
01607 
01611       static int 
01612       testProjectOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM, 
01613           const Teuchos::RCP< const MV >& S, 
01614           const Teuchos::RCP< const MV >& X1, 
01615           const Teuchos::RCP< const MV >& X2,
01616           const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
01617       {
01618   using Teuchos::Array;
01619   using Teuchos::null;
01620   using Teuchos::RCP;
01621   using Teuchos::rcp;
01622   using Teuchos::tuple;
01623 
01624   const scalar_type ONE = SCT::one();
01625   // We collect all the output in this string wrapper, and print
01626   // it at the end.
01627   std::ostringstream sout;
01628   // Total number of failed tests in this call of this routine.
01629   int numerr = 0;
01630 
01631   const int numRows = MVT::GetVecLength(*S);
01632   const int numCols = MVT::GetNumberVecs(*S);
01633   const int sizeS = MVT::GetNumberVecs(*S);
01634   const int sizeX1 = MVT::GetNumberVecs(*X1);
01635   const int sizeX2 = MVT::GetNumberVecs(*X2);
01636 
01637   // Relative tolerance against which all tests are performed.
01638   // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
01639   // The following bounds hold for all $m \times n$ matrices $A$:
01640   // \[
01641   // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
01642   // \]
01643   // where $r$ is the (column) rank of $A$.  We bound this above
01644   // by the number of columns in $A$.  
01645   //
01646   // Since we are measuring both of these quantitites in the
01647   // Frobenius norm instead, we scale all error tests by 
01648   // $\sqrt{n}$.
01649   //
01650   // A good heuristic is to scale the bound by the square root
01651   // of the number of floating-point operations.  One could
01652   // perhaps support this theoretically, since we are using
01653   // uniform random test problems.
01654   const magnitude_type fudgeFactor = 
01655     SMT::squareroot(magnitude_type(numRows) * 
01656         magnitude_type(numCols) * 
01657         magnitude_type(numCols));
01658   const magnitude_type TOL = SMT::eps() * fudgeFactor *
01659     SMT::squareroot(magnitude_type(numCols));
01660 
01661   // Absolute tolerance scaling: the Frobenius norm of the test
01662   // matrix S.  TOL*ATOL is the absolute tolerance for the
01663   // residual $\|A - Q*B\|_F$.
01664   const magnitude_type ATOL = frobeniusNorm (*S);
01665 
01666   sout << "The test matrix S has Frobenius norm " << ATOL 
01667        << ", and the relative error tolerance is TOL = " 
01668        << TOL << "." << endl;
01669 
01670 
01671   //
01672   // Output tests:
01673   //   <S_out,X1> = 0
01674   //   <S_out,X2> = 0
01675   //   S_in = S_out + X1 C1 + X2 C2
01676   // 
01677   // We will loop over an integer specifying the test combinations.
01678   // The bit pattern for the different tests is listed in parentheses.
01679   //
01680   // For the projectors, test the following combinations:
01681   // none              (00)
01682   // P_X1              (01)
01683   // P_X2              (10)
01684   // P_X1 P_X2         (11)
01685   // P_X2 P_X1         (11)
01686   // The latter two should be tested to give the same result.
01687   //
01688   // For each of these, we should test with C1 and C2:
01689   //
01690   // if hasM:
01691   // with and without MX1   (1--) 
01692   // with and without MX2  (1---) 
01693   // with and without MS  (1----) 
01694   //
01695   // As hasM controls the upper level bits, we need only run test
01696   // cases 0-3 if hasM==false.  Otherwise, we run test cases 0-31.
01697   //
01698 
01699   int numtests = 8;
01700 
01701   // test ortho error before orthonormalizing
01702   if (X1 != null) {
01703     magnitude_type err = OM->orthogError(*S,*X1);
01704     sout << "   || <S,X1> || before     : " << err << endl;
01705   }
01706   if (X2 != null) {
01707     magnitude_type err = OM->orthogError(*S,*X2);
01708     sout << "   || <S,X2> || before     : " << err << endl;
01709   }
01710 
01711   for (int t = 0; t < numtests; ++t) 
01712     {
01713       Array< RCP< const MV > > theX;
01714       Array< RCP< mat_type > > C;
01715       if ( (t && 3) == 0 ) {
01716         // neither X1 nor X2
01717         // C and theX are already empty
01718       }
01719       else if ( (t && 3) == 1 ) {
01720         // X1
01721         theX = tuple(X1);
01722         C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
01723       }
01724       else if ( (t && 3) == 2 ) {
01725         // X2
01726         theX = tuple(X2);
01727         C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
01728       }
01729       else {
01730         // X1 and X2, and the reverse.
01731         theX = tuple(X1,X2);
01732         C = tuple( rcp(new mat_type(sizeX1,sizeS)), 
01733        rcp(new mat_type(sizeX2,sizeS)) );
01734       }
01735 
01736       try {
01737         // call routine
01738         // if (t && 3) == 3, {
01739         //    call with reversed input: X2 X1
01740         // }
01741         // test all outputs for correctness
01742         // test all outputs for equivalence
01743 
01744         // here is where the outputs go
01745         Array< RCP< MV > > S_outs;
01746         Array< Array< RCP< mat_type > > > C_outs;
01747         RCP< MV > Scopy;
01748 
01749         // copies of S,MS
01750         Scopy = MVT::CloneCopy(*S);
01751         // randomize this data, it should be overwritten
01752         for (size_type i = 0; i < C.size(); ++i) {
01753     C[i]->random();
01754         }
01755         // Run test.
01756         // Note that Anasazi and Belos differ, among other places, 
01757         // in the order of arguments to project().
01758         OM->project(*Scopy,C,theX);
01759         // we allocate S and MS for each test, so we can save these as views
01760         // however, save copies of the C
01761         S_outs.push_back( Scopy );
01762         C_outs.push_back( Array< RCP< mat_type > >(0) );
01763         if (C.size() > 0) {
01764     C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
01765         }
01766         if (C.size() > 1) {
01767     C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
01768         }
01769 
01770         // do we run the reversed input?
01771         if ( (t && 3) == 3 ) {
01772     // copies of S,MS
01773     Scopy = MVT::CloneCopy(*S);
01774     // randomize this data, it should be overwritten
01775     for (size_type i = 0; i < C.size(); ++i) {
01776       C[i]->random();
01777     }
01778     // flip the inputs
01779     theX = tuple( theX[1], theX[0] );
01780     // Run test.
01781     // Note that Anasazi and Belos differ, among other places, 
01782     // in the order of arguments to project().
01783     OM->project(*Scopy,C,theX);
01784     // we allocate S and MS for each test, so we can save these as views
01785     // however, save copies of the C
01786     S_outs.push_back( Scopy );
01787     // we are in a special case: P_X1 and P_X2, so we know we applied 
01788     // two projectors, and therefore have two C[i]
01789     C_outs.push_back( Array<RCP<mat_type > >() );
01790     // reverse the Cs to compensate for the reverse projectors
01791     C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
01792     C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
01793     // flip the inputs back
01794     theX = tuple( theX[1], theX[0] );
01795         }
01796 
01797         // test all outputs for correctness
01798         for (size_type o = 0; o < S_outs.size(); ++o) {
01799     // S_in = X1*C1 + C2*C2 + S_out
01800     {
01801       RCP<MV> tmp = MVT::CloneCopy(*S_outs[o]);
01802       if (C_outs[o].size() > 0) {
01803         MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
01804         if (C_outs[o].size() > 1) {
01805           MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
01806         }
01807       }
01808       magnitude_type err = MVDiff(*tmp,*S);
01809       if (err > ATOL*TOL) {
01810         sout << "         vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv         tolerance exceeded! test failed!" << endl;
01811         numerr++;
01812       }
01813       sout << "  " << t << "|| S_in - X1*C1 - X2*C2 - S_out || : " << err << endl;
01814     }
01815     // <X1,S> == 0
01816     if (theX.size() > 0 && theX[0] != null) {
01817       magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
01818       if (err > TOL) {
01819         sout << "         vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv         tolerance exceeded! test failed!" << endl;
01820         numerr++;
01821       }
01822       sout << "  " << t << "|| <X[0],S> || after      : " << err << endl;
01823     }
01824     // <X2,S> == 0
01825     if (theX.size() > 1 && theX[1] != null) {
01826       magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
01827       if (err > TOL) {
01828         sout << "         vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv         tolerance exceeded! test failed!" << endl;
01829         numerr++;
01830       }
01831       sout << "  " << t << "|| <X[1],S> || after      : " << err << endl;
01832     }
01833         }
01834 
01835         // test all outputs for equivalence
01836         // check all combinations:
01837         //    output 0 == output 1
01838         //    output 0 == output 2
01839         //    output 1 == output 2
01840         for (size_type o1=0; o1<S_outs.size(); o1++) {
01841     for (size_type o2=o1+1; o2<S_outs.size(); o2++) {
01842       // don't need to check MS_outs because we check 
01843       //   S_outs and MS_outs = M*S_outs
01844       // don't need to check C_outs either
01845       //   
01846       // check that S_outs[o1] == S_outs[o2]
01847       magnitude_type err = MVDiff(*S_outs[o1],*S_outs[o2]);
01848       if (err > TOL) {
01849         sout << "    vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv         tolerance exceeded! test failed!" << endl;
01850         numerr++;
01851       }
01852     }
01853         }
01854 
01855       }
01856       catch (Belos::OrthoError& e) {
01857         sout << "   -------------------------------------------         project() threw exception" << endl;
01858         sout << "   Error: " << e.what() << endl;
01859         numerr++;
01860       }
01861     } // test for
01862 
01863   MsgType type = Debug;
01864   if (numerr>0) type = Errors;
01865   MyOM->stream(type) << sout.str();
01866   MyOM->stream(type) << endl;
01867 
01868   return numerr;
01869       }
01870 
01871 
01872     };
01873 
01874 
01875 
01876   } // namespace Test
01877 } // namespace Belos
01878 
01879 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines