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 (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef __TSQR_Test_FullTsqrTest_hpp
00043 #define __TSQR_Test_FullTsqrTest_hpp
00044 
00045 #include <Tsqr.hpp>
00046 #include <Tsqr_Random_NormalGenerator.hpp>
00047 #include <Tsqr_Random_GlobalMatrix.hpp>
00048 #include <Tsqr_TestSetup.hpp>
00049 //#include <TsqrFactory_SequentialTsqr.hpp>
00050 #include <Tsqr_GlobalVerify.hpp>
00051 #include <Tsqr_TeuchosMessenger.hpp>
00052 #include "Tsqr_TestUtils.hpp"
00053 #include <Teuchos_ScalarTraits.hpp>
00054 
00055 #include <iostream>
00056 #include <stdexcept>
00057 #include <string>
00058 
00059 namespace TSQR {
00060   namespace Test {
00061 
00065     class TsqrInaccurate : public std::exception {
00066     public:
00068       TsqrInaccurate (const std::string& msg) : msg_ (msg) {}
00069 
00071       const char* what() const throw() { return msg_.c_str(); }
00072 
00074       virtual ~TsqrInaccurate() throw() {}
00075 
00076     private:
00077       std::string msg_;
00078     };
00079 
00104     template<class Scalar>
00105     class FullTsqrVerifier {
00106     public:
00107       typedef Scalar scalar_type;
00108       typedef int ordinal_type;
00109       typedef SequentialTsqr<ordinal_type, scalar_type> node_tsqr_type;
00110       typedef DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
00111       typedef Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
00112       typedef Kokkos::SerialNode node_type;
00113 
00114     private:
00115 
00117       static Teuchos::RCP<tsqr_type>
00118       getTsqr (const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00119          const Teuchos::RCP<node_type>& node,
00120          const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
00121       {
00122   using Teuchos::ParameterList;
00123   using Teuchos::parameterList;
00124   using Teuchos::rcp_implicit_cast;
00125   using Teuchos::RCP;
00126   using Teuchos::rcp;
00127 
00128   const size_t cacheSizeHint = testParams->get<size_t> ("cacheSizeHint");
00129   //const int numTasks = testParams->get<int> ("numTasks");
00130 
00131   //RCP<ParameterList> tsqrParams = parameterList ("NodeTsqr");
00132   //tsqrParams->set ("Cache Size Hint", cacheSizeHint);
00133   //tsqrParams->set ("Num Tasks", numCores);
00134 
00135   // TODO (mfh 21 Oct 2011) Some node_tsqr_type classes need a
00136   // Kokkos Node instance.  SequentialTsqr doesn't, so this code
00137   // should be fine for now.
00138   (void) node;
00139   RCP<node_tsqr_type> seqTsqr = rcp (new node_tsqr_type (cacheSizeHint));
00140 
00141   RCP<TeuchosMessenger<scalar_type> > scalarMess =
00142     rcp (new TeuchosMessenger<scalar_type> (comm));
00143   RCP<MessengerBase<scalar_type> > scalarMessBase = 
00144     rcp_implicit_cast<MessengerBase<scalar_type> > (scalarMess);
00145   RCP<dist_tsqr_type> distTsqr = rcp (new dist_tsqr_type);
00146   distTsqr->init (scalarMessBase);
00147 
00148   return rcp (new tsqr_type (seqTsqr, distTsqr));
00149       }
00150 
00151     public:
00152 
00162       static void
00163       run (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00164      const Teuchos::RCP<node_type>& node,
00165      const Teuchos::RCP<Teuchos::ParameterList>& testParams,
00166      std::vector<int>& randomSeed)
00167       {
00168   using std::cerr;
00169   using std::cout;
00170   using std::endl;
00171   using Teuchos::arcp;
00172   using Teuchos::ParameterList;
00173   using Teuchos::parameterList;
00174   using Teuchos::RCP;
00175   using Teuchos::rcp;
00176   using Teuchos::rcp_const_cast;
00177   using Teuchos::rcp_implicit_cast;
00178   typedef Matrix<ordinal_type, scalar_type> matrix_type;
00179   typedef MatView<ordinal_type, scalar_type> view_type;
00180   typedef typename tsqr_type::FactorOutput factor_output_type;
00181 
00182   const int myRank = Teuchos::rank (*comm);
00183   const int numProcs = Teuchos::size (*comm);
00184 
00185   // Construct TSQR implementation instance.
00186   RCP<tsqr_type> tsqr = getTsqr (testParams, node, comm);
00187 
00188   // Fetch test parameters from the input parameter list.
00189   const ordinal_type numRowsLocal = testParams->get<ordinal_type> ("numRowsLocal");
00190   const ordinal_type numCols = testParams->get<ordinal_type> ("numCols"); 
00191   const int numCores = testParams->get<int> ("numCores");
00192   const bool contiguousCacheBlocks = testParams->get<bool> ("contiguousCacheBlocks");
00193   const bool testFactorExplicit = testParams->get<bool> ("testFactorExplicit");
00194   const bool testRankRevealing = testParams->get<bool> ("testRankRevealing");
00195   const bool debug = testParams->get<bool> ("debug");
00196 
00197   // Space for each node's local part of the test problem.
00198   // A_local, A_copy, and Q_local are distributed matrices, and
00199   // R is replicated on all processes sharing the communicator.
00200   matrix_type A_local (numRowsLocal, numCols);
00201   matrix_type A_copy (numRowsLocal, numCols);
00202   matrix_type Q_local (numRowsLocal, numCols);
00203   matrix_type R (numCols, numCols);
00204 
00205   // Start out by filling the test problem with zeros.
00206   typedef Teuchos::ScalarTraits<scalar_type> STS;
00207   A_local.fill (STS::zero());
00208   A_copy.fill (STS::zero());
00209   Q_local.fill (STS::zero());
00210   R.fill (STS::zero());
00211 
00212   // Create some reasonable singular values for the test problem:
00213   // 1, 1/2, 1/4, 1/8, ...
00214   typedef typename STS::magnitudeType magnitude_type;
00215   std::vector<magnitude_type> singularValues (numCols);
00216   typedef Teuchos::ScalarTraits<magnitude_type> STM;
00217   {
00218     const magnitude_type scalingFactor = STM::one() + STM::one();
00219     magnitude_type curVal = STM::one();
00220     typedef typename std::vector<magnitude_type>::iterator iter_type;
00221     for (iter_type it = singularValues.begin(); 
00222          it != singularValues.end(); ++it)
00223       {
00224         *it = curVal;
00225         curVal = curVal / scalingFactor;
00226       }
00227   }
00228 
00229   // Construct a normal(0,1) pseudorandom number generator with
00230   // the given random seed.
00231   using TSQR::Random::NormalGenerator;
00232   typedef NormalGenerator<ordinal_type, scalar_type> generator_type;
00233   generator_type gen (randomSeed);
00234 
00235   // We need a Messenger for Ordinal-type data, so that we can
00236   // build a global random test matrix.
00237   RCP<MessengerBase<ordinal_type> > ordinalMessenger = 
00238     rcp_implicit_cast<MessengerBase<ordinal_type> > (rcp (new TeuchosMessenger<ordinal_type> (comm)));
00239 
00240   // We also need a Messenger for Scalar-type data.  The TSQR
00241   // implementation already constructed one, but it's OK to
00242   // construct another one; TeuchosMessenger is just a thin
00243   // wrapper over the Teuchos::Comm object.
00244   RCP<MessengerBase<scalar_type> > scalarMessenger =
00245     rcp_implicit_cast<MessengerBase<scalar_type> > (rcp (new TeuchosMessenger<scalar_type> (comm)));
00246 
00247   {
00248     // Generate a global distributed matrix (whose part local to
00249     // this node is in A_local) with the given singular values.
00250     // This part has O(P) communication for P MPI processes.
00251     using TSQR::Random::randomGlobalMatrix;
00252     // Help the C++ compiler with type inference.
00253     view_type A_local_view (A_local.nrows(), A_local.ncols(), A_local.get(), A_local.lda());
00254     const magnitude_type* const singVals = (numCols == 0) ? NULL : &singularValues[0];
00255     randomGlobalMatrix<view_type, generator_type> (&gen, A_local_view, singVals,
00256                ordinalMessenger.getRawPtr(),
00257                scalarMessenger.getRawPtr());
00258   }
00259   // Save the pseudorandom number generator's seed for any later
00260   // tests.  The generator keeps its own copy of the seed and
00261   // updates it internally, so we have to ask for its copy.
00262   gen.getSeed (randomSeed);
00263 
00264   // If specified in the test parameters, rearrange cache blocks
00265   // in the copy.  Otherwise, just copy the test problem into
00266   // A_copy.  The factorization overwrites the input matrix, so
00267   // we have to make a copy in order to validate the final
00268   // result.
00269   if (contiguousCacheBlocks)
00270     {
00271       tsqr->cache_block (numRowsLocal, numCols, A_copy.get(), 
00272              A_local.get(), A_local.lda());
00273       if (debug)
00274         {
00275     Teuchos::barrier (*comm);
00276     if (myRank == 0)
00277       cerr << "-- Finished Tsqr::cache_block" << endl;
00278         }
00279     }
00280   else
00281     A_copy.copy (A_local);
00282 
00283   // "factorExplicit" is an alternate, hopefully faster way of
00284   // factoring the matrix, when only the explicit Q factor is
00285   // wanted.
00286   typedef Kokkos::MultiVector<scalar_type, node_type> KMV;
00287   if (testFactorExplicit)
00288     {
00289       KMV A_copy_view (node);
00290       A_copy_view.initializeValues (static_cast<size_t> (A_copy.nrows()), 
00291             static_cast<size_t> (A_copy.ncols()), 
00292             arcp (A_copy.get(), 0, A_copy.nrows()*A_copy.ncols(), false), // non-owning ArrayRCP
00293             static_cast<size_t> (A_copy.lda()));
00294       KMV Q_view (node);
00295       Q_view.initializeValues (static_cast<size_t> (Q_local.nrows()), 
00296              static_cast<size_t> (Q_local.ncols()),
00297              arcp (Q_local.get(), 0, Q_local.nrows()*Q_local.ncols(), false), // non-owning ArrayRCP
00298              static_cast<size_t> (Q_local.lda()));
00299       Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> 
00300         R_view (Teuchos::View, R.get(), R.lda(), R.nrows(), R.ncols());
00301 
00302       tsqr->factorExplicit (A_copy_view, Q_view, R_view,
00303           contiguousCacheBlocks);
00304       if (debug)
00305         {
00306     Teuchos::barrier (*comm);
00307     if (myRank == 0)
00308       cerr << "-- Finished Tsqr::factorExplicit" << endl;
00309         }
00310     }
00311   else
00312     {
00313       // Factor the (copy of the) matrix.
00314       factor_output_type factorOutput = 
00315         tsqr->factor (numRowsLocal, numCols, A_copy.get(), A_copy.lda(), 
00316           R.get(), R.lda(), contiguousCacheBlocks);
00317       if (debug)
00318         {
00319     Teuchos::barrier (*comm);
00320     if (myRank == 0)
00321       cerr << "-- Finished Tsqr::factor" << endl;
00322         }
00323       // Compute the explicit Q factor in Q_local.
00324       tsqr->explicit_Q (numRowsLocal, numCols, A_copy.get(), A_copy.lda(), 
00325             factorOutput, numCols, Q_local.get(), Q_local.lda(),
00326             contiguousCacheBlocks);
00327       if (debug)
00328         {
00329     Teuchos::barrier (*comm);
00330     if (myRank == 0)
00331       cerr << "-- Finished Tsqr::explicit_Q" << endl;
00332         }
00333     }
00334 
00335   // Optionally, test rank-revealing capability.  We do this
00336   // before un-cache-blocking the explicit Q factor, since
00337   // revealRank can work with contiguous cache blocks, and
00338   // modifies the Q factor if the matrix doesn't have full
00339   // column rank.
00340   if (testRankRevealing)
00341     {
00342       KMV Q_view (node);
00343       Q_view.initializeValues (static_cast<size_t> (Q_local.nrows()), 
00344              static_cast<size_t> (Q_local.ncols()),
00345              arcp (Q_local.get(), 0, Q_local.nrows()*Q_local.ncols(), false), // non-owning ArrayRCP
00346              static_cast<size_t> (Q_local.lda()));
00347       Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> 
00348         R_view (Teuchos::View, R.get(), R.lda(), R.nrows(), R.ncols());
00349       // If 2^{# columns} > machine precision, then our choice
00350       // of singular values will make the smallest singular
00351       // value < machine precision.  In that case, the SVD can't
00352       // promise it will distinguish between tiny and zero.  If
00353       // the number of columns is less than that, we can use a
00354       // tolerance of zero to test the purported rank with the
00355       // actual numerical rank.
00356       const magnitude_type tol = STM::zero();
00357       const ordinal_type rank = 
00358         tsqr->revealRank (Q_view, R_view, tol, contiguousCacheBlocks);
00359 
00360       magnitude_type two_to_the_numCols = STM::one();
00361       for (int k = 0; k < numCols; ++k)
00362         {
00363     const magnitude_type two = STM::one() + STM::one();
00364     two_to_the_numCols *= two;
00365         }
00366       // Throw in a factor of 10, just for more tolerance of
00367       // rounding error (so the test only fails if something is
00368       // really broken).
00369       if (two_to_the_numCols > magnitude_type(10) * STM::eps())
00370         {
00371     TEUCHOS_TEST_FOR_EXCEPTION(rank != numCols, std::logic_error,
00372            "The matrix of " << numCols << " columns "
00373            "should have full numerical rank, but Tsqr "
00374            "reports that it has rank " << rank << ".  "
00375            "Please report this bug to the Kokkos "
00376            "developers.");
00377     if (debug)
00378       {
00379         Teuchos::barrier (*comm);
00380         if (myRank == 0)
00381           cerr << "-- Tested rank-revealing capability" << endl;
00382       }
00383         }
00384       else
00385         {
00386     if (debug)
00387       {
00388         Teuchos::barrier (*comm);
00389         if (myRank == 0)
00390           cerr << "-- Not testing rank-revealing capability; too many columns" << endl;
00391       }
00392         }
00393     }
00394   // "Un"-cache-block the output, if contiguous cache blocks
00395   // were used.  This is only necessary because global_verify()
00396   // doesn't currently support contiguous cache blocks.
00397   if (contiguousCacheBlocks)
00398     {
00399       // We can use A_copy as scratch space for
00400       // un-cache-blocking Q_local, since we're done using
00401       // A_copy for other things.
00402       tsqr->un_cache_block (numRowsLocal, numCols, A_copy.get(), 
00403           A_copy.lda(), Q_local.get());
00404       // Overwrite Q_local with the un-cache-blocked Q factor.
00405       Q_local.copy (A_copy);
00406       if (debug)
00407         {
00408     Teuchos::barrier (*comm);
00409     if (myRank == 0)
00410       cerr << "-- Finished Tsqr::un_cache_block" << endl;
00411         }
00412     }
00413 
00414   // Test accuracy of the factorization.
00415   const std::vector<magnitude_type> results = 
00416     global_verify (numRowsLocal, numCols, A_local.get(), A_local.lda(),
00417        Q_local.get(), Q_local.lda(), R.get(), R.lda(), 
00418        scalarMessenger.getRawPtr());
00419   if (debug)
00420     {
00421       Teuchos::barrier (*comm);
00422       if (myRank == 0)
00423         cerr << "-- Finished global_verify" << endl;
00424     }
00425 
00426   // Print the results on Proc 0.
00427   if (myRank == 0)
00428     {
00429       if (testParams->get<bool> ("printFieldNames"))
00430         {
00431     cout << "%"
00432          << "method"
00433          << ",scalarType"
00434          << ",numRowsLocal"
00435          << ",numCols"
00436          << ",numProcs"
00437          << ",numCores"
00438          << ",cacheSizeHint"
00439          << ",contiguousCacheBlocks"
00440          << ",absFrobResid"
00441          << ",absFrobOrthog"
00442          << ",frobA" << endl;
00443     // We don't need to print field names again for the other tests, so set the test parameters accordingly.
00444     testParams->set ("printFieldNames", false);
00445         }
00446       if (testParams->get<bool> ("printResults"))
00447         {
00448     cout << "Tsqr"
00449          << "," << Teuchos::TypeNameTraits<scalar_type>::name()
00450          << "," << numRowsLocal
00451          << "," << numCols
00452          << "," << numProcs
00453          << "," << numCores
00454          << "," << tsqr->cache_size_hint()
00455          << "," << contiguousCacheBlocks 
00456          << "," << results[0] 
00457          << "," << results[1]
00458          << "," << results[2]
00459          << endl;
00460         }
00461     } // if (myRank == 0)
00462 
00463   // If requested, check accuracy and fail if results are not
00464   // sufficiently accurate.
00465   if (testParams->get<bool> ("failIfInaccurate"))
00466     {
00467       // Avoid overflow of the local Ordinal type, by casting
00468       // first to a floating-point type.
00469       const magnitude_type dimsProd = magnitude_type(numRowsLocal) * 
00470         magnitude_type(numProcs) * magnitude_type(numCols*numCols);
00471 
00472       // Relative residual error is ||A-Q*R|| / ||A||, or just
00473       // ||A-Q*R|| if ||A|| == 0.  (The result had better be
00474       // zero in the latter case.)  A reasonable error bound
00475       // should incorporate the dimensions of the matrix, since
00476       // this indicates the amount of rounding error.  Square
00477       // root of the matrix dimensions is an old heuristic from
00478       // Wilkinson or perhaps even an earlier source.  We
00479       // include a factor of 10 so that the test won't fail
00480       // unless there is a really good reason.
00481       const magnitude_type relResidBound = 
00482         magnitude_type(10) * STM::squareroot(dimsProd) * STM::eps();
00483 
00484       // Orthogonality of the matrix should not depend on the
00485       // matrix dimensions, if we measure in the 2-norm.
00486       // However, we are measuring in the Frobenius norm, so
00487       // it's appropriate to multiply eps by the number of
00488       // entries in the matrix for which we compute the
00489       // Frobenius norm.  We include a factor of 10 for the same
00490       // reason as mentioned above.
00491       const magnitude_type orthoBound = 
00492         magnitude_type(10*numCols*numCols) * STM::eps();
00493 
00494       // Avoid division by zero.
00495       const magnitude_type relResidError = 
00496         results[0] / (results[2] == STM::zero() ? STM::one() : results[2]);
00497       TEUCHOS_TEST_FOR_EXCEPTION(relResidError > relResidBound, TsqrInaccurate,
00498              "Full Tsqr (SequentialTsqr + DistTsqr) has an "
00499              "inaccurate relative residual ||A - QR||_F"
00500              << (results[2] == STM::zero() ? " / ||A||_F" : "")
00501              << " = " << relResidError << ", which is greater"
00502              " than the bound " << relResidBound << " by a "
00503              "factor of " << relResidError / relResidBound
00504              << ".");
00505       const magnitude_type orthoError = results[1];
00506       TEUCHOS_TEST_FOR_EXCEPTION(orthoError > orthoBound, TsqrInaccurate,
00507              "Full Tsqr (SequentialTsqr + DistTsqr) has an "
00508              "inaccurate orthogonality measure ||I - Q^* Q||"
00509              "_F" << results[1] << " = " << orthoError 
00510              << ", which is greater than the bound " 
00511              << orthoBound << " by a factor of " 
00512              << orthoError / orthoBound << ".");
00513     } // if (the tests should fail on inaccuracy)
00514       }
00515     };
00516 
00536     template<class TypeListType>
00537     class FullTsqrVerifierCallerImpl {
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<node_type>& node,
00544      const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00545      std::vector<int>& randomSeed);
00546     };
00547 
00548     //
00549     // Partial specialization for Cons<CarType, CdrType>.
00550     //
00551     template<class CarType, class CdrType>
00552     class FullTsqrVerifierCallerImpl<TSQR::Test::Cons<CarType, CdrType> > {
00553     public:
00554       typedef Kokkos::SerialNode node_type; 
00555 
00556       static void
00557       run (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00558      const Teuchos::RCP<node_type>& node,
00559      const Teuchos::RCP<Teuchos::ParameterList>& testParams, 
00560      std::vector<int>& randomSeed)
00561       {
00562   typedef CarType car_type;
00563   typedef CdrType cdr_type;
00564   FullTsqrVerifier<car_type>::run (comm, node, testParams, randomSeed);
00565   FullTsqrVerifierCallerImpl<cdr_type>::run (comm, node, testParams, randomSeed);
00566       }
00567     };
00568 
00569     //
00570     // Full specialization for NullCons.
00571     //
00572     template<>
00573     class FullTsqrVerifierCallerImpl<TSQR::Test::NullCons> {
00574     public:
00575       typedef Kokkos::SerialNode node_type; 
00576 
00577       static void
00578       run (const Teuchos::RCP<const Teuchos::Comm<int> >&,
00579      const Teuchos::RCP<node_type>&,
00580      const Teuchos::RCP<Teuchos::ParameterList>&, 
00581      std::vector<int>&)
00582       {
00583   // We're at the end of the type list, so do nothing.
00584       }
00585     };
00586 
00594     class FullTsqrVerifierCaller {
00595     public:
00598       typedef Kokkos::SerialNode node_type;
00599 
00609       typedef int ordinal_type;
00610 
00616       Teuchos::RCP<const Teuchos::ParameterList>
00617       getValidParameterList () const
00618       {
00619   using Teuchos::ParameterList;
00620   using Teuchos::parameterList;
00621   using Teuchos::RCP;
00622 
00623   RCP<ParameterList> plist = parameterList ("FullTsqrVerifier");
00624 
00625   const size_t cacheSizeHint = 0;
00626   const int numCores = 1;
00627   const ordinal_type numRowsLocal = 100;
00628   const ordinal_type numCols = 10;
00629   const bool contiguousCacheBlocks = false;
00630   const bool testFactorExplicit = true;
00631   const bool testRankRevealing = true;
00632   const bool printFieldNames = true;
00633   const bool printResults = true;
00634   const bool failIfInaccurate = true;
00635   const bool debug = false;
00636 
00637   // Parameters for configuring Tsqr itself.
00638   plist->set ("cacheSizeHint", cacheSizeHint, 
00639         "Cache size hint in bytes.  "
00640         "Zero means TSQR picks a reasonable default.");
00641   plist->set ("numCores", numCores,
00642         "Number of partition(s) to use for TbbTsqr (if "
00643         "applicable).  Must be a positive integer.");
00644 
00645   // Parameters for testing Tsqr.
00646   plist->set ("numRowsLocal", numRowsLocal, 
00647         "Number of rows per (MPI) process in the test matrix.  "
00648         "Must be >= the number of columns.");
00649   plist->set ("numCols", numCols, 
00650         "Number of columns in the test matrix.");
00651   plist->set ("contiguousCacheBlocks", contiguousCacheBlocks, 
00652         "Whether to test the factorization with contiguously "
00653         "stored cache blocks.");
00654   plist->set ("testFactorExplicit", testFactorExplicit, 
00655         "Whether to test TSQR's factorExplicit() (a hopefully "
00656         "faster path than calling factor() and explicit_Q() in "
00657         "sequence).");
00658   plist->set ("testRankRevealing", testRankRevealing, 
00659         "Whether to test TSQR's rank-revealing capability.");
00660   plist->set ("printFieldNames", printFieldNames, 
00661         "Whether to print field names (this is only done once, "
00662         "for all Scalar types tested).");
00663   plist->set ("printResults", printResults, 
00664         "Whether to print test results.");
00665   plist->set ("failIfInaccurate", failIfInaccurate,
00666         "Whether to fail the test if the factorization "
00667         "is not sufficiently accurate.");
00668   plist->set ("debug", debug, 
00669         "Whether to print debugging output.");
00670   return plist;
00671       }
00672 
00684       template<class TypeListType>
00685       void 
00686       run (const Teuchos::RCP<Teuchos::ParameterList>& testParams)
00687       {
00688   // Using a class with a static method is a way to implement
00689   // "partial specialization of function templates" (which by
00690   // itself is not allowed in C++).
00691   typedef FullTsqrVerifierCallerImpl<TypeListType> impl_type;
00692   impl_type::run (comm_, node_, testParams, randomSeed_);
00693       }
00694 
00695 
00710       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00711             const Teuchos::RCP<node_type>& node,
00712             const std::vector<int>& randomSeed) :
00713   comm_ (comm), node_ (node), randomSeed_ (validateRandomSeed (randomSeed))
00714       {}
00715 
00723       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm) :
00724   comm_ (comm),
00725   node_ (getNode<node_type> (getValidNodeParameters<node_type> ())),
00726   randomSeed_ (defaultRandomSeed ())
00727       {}
00728 
00738       FullTsqrVerifierCaller (const Teuchos::RCP<const Teuchos::Comm<int> >& comm,
00739             const Teuchos::RCP<node_type>& node) :
00740   comm_ (comm),
00741   node_ (node),
00742   randomSeed_ (defaultRandomSeed ())
00743       {}
00744 
00746       static std::vector<int> 
00747       validateRandomSeed (const std::vector<int>& seed) 
00748       {
00749   TEUCHOS_TEST_FOR_EXCEPTION(seed.size() < 4, std::invalid_argument,
00750          "Invalid random seed: Need an array of four integers.");
00751   for (std::vector<int>::size_type k = 0; k < seed.size(); ++k)
00752     {
00753       TEUCHOS_TEST_FOR_EXCEPTION(seed[k] < 0 || seed[k] > 4095,
00754              std::invalid_argument,
00755              "Invalid random seed: Each of the four integers must be in [0, 4095].");
00756     }
00757   TEUCHOS_TEST_FOR_EXCEPTION(seed[3] % 2 != 1, std::invalid_argument,
00758          "Invalid random seed: The last of the four integers must be odd.");
00759   return seed;
00760       }
00761 
00763       static std::vector<int> 
00764       defaultRandomSeed () 
00765       {
00766   std::vector<int> seed (4);
00767   seed[0] = 0;
00768   seed[1] = 0;
00769   seed[2] = 0;
00770   seed[3] = 1;
00771   return seed;
00772       }
00773 
00774     private:
00779       Teuchos::RCP<const Teuchos::Comm<int> > comm_;
00780 
00782       Teuchos::RCP<node_type> node_;
00783 
00789       std::vector<int> randomSeed_;
00790     };
00791 
00792   } // namespace Test
00793 } // namespace TSQR
00794 
00795 #endif // __TSQR_Test_FullTsqrTest_hpp
00796 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends