Tpetra parallel linear algebra  Version of the Day
Tpetra_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 __Tpetra_TsqrAdaptor_hpp
43 #define __Tpetra_TsqrAdaptor_hpp
44 
48 
49 #include <Tpetra_ConfigDefs.hpp>
50 
51 #ifdef HAVE_TPETRA_TSQR
52 # include <Tsqr_NodeTsqrFactory.hpp> // create intranode TSQR object
53 # include <Tsqr.hpp> // full (internode + intranode) TSQR
54 # include <Tsqr_DistTsqr.hpp> // internode TSQR
55 // Subclass of TSQR::MessengerBase, implemented using Teuchos
56 // communicator template helper functions
57 # include <Tsqr_TeuchosMessenger.hpp>
58 # include <Tpetra_MultiVector.hpp>
59 # include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
60 # include <stdexcept>
61 
62 
63 namespace Tpetra {
64 
87  template<class MV>
88  class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
89  public:
90  typedef typename MV::scalar_type scalar_type;
91  typedef typename MV::local_ordinal_type ordinal_type;
92  typedef typename MV::node_type node_type;
93  typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
94  typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
95 
96  private:
97  //typedef TSQR::MatView<ordinal_type, scalar_type> matview_type;
98  typedef TSQR::NodeTsqrFactory<node_type, scalar_type, ordinal_type> node_tsqr_factory_type;
99  typedef typename node_tsqr_factory_type::node_tsqr_type node_tsqr_type;
100  typedef TSQR::DistTsqr<ordinal_type, scalar_type> dist_tsqr_type;
101  typedef TSQR::Tsqr<ordinal_type, scalar_type, node_tsqr_type> tsqr_type;
102 
103  public:
110  TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
111  nodeTsqr_ (new node_tsqr_type),
112  distTsqr_ (new dist_tsqr_type),
113  tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
114  ready_ (false)
115  {
116  setParameterList (plist);
117  }
118 
120  TsqrAdaptor () :
121  nodeTsqr_ (new node_tsqr_type),
122  distTsqr_ (new dist_tsqr_type),
123  tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
124  ready_ (false)
125  {
126  setParameterList (Teuchos::null);
127  }
128 
130  Teuchos::RCP<const Teuchos::ParameterList>
131  getValidParameters () const
132  {
133  using Teuchos::RCP;
134  using Teuchos::rcp;
135  using Teuchos::ParameterList;
136  using Teuchos::parameterList;
137 
138  if (defaultParams_.is_null()) {
139  RCP<ParameterList> params = parameterList ("TSQR implementation");
140  params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
141  params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
142  defaultParams_ = params;
143  }
144  return defaultParams_;
145  }
146 
172  void
173  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
174  {
175  using Teuchos::ParameterList;
176  using Teuchos::parameterList;
177  using Teuchos::RCP;
178  using Teuchos::sublist;
179 
180  RCP<ParameterList> params = plist.is_null() ?
181  parameterList (*getValidParameters ()) : plist;
182  nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
183  distTsqr_->setParameterList (sublist (params, "DistTsqr"));
184 
185  this->setMyParamList (params);
186  }
187 
209  void
210  factorExplicit (MV& A,
211  MV& Q,
212  dense_matrix_type& R,
213  const bool forceNonnegativeDiagonal=false)
214  {
215  typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
216 
217  prepareTsqr (Q); // Finish initializing TSQR.
218  KMV A_view = getNonConstView (A);
219  KMV Q_view = getNonConstView (Q);
220  tsqr_->factorExplicit (A_view, Q_view, R, false,
221  forceNonnegativeDiagonal);
222  }
223 
254  int
255  revealRank (MV& Q,
256  dense_matrix_type& R,
257  const magnitude_type& tol)
258  {
259  typedef KokkosClassic::MultiVector<scalar_type, node_type> KMV;
260 
261  prepareTsqr (Q); // Finish initializing TSQR.
262 
263  // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
264  // to make sure it is the same communicator as the one we are
265  // using in our dist_tsqr_type implementation.
266  KMV Q_view = getNonConstView (Q);
267  return tsqr_->revealRank (Q_view, R, tol, false);
268  }
269 
270  private:
272  Teuchos::RCP<node_tsqr_type> nodeTsqr_;
273 
275  Teuchos::RCP<dist_tsqr_type> distTsqr_;
276 
278  Teuchos::RCP<tsqr_type> tsqr_;
279 
281  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
282 
284  bool ready_;
285 
306  void
307  prepareTsqr (const MV& mv)
308  {
309  if (! ready_) {
310  prepareDistTsqr (mv);
311  prepareNodeTsqr (mv);
312  ready_ = true;
313  }
314  }
315 
319  void
320  prepareNodeTsqr (const MV& mv)
321  {
322  node_tsqr_factory_type::prepareNodeTsqr (nodeTsqr_, mv.getMap()->getNode());
323  }
324 
331  void
332  prepareDistTsqr (const MV& mv)
333  {
334  using Teuchos::RCP;
335  using Teuchos::rcp_implicit_cast;
336  typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
337  typedef TSQR::MessengerBase<scalar_type> base_mess_type;
338 
339  RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
340  RCP<mess_type> mess (new mess_type (comm));
341  RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
342  distTsqr_->init (messBase);
343  }
344 
357  static KokkosClassic::MultiVector<scalar_type, node_type>
358  getNonConstView (MV& A)
359  {
360  // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
361  // storage of A uses nonconstant stride internally. We would
362  // have to copy and pack into a matrix with constant stride, and
363  // then unpack on exit. For now we choose just to raise an
364  // exception.
365  TEUCHOS_TEST_FOR_EXCEPTION(
366  ! A.isConstantStride(), std::invalid_argument,
367  "Tpetra::TsqrAdaptor::getNonConstView: TSQR does not currently "
368  "support Tpetra::MultiVector inputs that do not have constant "
369  "stride.");
370  return A.getLocalMV ();
371  }
372  };
373 
374 } // namespace Tpetra
375 
376 #endif // HAVE_TPETRA_TSQR
377 
378 #endif // __Tpetra_TsqrAdaptor_hpp
379 
Namespace Tpetra contains the class and methods constituting the Tpetra library.
KokkosClassic::DefaultNode::DefaultNodeType node_type
Default value of Node template parameter.
double scalar_type
Default value of Scalar template parameter.