Tpetra parallel linear algebra  Version of the Day
Epetra_TsqrAdaptor.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Tpetra: Templated Linear Algebra Services Package
5 // Copyright (2008) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 // @HEADER
41 
42 #ifndef __Epetra_TsqrAdaptor_hpp
43 #define __Epetra_TsqrAdaptor_hpp
44 
58 
59 #include <Tpetra_ConfigDefs.hpp>
60 
61 #if defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
62 
63 #include <Kokkos_DefaultNode.hpp> // Include minimal Kokkos Node types
64 #include <Tsqr_NodeTsqrFactory.hpp> // create intranode TSQR object
65 #include <Tsqr.hpp> // full (internode + intranode) TSQR
66 #include <Tsqr_DistTsqr.hpp> // internode TSQR
67 #include <Epetra_Comm.h>
68 // Subclass of TSQR::MessengerBase, implemented using Teuchos
69 // communicator template helper functions
70 #include <Epetra_TsqrMessenger.hpp>
71 #include <Epetra_MultiVector.h>
72 #include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
73 #include <stdexcept>
74 
75 
76 namespace Epetra {
77 
102  class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
103  public:
104  typedef Epetra_MultiVector MV;
105 
112  typedef double scalar_type;
113 
120  typedef int ordinal_type;
121 
127 
136  typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
137 
143  typedef double magnitude_type;
144 
145  private:
146  typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
147  typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
148  // Don't need a "typename" here, because there are no template
149  // parameters involved in the type definition.
150  typedef node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
151  typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
152  typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
153 
154  public:
161  TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
162  nodeTsqr_ (new node_tsqr_type),
163  distTsqr_ (new dist_tsqr_type),
164  tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
165  ready_ (false)
166  {
167  setParameterList (plist);
168  }
169 
171  TsqrAdaptor () :
172  nodeTsqr_ (new node_tsqr_type),
173  distTsqr_ (new dist_tsqr_type),
174  tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
175  ready_ (false)
176  {
177  setParameterList (Teuchos::null);
178  }
179 
180  Teuchos::RCP<const Teuchos::ParameterList>
181  getValidParameters () const
182  {
183  using Teuchos::RCP;
184  using Teuchos::rcp;
185  using Teuchos::ParameterList;
186  using Teuchos::parameterList;
187 
188  if (defaultParams_.is_null()) {
189  RCP<ParameterList> params = parameterList ("TSQR implementation");
190  params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
191  params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
192  defaultParams_ = params;
193  }
194  return defaultParams_;
195  }
196 
197  void
198  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
199  {
200  using Teuchos::ParameterList;
201  using Teuchos::parameterList;
202  using Teuchos::RCP;
203  using Teuchos::sublist;
204 
205  RCP<ParameterList> params = plist.is_null() ?
206  parameterList (*getValidParameters ()) : plist;
207  nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
208  distTsqr_->setParameterList (sublist (params, "DistTsqr"));
209 
210  this->setMyParamList (params);
211  }
212 
234  void
235  factorExplicit (MV& A,
236  MV& Q,
237  dense_matrix_type& R,
238  const bool forceNonnegativeDiagonal=false)
239  {
240  typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
241 
242  prepareTsqr (Q); // Finish initializing TSQR.
243  KMV A_view = getNonConstView (A);
244  KMV Q_view = getNonConstView (Q);
245  tsqr_->factorExplicit (A_view, Q_view, R, false,
246  forceNonnegativeDiagonal);
247  }
248 
279  int
280  revealRank (MV& Q,
281  dense_matrix_type& R,
282  const magnitude_type& tol)
283  {
284  typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
285 
286  prepareTsqr (Q); // Finish initializing TSQR.
287 
288  // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
289  // sure it is the same communicator as the one we are using in
290  // our dist_tsqr_type implementation.
291  KMV Q_view = getNonConstView (Q);
292  return tsqr_->revealRank (Q_view, R, tol, false);
293  }
294 
295  private:
297  Teuchos::RCP<node_tsqr_type> nodeTsqr_;
298 
300  Teuchos::RCP<dist_tsqr_type> distTsqr_;
301 
303  Teuchos::RCP<tsqr_type> tsqr_;
304 
306  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
307 
309  bool ready_;
310 
329  void
330  prepareTsqr (const MV& mv)
331  {
332  if (! ready_) {
333  prepareDistTsqr (mv);
334  prepareNodeTsqr (mv);
335  ready_ = true;
336  }
337  }
338 
342  void
343  prepareNodeTsqr (const MV& mv)
344  {
345  (void) mv; // Epetra objects don't have a Kokkos Node.
346 
347  // Create Node with empty ParameterList.
348  Teuchos::ParameterList plist;
349  Teuchos::RCP<node_type> node (new node_type (plist));
350  node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, node);
351  }
352 
359  void
360  prepareDistTsqr (const MV& mv)
361  {
362  using Teuchos::RCP;
363  using Teuchos::rcp;
364  using TSQR::Epetra::makeTsqrMessenger;
365  typedef TSQR::MessengerBase<scalar_type> base_mess_type;
366 
367  // If mv falls out of scope, its Epetra_Comm may become invalid.
368  // Thus, we clone the input Epetra_Comm, so that the messenger
369  // owns the object.
370  RCP<const Epetra_Comm> comm = rcp (mv.Comm().Clone());
371  RCP<base_mess_type> messBase = makeTsqrMessenger<scalar_type> (comm);
372  distTsqr_->init (messBase);
373  }
374 
387  static KokkosClassic::MultiVector<scalar_type, node_type>
388  getNonConstView (MV& A)
389  {
390  // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
391  // storage of A uses nonconstant stride internally. We would
392  // have to copy and pack into a matrix with constant stride, and
393  // then unpack on exit. For now we choose just to raise an
394  // exception.
395  TEUCHOS_TEST_FOR_EXCEPTION(! A.ConstantStride(), std::invalid_argument,
396  "TSQR does not currently support Epetra_MultiVector "
397  "inputs that do not have constant stride.");
398  const int numRows = A.MyLength();
399  const int numCols = A.NumVectors();
400  const int stride = A.Stride();
401  // A_ptr does _not_ own the data. TSQR only operates within the
402  // scope of the multivector objects on which it operates, so it
403  // doesn't need ownership of the data.
404  Teuchos::ArrayRCP<double> A_ptr (A.Values(), 0, numRows*stride, false);
405 
406  typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
407  // KMV objects want a Kokkos Node instance. Epetra objects
408  // don't have a Kokkos Node, so we make a temporary node just
409  // for the KMV.
410  //
411  // Create Node with empty ParameterList.
412  Teuchos::ParameterList plist;
413  Teuchos::RCP<node_type> node (new node_type (plist));
414  KMV A_kmv (node);
415  A_kmv.initializeValues (numRows, numCols, A_ptr, stride);
416  return A_kmv;
417  }
418  };
419 
420 } // namespace Epetra
421 
422 #endif // defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
423 
424 #endif // __Epetra_TsqrAdaptor_hpp
425 
KokkosClassic::DefaultNode::DefaultNodeType node_type
Default value of Node template parameter.
Function for wrapping Epetra_Comm in a communicator wrapper that TSQR can use.
double scalar_type
Default value of Scalar template parameter.