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