Anasazi Version of the Day
TsqrTpetraTest.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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_Trilinos_TsqrTpetraTest_hpp
00030 #define __TSQR_Trilinos_TsqrTpetraTest_hpp
00031 
00032 #include "Teuchos_Tuple.hpp"
00033 #include "Teuchos_ParameterList.hpp"
00034 
00035 #include "Kokkos_ConfigDefs.hpp"
00036 #include "Kokkos_NodeHelpers.hpp"
00037 #include "Kokkos_SerialNode.hpp"
00038 #ifdef HAVE_KOKKOSCLASSIC_TBB
00039 #include "Kokkos_TBBNode.hpp"
00040 #endif // HAVE_KOKKOSCLASSIC_TBB
00041 
00042 #include "Teuchos_Time.hpp"
00043 #include "AnasaziTpetraAdapter.hpp" // sic (not "-or")
00044 #include "Tsqr_Random_NormalGenerator.hpp"
00045 
00046 #include <sstream>
00047 #include <stdexcept>
00048 #include <vector>
00049 
00052 
00053 namespace TSQR { 
00054   namespace Trilinos { 
00055     namespace Test {
00056       using Teuchos::RCP;
00057       using Teuchos::Tuple;
00058 
00059       template< class S, class LO, class GO, class Node >
00060       class TpetraTsqrTest {
00061       public:
00062   typedef S scalar_type;
00063   typedef LO local_ordinal_type;
00064   typedef GO global_ordinal_type;
00065   typedef Node node_type;
00066 
00067   typedef typename TSQR::ScalarTraits< S >::magnitude_type magnitude_type;
00068   typedef TSQR::Trilinos::TsqrTpetraAdaptor< S, LO, GO, Node > adaptor_type;
00069 
00070   TpetraTsqrTest (const Tpetra::global_size_t nrowsGlobal,
00071       const size_t ncols,
00072       const Teuchos::RCP< const Teuchos::Comm<int> >& comm,
00073       const Teuchos::RCP< Node >& node,
00074       const Teuchos::ParameterList& params) :
00075     results_ (magnitude_type(0), magnitude_type(0))
00076   {
00077     using Teuchos::Tuple;
00078     using Teuchos::Exceptions::InvalidParameter;
00079     typedef typename adaptor_type::factor_output_type factor_output_type;
00080     typedef Teuchos::SerialDenseMatrix< LO, S > matrix_type;
00081 
00082     bool contiguousCacheBlocks = false;
00083     try {
00084       contiguousCacheBlocks = params.get<bool>("contiguousCacheBlocks");
00085     } catch (InvalidParameter&) {
00086       contiguousCacheBlocks = false;
00087     }
00088 
00089     triple_type testProblem = 
00090       makeTestProblem (nrowsGlobal, ncols, comm, node, params);
00091     // A is already filled in with the test problem.  A_copy and
00092     // Q are just multivectors with the same layout as A.
00093     // A_copy will be used for temporary storage, and Q will
00094     // store the (explicitly represented) Q factor output.  R
00095     // will store the R factor output.
00096     RCP< MV > A = testProblem[0]; 
00097     RCP< MV > A_copy = testProblem[1];
00098     RCP< MV > Q = testProblem[2];
00099     matrix_type R (ncols, ncols);
00100 
00101     // Adaptor uses one of the multivectors only to reference
00102     // the underlying communicator object.
00103     adaptor_type adaptor (*A, params);
00104     if (contiguousCacheBlocks)
00105       adaptor.cacheBlock (*A, *A_copy);
00106 
00107     factor_output_type factorOutput = 
00108       adaptor.factor (*A_copy, R, contiguousCacheBlocks);
00109     adaptor.explicitQ (*A_copy, factorOutput, *Q, contiguousCacheBlocks);
00110     if (contiguousCacheBlocks)
00111       {
00112         // Use A_copy as temporary storage for un-cache-blocking
00113         // Q.  Tpetra::MultiVector objects copy deeply.
00114         *A_copy = *Q;
00115         adaptor.unCacheBlock (*A_copy, *Q);
00116       }
00117     results_ = adaptor.verify (*A, *Q, R);
00118   }
00119 
00122   std::pair< magnitude_type, magnitude_type > 
00123   getResults() const 
00124   { 
00125     return results_;
00126   }
00127 
00128       private:
00129   typedef Tpetra::MultiVector< S, LO, GO, Node >   MV;
00130   typedef Teuchos::Tuple< RCP< MV >, 3 >           triple_type;
00131   typedef Teuchos::RCP< const Teuchos::Comm<int> > comm_ptr;
00132   typedef Tpetra::Map< LO, GO, Node >              map_type;
00133   typedef Teuchos::RCP< const map_type >           map_ptr;
00134   typedef TSQR::Random::NormalGenerator< LO, S >   normalgen_type;
00135   typedef Teuchos::RCP< Node >                     node_ptr;
00136 
00141   std::pair< magnitude_type, magnitude_type > results_;
00142 
00151   static map_ptr
00152   makeMap (const Tpetra::global_size_t nrowsGlobal,
00153      const comm_ptr& comm,
00154      const node_ptr& node)
00155   {
00156     using Tpetra::createUniformContigMapWithNode;
00157     return createUniformContigMapWithNode< LO, GO, Node > (nrowsGlobal, 
00158                  comm, node);
00159   }
00160 
00163   static RCP< MV >
00164   makeMultiVector (const map_ptr& map,
00165        const size_t ncols)
00166   {
00167     // "false" means "don't fill with zeros"; we'll fill it in
00168     // fillTpetraMultiVector().
00169     return Teuchos::rcp (new MV (map, ncols, false));
00170   }
00171 
00174   static void
00175   fillMultiVector (const RCP< MV >& mv,
00176        const RCP< normalgen_type >& pGen)
00177   {
00178     using TSQR::Trilinos::TpetraRandomizer;
00179     typedef TpetraRandomizer< S, LO, GO, Node, normalgen_type > randomizer_type;
00180 
00181     const LO ncols = mv->getNumVectors();
00182     std::vector< S > singular_values (ncols);
00183     if (ncols > 0)
00184       {
00185         singular_values[0] = S(1);
00186         for (LO k = 1; k < ncols; ++k)
00187     singular_values[k] = singular_values[k-1] / S(2);
00188       }
00189     randomizer_type randomizer (*mv, pGen);
00190     randomizer.randomMultiVector (*mv, &singular_values[0]);
00191   }
00192 
00193   static triple_type
00194   makeTestProblem (const Tpetra::global_size_t nrowsGlobal,
00195        const size_t ncols,
00196        const comm_ptr& comm,
00197        const node_ptr& node,
00198        const Teuchos::ParameterList& params)
00199   {
00200     using TSQR::Trilinos::TpetraMessenger;
00201     using TSQR::MessengerBase;
00202 
00203     map_ptr map = makeMap (nrowsGlobal, comm, node);
00204     RCP< MV > A = makeMultiVector (map, ncols);
00205     RCP< MV > A_copy = makeMultiVector (map, ncols);
00206     RCP< MV > Q = makeMultiVector (map, ncols);
00207 
00208     // Fill A with the random test problem
00209     RCP< normalgen_type > pGen (new normalgen_type);
00210     fillMultiVector (A, pGen);
00211 
00212     return Teuchos::tuple (A, A_copy, Q);
00213   }
00214       };
00215 
00216     } // namespace Test
00217   } // namespace Trilinos
00218 } // namespace TSQR
00219 
00220 #endif // __TSQR_Trilinos_TsqrTpetraTest_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends