Anasazi  Version of the Day
TsqrTpetraTest.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2010) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef __TSQR_Trilinos_TsqrTpetraTest_hpp
30 #define __TSQR_Trilinos_TsqrTpetraTest_hpp
31 
32 #include "AnasaziTpetraAdapter.hpp" // sic (not "-or")
33 #include "Tsqr_Random_NormalGenerator.hpp"
35 #include "Teuchos_Tuple.hpp"
36 #include "Teuchos_Time.hpp"
37 #include <sstream>
38 #include <stdexcept>
39 #include <vector>
40 
41 
42 namespace TSQR {
43  namespace Trilinos {
44  namespace Test {
45  using Teuchos::RCP;
46  using Teuchos::Tuple;
47 
48  template< class S, class LO, class GO, class Node >
49  class TpetraTsqrTest {
50  public:
51  typedef S scalar_type;
52  typedef LO local_ordinal_type;
53  typedef GO global_ordinal_type;
54  typedef Node node_type;
55 
56  typedef typename TSQR::ScalarTraits< S >::magnitude_type magnitude_type;
57  typedef TSQR::Trilinos::TsqrTpetraAdaptor< S, LO, GO, Node > adaptor_type;
58 
59  TpetraTsqrTest (const Tpetra::global_size_t nrowsGlobal,
60  const size_t ncols,
61  const Teuchos::RCP< const Teuchos::Comm<int> >& comm,
62  const Teuchos::RCP< Node >& node,
63  const Teuchos::ParameterList& params) :
64  results_ (magnitude_type(0), magnitude_type(0))
65  {
66  using Teuchos::Tuple;
68  typedef typename adaptor_type::factor_output_type factor_output_type;
69  typedef Teuchos::SerialDenseMatrix< LO, S > matrix_type;
70 
71  bool contiguousCacheBlocks = false;
72  try {
73  contiguousCacheBlocks = params.get<bool>("contiguousCacheBlocks");
74  } catch (InvalidParameter&) {
75  contiguousCacheBlocks = false;
76  }
77 
78  triple_type testProblem =
79  makeTestProblem (nrowsGlobal, ncols, comm, node, params);
80  // A is already filled in with the test problem. A_copy and
81  // Q are just multivectors with the same layout as A.
82  // A_copy will be used for temporary storage, and Q will
83  // store the (explicitly represented) Q factor output. R
84  // will store the R factor output.
85  RCP< MV > A = testProblem[0];
86  RCP< MV > A_copy = testProblem[1];
87  RCP< MV > Q = testProblem[2];
88  matrix_type R (ncols, ncols);
89 
90  // Adaptor uses one of the multivectors only to reference
91  // the underlying communicator object.
92  adaptor_type adaptor (*A, params);
93  if (contiguousCacheBlocks)
94  adaptor.cacheBlock (*A, *A_copy);
95 
96  factor_output_type factorOutput =
97  adaptor.factor (*A_copy, R, contiguousCacheBlocks);
98  adaptor.explicitQ (*A_copy, factorOutput, *Q, contiguousCacheBlocks);
99  if (contiguousCacheBlocks)
100  {
101  // Use A_copy as temporary storage for un-cache-blocking
102  // Q. Tpetra::MultiVector objects copy deeply.
103  *A_copy = *Q;
104  adaptor.unCacheBlock (*A_copy, *Q);
105  }
106  results_ = adaptor.verify (*A, *Q, R);
107  }
108 
111  std::pair< magnitude_type, magnitude_type >
112  getResults() const
113  {
114  return results_;
115  }
116 
117  private:
118  typedef Tpetra::MultiVector< S, LO, GO, Node > MV;
119  typedef Teuchos::Tuple< RCP< MV >, 3 > triple_type;
120  typedef Teuchos::RCP< const Teuchos::Comm<int> > comm_ptr;
121  typedef Tpetra::Map< LO, GO, Node > map_type;
122  typedef Teuchos::RCP< const map_type > map_ptr;
123  typedef TSQR::Random::NormalGenerator< LO, S > normalgen_type;
124  typedef Teuchos::RCP< Node > node_ptr;
125 
130  std::pair< magnitude_type, magnitude_type > results_;
131 
140  static map_ptr
141  makeMap (const Tpetra::global_size_t nrowsGlobal,
142  const comm_ptr& comm,
143  const node_ptr& node)
144  {
145  using Tpetra::createUniformContigMapWithNode;
146  return createUniformContigMapWithNode< LO, GO, Node > (nrowsGlobal,
147  comm, node);
148  }
149 
152  static RCP< MV >
153  makeMultiVector (const map_ptr& map,
154  const size_t ncols)
155  {
156  // "false" means "don't fill with zeros"; we'll fill it in
157  // fillTpetraMultiVector().
158  return Teuchos::rcp (new MV (map, ncols, false));
159  }
160 
163  static void
164  fillMultiVector (const RCP< MV >& mv,
165  const RCP< normalgen_type >& pGen)
166  {
167  using TSQR::Trilinos::TpetraRandomizer;
168  typedef TpetraRandomizer< S, LO, GO, Node, normalgen_type > randomizer_type;
169 
170  const LO ncols = mv->getNumVectors();
171  std::vector< S > singular_values (ncols);
172  if (ncols > 0)
173  {
174  singular_values[0] = S(1);
175  for (LO k = 1; k < ncols; ++k)
176  singular_values[k] = singular_values[k-1] / S(2);
177  }
178  randomizer_type randomizer (*mv, pGen);
179  randomizer.randomMultiVector (*mv, &singular_values[0]);
180  }
181 
182  static triple_type
183  makeTestProblem (const Tpetra::global_size_t nrowsGlobal,
184  const size_t ncols,
185  const comm_ptr& comm,
186  const node_ptr& node,
187  const Teuchos::ParameterList& params)
188  {
189  using TSQR::Trilinos::TpetraMessenger;
190  using TSQR::MessengerBase;
191 
192  map_ptr map = makeMap (nrowsGlobal, comm, node);
193  RCP< MV > A = makeMultiVector (map, ncols);
194  RCP< MV > A_copy = makeMultiVector (map, ncols);
195  RCP< MV > Q = makeMultiVector (map, ncols);
196 
197  // Fill A with the random test problem
198  RCP< normalgen_type > pGen (new normalgen_type);
199  fillMultiVector (A, pGen);
200 
201  return Teuchos::tuple (A, A_copy, Q);
202  }
203  };
204 
205  } // namespace Test
206  } // namespace Trilinos
207 } // namespace TSQR
208 
209 #endif // __TSQR_Trilinos_TsqrTpetraTest_hpp
T & get(const std::string &name, T def_value)
Partial specialization of Anasazi::MultiVecTraits and Anasazi::OperatorTraits for Tpetra objects...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)