Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_FullTsqrTest.hpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2009) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ************************************************************************
00027 //@HEADER
00028 
00029 #ifndef __TSQR_Test_FullTsqrTest_hpp
00030 #define __TSQR_Test_FullTsqrTest_hpp
00031 
00032 #include <Tsqr.hpp>
00033 #include <Tsqr_Random_NormalGenerator.hpp>
00034 #include <Tsqr_Random_GlobalMatrix.hpp>
00035 #include <Tsqr_TestSetup.hpp>
00036 //#include <TsqrFactory_SequentialTsqr.hpp>
00037 #include <Tsqr_GlobalVerify.hpp>
00038 #include <Tsqr_TeuchosMessenger.hpp>
00039 #include "Tsqr_TestUtils.hpp"
00040 #include <Teuchos_ScalarTraits.hpp>
00041 
00042 #include <iostream>
00043 #include <stdexcept>
00044 #include <string>
00045 
00046 namespace TSQR {
00047   namespace Test {
00048 
00052     class TsqrInaccurate : public std::exception {
00053     public:
00055       TsqrInaccurate (const std::string& msg) : msg_ (msg) {}
00056 
00058       const char* what() const throw() { return msg_.c_str(); }
00059 
00061       virtual ~TsqrInaccurate() throw() {}
00062 
00063     private:
00064       std::string msg_;
00065     };
00066 
00091     template<class Scalar>
00092     class FullTsqrVerifier {
00093     public:
00094       typedef Scalar scalar_type;
00095       typedef int ordinal_type;
00096       typedef SequentialTsqr<ordinal_type, scalar_type> node_tsqr_type;
00097       typedef DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
00098       typedef Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
00099       typedef Kokkos::SerialNode node_type;
00100 
00101     private:
00102 
00104       static Teuchos::RCP<tsqr_type>
00105       getTsqr (const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00106          const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
00107       {
00108   using Teuchos::ParameterList;
00109   using Teuchos::parameterList;
00110   using Teuchos::rcp_implicit_cast;
00111   using Teuchos::RCP;
00112   using Teuchos::rcp;
00113 
00114   const size_t cacheSizeHint = testParams->get<size_t> ("cacheSizeHint");
00115   //const int numCores = testParams->get<int> ("numCores");
00116 
00117   //RCP<ParameterList> tsqrParams = parameterList ("Intranode TSQR");
00118   //tsqrParams->set ("cacheSizeHint", cacheSizeHint);
00119   //tsqrParams->set ("numCores", numCores);
00120 
00121   RCP<node_tsqr_type> seqTsqr = rcp (new node_tsqr_type (cacheSizeHint));
00122 
00123   RCP<TeuchosMessenger<scalar_type> > scalarMess =
00124     rcp (new TeuchosMessenger<scalar_type> (comm));
00125   RCP<MessengerBase<scalar_type> > scalarMessBase = 
00126     rcp_implicit_cast<MessengerBase<scalar_type> > (scalarMess);
00127   RCP<dist_tsqr_type> distTsqr = 
00128     rcp (new dist_tsqr_type (scalarMessBase));
00129 
00130   return rcp (new tsqr_type (seqTsqr, distTsqr));
00131       }
00132 
00133     public:
00134 
00144       static void
00145       run (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00146      const Teuchos::RCP<const node_type>& node,
00147      const Teuchos::RCP<Teuchos::ParameterList>& testParams,
00148      std::vector<int>& randomSeed)
00149       {
00150   using std::cerr;
00151   using std::cout;
00152   using std::endl;
00153   using Teuchos::arcp;
00154   using Teuchos::ParameterList;
00155   using Teuchos::parameterList;
00156   using Teuchos::RCP;
00157   using Teuchos::rcp;
00158   using Teuchos::rcp_const_cast;
00159   using Teuchos::rcp_implicit_cast;
00160   typedef Matrix<ordinal_type, scalar_type> matrix_type;
00161   typedef MatView<ordinal_type, scalar_type> view_type;
00162   typedef typename tsqr_type::FactorOutput factor_output_type;
00163 
00164   const int myRank = Teuchos::rank (*comm);
00165   const int numProcs = Teuchos::size (*comm);
00166 
00167   // Construct TSQR implementation instance.
00168   RCP<tsqr_type> tsqr = getTsqr (testParams, comm);
00169 
00170   // Fetch test parameters from the input parameter list.
00171   const ordinal_type numRowsLocal = testParams->get<ordinal_type> ("numRowsLocal");
00172   const ordinal_type numCols = testParams->get<ordinal_type> ("numCols"); 
00173   const int numCores = testParams->get<int> ("numCores");
00174   const bool contiguousCacheBlocks = testParams->get<bool> ("contiguousCacheBlocks");
00175   const bool testFactorExplicit = testParams->get<bool> ("testFactorExplicit");
00176   const bool testRankRevealing = testParams->get<bool> ("testRankRevealing");
00177   const bool debug = testParams->get<bool> ("debug");
00178 
00179   // Space for each node's local part of the test problem.
00180   // A_local, A_copy, and Q_local are distributed matrices, and
00181   // R is replicated on all processes sharing the communicator.
00182   matrix_type A_local (numRowsLocal, numCols);
00183   matrix_type A_copy (numRowsLocal, numCols);
00184   matrix_type Q_local (numRowsLocal, numCols);
00185   matrix_type R (numCols, numCols);
00186 
00187   // Start out by filling the test problem with zeros.
00188   typedef Teuchos::ScalarTraits<scalar_type> STS;
00189   A_local.fill (STS::zero());
00190   A_copy.fill (STS::zero());
00191   Q_local.fill (STS::zero());
00192   R.fill (STS::zero());
00193 
00194   // Create some reasonable singular values for the test problem:
00195   // 1, 1/2, 1/4, 1/8, ...
00196   typedef typename STS::magnitudeType magnitude_type;
00197   std::vector<magnitude_type> singularValues (numCols);
00198   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00199   {
00200     const magnitude_type scalingFactor = STM::one() + STM::one();
00201     magnitude_type curVal = STM::one();
00202     typedef typename std::vector<magnitude_type>::iterator iter_type;
00203     for (iter_type it = singularValues.begin(); 
00204          it != singularValues.end(); ++it)
00205       {
00206         *it = curVal;
00207         curVal = curVal / scalingFactor;
00208       }
00209   }
00210 
00211   // Construct a normal(0,1) pseudorandom number generator with
00212   // the given random seed.
00213   using TSQR::Random::NormalGenerator;
00214   typedef NormalGenerator<ordinal_type, scalar_type> generator_type;
00215   generator_type gen (randomSeed);
00216 
00217   // We need a Messenger for Ordinal-type data, so that we can
00218   // build a global random test matrix.
00219   RCP<MessengerBase<ordinal_type> > ordinalMessenger = 
00220     rcp_implicit_cast<MessengerBase<ordinal_type> > (rcp (new TeuchosMessenger<ordinal_type> (comm)));
00221 
00222   // We also need a Messenger for Scalar-type data.  The TSQR
00223   // implementation already constructed one, but it's OK to
00224   // construct another one; TeuchosMessenger is just a thin
00225   // wrapper over the Teuchos::Comm object.
00226   RCP<MessengerBase<scalar_type> > scalarMessenger =
00227     rcp_implicit_cast<MessengerBase<scalar_type> > (rcp (new TeuchosMessenger<scalar_type> (comm)));
00228 
00229   {
00230     // Generate a global distributed matrix (whose part local to
00231     // this node is in A_local) with the given singular values.
00232     // This part has O(P) communication for P MPI processes.
00233     using TSQR::Random::randomGlobalMatrix;
00234     // Help the C++ compiler with type inference.
00235     view_type A_local_view (A_local.nrows(), A_local.ncols(), A_local.get(), A_local.lda());
00236     const magnitude_type* const singVals = (numCols == 0) ? NULL : &singularValues[0];
00237     randomGlobalMatrix<view_type, generator_type> (&gen, A_local_view, singVals,
00238                ordinalMessenger.getRawPtr(),
00239                scalarMessenger.getRawPtr());
00240   }
00241   // Save the pseudorandom number generator's seed for any later
00242   // tests.  The generator keeps its own copy of the seed and
00243   // updates it internally, so we have to ask for its copy.
00244   gen.getSeed (randomSeed);
00245 
00246   // If specified in the test parameters, rearrange cache blocks
00247   // in the copy.  Otherwise, just copy the test problem into
00248   // A_copy.  The factorization overwrites the input matrix, so
00249   // we have to make a copy in order to validate the final
00250   // result.
00251   if (contiguousCacheBlocks)
00252     {
00253       tsqr->cache_block (numRowsLocal, numCols, A_copy.get(), 
00254              A_local.get(), A_local.lda());
00255       if (debug)
00256         {
00257     Teuchos::barrier (*comm);
00258     if (myRank == 0)
00259       cerr << "-- Finished Tsqr::cache_block" << endl;
00260         }
00261     }
00262   else
00263     A_copy.copy (A_local);
00264 
00265   // "factorExplicit" is an alternate, hopefully faster way of
00266   // factoring the matrix, when only the explicit Q factor is
00267   // wanted.
00268   typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00269   if (testFactorExplicit)
00270     {
00271       // Kokkos::MultiVector wants a non-const Node, for some reason.
00272       KMV A_copy_view (rcp_const_cast<node_type> (node));
00273       A_copy_view.initializeValues (static_cast<size_t> (A_copy.nrows()), 
00274             static_cast<size_t> (A_copy.ncols()), 
00275             arcp (A_copy.get(), 0, A_copy.nrows()*A_copy.ncols(), false), // non-owning ArrayRCP
00276             static_cast<size_t> (A_copy.lda()));
00277       // Kokkos::MultiVector wants a non-const Node, for some reason.
00278       KMV Q_view (rcp_const_cast<node_type> (node));
00279       Q_view.initializeValues (static_cast<size_t> (Q_local.nrows()), 
00280              static_cast<size_t> (Q_local.ncols()),
00281              arcp (Q_local.get(), 0, Q_local.nrows()*Q_local.ncols(), false), // non-owning ArrayRCP
00282              static_cast<size_t> (Q_local.lda()));
00283       Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> 
00284         R_view (Teuchos::View, R.get(), R.lda(), R.nrows(), R.ncols());
00285 
00286       tsqr->factorExplicit (A_copy_view, Q_view, R_view,
00287           contiguousCacheBlocks);
00288       if (debug)
00289         {
00290     Teuchos::barrier (*comm);
00291     if (myRank == 0)
00292       cerr << "-- Finished Tsqr::factorExplicit" << endl;
00293         }
00294     }
00295   else
00296     {
00297       // Factor the (copy of the) matrix.
00298       factor_output_type factorOutput = 
00299         tsqr->factor (numRowsLocal, numCols, A_copy.get(), A_copy.lda(), 
00300           R.get(), R.lda(), contiguousCacheBlocks);
00301       if (debug)
00302         {
00303     Teuchos::barrier (*comm);
00304     if (myRank == 0)
00305       cerr << "-- Finished Tsqr::factor" << endl;
00306         }
00307       // Compute the explicit Q factor in Q_local.
00308       tsqr->explicit_Q (numRowsLocal, numCols, A_copy.get(), A_copy.lda(), 
00309             factorOutput, numCols, Q_local.get(), Q_local.lda(),
00310             contiguousCacheBlocks);
00311       if (debug)
00312         {
00313     Teuchos::barrier (*comm);
00314     if (myRank == 0)
00315       cerr << "-- Finished Tsqr::explicit_Q" << endl;
00316         }
00317     }
00318 
00319   // Optionally, test rank-revealing capability.  We do this
00320   // before un-cache-blocking the explicit Q factor, since
00321   // revealRank can work with contiguous cache blocks, and
00322   // modifies the Q factor if the matrix doesn't have full
00323   // column rank.
00324   if (testRankRevealing)
00325     {
00326       // Kokkos::MultiVector wants a non-const Node, for some reason.
00327       KMV Q_view (rcp_const_cast<node_type> (node));
00328       Q_view.initializeValues (static_cast<size_t> (Q_local.nrows()), 
00329              static_cast<size_t> (Q_local.ncols()),
00330              arcp (Q_local.get(), 0, Q_local.nrows()*Q_local.ncols(), false), // non-owning ArrayRCP
00331              static_cast<size_t> (Q_local.lda()));
00332       Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> 
00333         R_view (Teuchos::View, R.get(), R.lda(), R.nrows(), R.ncols());
00334       // If 2^{# columns} > machine precision, then our choice
00335       // of singular values will make the smallest singular
00336       // value < machine precision.  In that case, the SVD can't
00337       // promise it will distinguish between tiny and zero.  If
00338       // the number of columns is less than that, we can use a
00339       // tolerance of zero to test the purported rank with the
00340       // actual numerical rank.
00341       const magnitude_type tol = STM::zero();
00342       const ordinal_type rank = 
00343         tsqr->revealRank (Q_view, R_view, tol, contiguousCacheBlocks);
00344 
00345       magnitude_type two_to_the_numCols = STM::one();
00346       for (int k = 0; k < numCols; ++k)
00347         {
00348     const magnitude_type two = STM::one() + STM::one();
00349     two_to_the_numCols *= two;
00350         }
00351       // Throw in a factor of 10, just for more tolerance of
00352       // rounding error (so the test only fails if something is
00353       // really broken).
00354       if (two_to_the_numCols > magnitude_type(10) * STM::eps())
00355         {
00356     TEST_FOR_EXCEPTION(rank != numCols, std::logic_error,
00357            "The matrix of " << numCols << " columns "
00358            "should have full numerical rank, but Tsqr "
00359            "reports that it has rank " << rank << ".  "
00360            "Please report this bug to the Kokkos "
00361            "developers.");
00362     if (debug)
00363       {
00364         Teuchos::barrier (*comm);
00365         if (myRank == 0)
00366           cerr << "-- Tested rank-revealing capability" << endl;
00367       }
00368         }
00369       else
00370         {
00371     if (debug)
00372       {
00373         Teuchos::barrier (*comm);
00374         if (myRank == 0)
00375           cerr << "-- Not testing rank-revealing capability; too many columns" << endl;
00376       }
00377         }
00378     }
00379   // "Un"-cache-block the output, if contiguous cache blocks
00380   // were used.  This is only necessary because global_verify()
00381   // doesn't currently support contiguous cache blocks.
00382   if (contiguousCacheBlocks)
00383     {
00384       // We can use A_copy as scratch space for
00385       // un-cache-blocking Q_local, since we're done using
00386       // A_copy for other things.
00387       tsqr->un_cache_block (numRowsLocal, numCols, A_copy.get(), 
00388           A_copy.lda(), Q_local.get());
00389       // Overwrite Q_local with the un-cache-blocked Q factor.
00390       Q_local.copy (A_copy);
00391       if (debug)
00392         {
00393     Teuchos::barrier (*comm);
00394     if (myRank == 0)
00395       cerr << "-- Finished Tsqr::un_cache_block" << endl;
00396         }
00397     }
00398 
00399   // Test accuracy of the factorization.
00400   const std::vector<magnitude_type> results = 
00401     global_verify (numRowsLocal, numCols, A_local.get(), A_local.lda(),
00402        Q_local.get(), Q_local.lda(), R.get(), R.lda(), 
00403        scalarMessenger.getRawPtr());
00404   if (debug)
00405     {
00406       Teuchos::barrier (*comm);
00407       if (myRank == 0)
00408         cerr << "-- Finished global_verify" << endl;
00409     }
00410 
00411   // Print the results on Proc 0.
00412   if (myRank == 0)
00413     {
00414       if (testParams->get<bool> ("printFieldNames"))
00415         {
00416     cout << "%"
00417          << "method"
00418          << ",scalarType"
00419          << ",numRowsLocal"
00420          << ",numCols"
00421          << ",numProcs"
00422          << ",numCores"
00423          << ",cacheSizeHint"
00424          << ",contiguousCacheBlocks"
00425          << ",absFrobResid"
00426          << ",absFrobOrthog"
00427          << ",frobA" << endl;
00428     // We don't need to print field names again for the other tests, so set the test parameters accordingly.
00429     testParams->set ("printFieldNames", false);
00430         }
00431       if (testParams->get<bool> ("printResults"))
00432         {
00433     cout << "Tsqr"
00434          << "," << Teuchos::TypeNameTraits<scalar_type>::name()
00435          << "," << numRowsLocal
00436          << "," << numCols
00437          << "," << numProcs
00438          << "," << numCores
00439          << "," << tsqr->cache_size_hint()
00440          << "," << contiguousCacheBlocks 
00441          << "," << results[0] 
00442          << "," << results[1]
00443          << "," << results[2]
00444          << endl;
00445         }
00446     } // if (myRank == 0)
00447 
00448   // If requested, check accuracy and fail if results are not
00449   // sufficiently accurate.
00450   if (testParams->get<bool> ("failIfInaccurate"))
00451     {
00452       // Avoid overflow of the local Ordinal type, by casting
00453       // first to a floating-point type.
00454       const magnitude_type dimsProd = magnitude_type(numRowsLocal) * 
00455         magnitude_type(numProcs) * magnitude_type(numCols*numCols);
00456 
00457       // Relative residual error is ||A-Q*R|| / ||A||, or just
00458       // ||A-Q*R|| if ||A|| == 0.  (The result had better be
00459       // zero in the latter case.)  A reasonable error bound
00460       // should incorporate the dimensions of the matrix, since
00461       // this indicates the amount of rounding error.  Square
00462       // root of the matrix dimensions is an old heuristic from
00463       // Wilkinson or perhaps even an earlier source.  We
00464       // include a factor of 10 so that the test won't fail
00465       // unless there is a really good reason.
00466       const magnitude_type relResidBound = 
00467         magnitude_type(10) * STM::squareroot(dimsProd) * STM::eps();
00468 
00469       // Orthogonality of the matrix should not depend on the
00470       // matrix dimensions, if we measure in the 2-norm.
00471       // However, we are measuring in the Frobenius norm, so
00472       // it's appropriate to multiply eps by the number of
00473       // entries in the matrix for which we compute the
00474       // Frobenius norm.  We include a factor of 10 for the same
00475       // reason as mentioned above.
00476       const magnitude_type orthoBound = 
00477         magnitude_type(10*numCols*numCols) * STM::eps();
00478 
00479       // Avoid division by zero.
00480       const magnitude_type relResidError = 
00481         results[0] / (results[2] == STM::zero() ? STM::one() : results[2]);
00482       TEST_FOR_EXCEPTION(relResidError > relResidBound, TsqrInaccurate,
00483              "Full Tsqr (SequentialTsqr + DistTsqr) has an "
00484              "inaccurate relative residual ||A - QR||_F"
00485              << (results[2] == STM::zero() ? " / ||A||_F" : "")
00486              << " = " << relResidError << ", which is greater"
00487              " than the bound " << relResidBound << " by a "
00488              "factor of " << relResidError / relResidBound
00489              << ".");
00490       const magnitude_type orthoError = results[1];
00491       TEST_FOR_EXCEPTION(orthoError > orthoBound, TsqrInaccurate,
00492              "Full Tsqr (SequentialTsqr + DistTsqr) has an "
00493              "inaccurate orthogonality measure ||I - Q^* Q||"
00494              "_F" << results[1] << " = " << orthoError 
00495              << ", which is greater than the bound " 
00496              << orthoBound << " by a factor of " 
00497              << orthoError / orthoBound << ".");
00498     } // if (the tests should fail on inaccuracy)
00499       }
00500     };
00501 
00521     template<class TypeListType>
00522     class FullTsqrVerifierCallerImpl {
00523     public:
00524       typedef Kokkos::SerialNode node_type; 
00525 
00526       static void
00527       run (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00528      const Teuchos::RCP<const node_type>& node,
00529      const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00530      std::vector<int>& randomSeed);
00531     };
00532 
00533     //
00534     // Partial specialization for Cons<CarType, CdrType>.
00535     //
00536     template<class CarType, class CdrType>
00537     class FullTsqrVerifierCallerImpl<TSQR::Test::Cons<CarType, CdrType> > {
00538     public:
00539       typedef Kokkos::SerialNode node_type; 
00540 
00541       static void
00542       run (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00543      const Teuchos::RCP<const node_type>& node,
00544      const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00545      std::vector<int>& randomSeed)
00546       {
00547   typedef CarType car_type;
00548   typedef CdrType cdr_type;
00549   FullTsqrVerifier<car_type>::run (comm, node, testParams, randomSeed);
00550   FullTsqrVerifierCallerImpl<cdr_type>::run (comm, node, testParams, randomSeed);
00551       }
00552     };
00553 
00554     //
00555     // Full specialization for NullCons.
00556     //
00557     template<>
00558     class FullTsqrVerifierCallerImpl<TSQR::Test::NullCons> {
00559     public:
00560       typedef Kokkos::SerialNode node_type; 
00561 
00562       static void
00563       run (const Teuchos::RCP<const Teuchos::Comm<int> >&,
00564      const Teuchos::RCP<const node_type>&,
00565      const Teuchos::RCP<Teuchos::ParameterList>&, 
00566      std::vector<int>&)
00567       {
00568   // We're at the end of the type list, so do nothing.
00569       }
00570     };
00571 
00579     class FullTsqrVerifierCaller {
00580     public:
00583       typedef Kokkos::SerialNode node_type;
00584 
00594       typedef int ordinal_type;
00595 
00601       Teuchos::RCP<const Teuchos::ParameterList>
00602       getValidParameterList () const
00603       {
00604   using Teuchos::ParameterList;
00605   using Teuchos::parameterList;
00606   using Teuchos::RCP;
00607 
00608   RCP<ParameterList> plist = parameterList ("FullTsqrVerifier");
00609 
00610   const size_t cacheSizeHint = 0;
00611   const int numCores = 1;
00612   const ordinal_type numRowsLocal = 100;
00613   const ordinal_type numCols = 10;
00614   const bool contiguousCacheBlocks = false;
00615   const bool testFactorExplicit = true;
00616   const bool testRankRevealing = true;
00617   const bool printFieldNames = true;
00618   const bool printResults = true;
00619   const bool failIfInaccurate = true;
00620   const bool debug = false;
00621 
00622   // Parameters for configuring Tsqr itself.
00623   plist->set ("cacheSizeHint", cacheSizeHint, 
00624         "Cache size hint in bytes.  "
00625         "Zero means TSQR picks a reasonable default.");
00626   plist->set ("numCores", numCores,
00627         "Number of partition(s) to use for TbbTsqr (if "
00628         "applicable).  Must be a positive integer.");
00629 
00630   // Parameters for testing Tsqr.
00631   plist->set ("numRowsLocal", numRowsLocal, 
00632         "Number of rows per (MPI) process in the test matrix.  "
00633         "Must be >= the number of columns.");
00634   plist->set ("numCols", numCols, 
00635         "Number of columns in the test matrix.");
00636   plist->set ("contiguousCacheBlocks", contiguousCacheBlocks, 
00637         "Whether to test the factorization with contiguously "
00638         "stored cache blocks.");
00639   plist->set ("testFactorExplicit", testFactorExplicit, 
00640         "Whether to test TSQR's factorExplicit() (a hopefully "
00641         "faster path than calling factor() and explicit_Q() in "
00642         "sequence).");
00643   plist->set ("testRankRevealing", testRankRevealing, 
00644         "Whether to test TSQR's rank-revealing capability.");
00645   plist->set ("printFieldNames", printFieldNames, 
00646         "Whether to print field names (this is only done once, "
00647         "for all Scalar types tested).");
00648   plist->set ("printResults", printResults, 
00649         "Whether to print test results.");
00650   plist->set ("failIfInaccurate", failIfInaccurate,
00651         "Whether to fail the test if the factorization "
00652         "is not sufficiently accurate.");
00653   plist->set ("debug", debug, 
00654         "Whether to print debugging output.");
00655   return plist;
00656       }
00657 
00669       template<class TypeListType>
00670       void 
00671       run (const Teuchos::RCP<Teuchos::ParameterList>& testParams)
00672       {
00673   // Using a class with a static method is a way to implement
00674   // "partial specialization of function templates" (which by
00675   // itself is not allowed in C++).
00676   typedef FullTsqrVerifierCallerImpl<TypeListType> impl_type;
00677   impl_type::run (comm_, node_, testParams, randomSeed_);
00678       }
00679 
00680 
00695       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00696             const Teuchos::RCP<const node_type>& node,
00697             const std::vector<int>& randomSeed) :
00698   comm_ (comm), node_ (node), randomSeed_ (validateRandomSeed (randomSeed))
00699       {}
00700 
00708       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm) :
00709   comm_ (comm),
00710   node_ (getNode<node_type> (getValidNodeParameters<node_type> ())),
00711   randomSeed_ (defaultRandomSeed ())
00712       {}
00713 
00723       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00724             const Teuchos::RCP<const node_type>& node) :
00725   comm_ (comm),
00726   node_ (node),
00727   randomSeed_ (defaultRandomSeed ())
00728       {}
00729 
00731       static std::vector<int> 
00732       validateRandomSeed (const std::vector<int>& seed) 
00733       {
00734   TEST_FOR_EXCEPTION(seed.size() < 4, std::invalid_argument,
00735          "Invalid random seed: Need an array of four integers.");
00736   for (std::vector<int>::size_type k = 0; k < seed.size(); ++k)
00737     {
00738       TEST_FOR_EXCEPTION(seed[k] < 0 || seed[k] > 4095,
00739              std::invalid_argument,
00740              "Invalid random seed: Each of the four integers must be in [0, 4095].");
00741     }
00742   TEST_FOR_EXCEPTION(seed[3] % 2 != 1, std::invalid_argument,
00743          "Invalid random seed: The last of the four integers must be odd.");
00744   return seed;
00745       }
00746 
00748       static std::vector<int> 
00749       defaultRandomSeed () 
00750       {
00751   std::vector<int> seed (4);
00752   seed[0] = 0;
00753   seed[1] = 0;
00754   seed[2] = 0;
00755   seed[3] = 1;
00756   return seed;
00757       }
00758 
00759     private:
00764       Teuchos::RCP<const Teuchos::Comm<int> > comm_;
00765 
00767       Teuchos::RCP<const node_type> node_;
00768 
00774       std::vector<int> randomSeed_;
00775     };
00776 
00777   } // namespace Test
00778 } // namespace TSQR
00779 
00780 #endif // __TSQR_Test_FullTsqrTest_hpp
00781 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends