Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Epetra_TsqrAdaptor.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Tpetra: Templated Linear Algebra Services Package
4 //
5 // Copyright 2008 NTESS and the Tpetra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef EPETRA_TSQRADAPTOR_HPP
11 #define EPETRA_TSQRADAPTOR_HPP
12 
24 
25 #include "Tpetra_ConfigDefs.hpp"
26 
27 #if defined(TPETRA_ENABLE_DEPRECATED_CODE)
28 #if defined(TPETRA_DEPRECATED_DECLARATIONS)
29 #warning This file is deprecated due to Epetra removal and will be removed
30 #endif
31 #else
32 #error This file is deprecated due to Epetra removal and will be removed
33 #endif
34 
35 #if defined(TPETRA_ENABLE_DEPRECATED_CODE) && defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
36 
37 #include "Tsqr_NodeTsqrFactory.hpp" // create intranode TSQR object
38 #include "Tsqr.hpp" // full (internode + intranode) TSQR
39 #include "Tsqr_DistTsqr.hpp" // internode TSQR
40 #include "Epetra_Comm.h"
41 // Subclass of TSQR::MessengerBase, implemented using Teuchos
42 // communicator template helper functions
43 #include "Epetra_TsqrMessenger.hpp"
44 #include "Epetra_MultiVector.h"
45 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
46 #include <stdexcept>
47 
48 namespace Epetra {
49 
74 TPETRA_DEPRECATED_MSG("epetra removal")
75 class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
76  public:
77  typedef Epetra_MultiVector MV;
78 
85  typedef double scalar_type;
86 
93  typedef int ordinal_type;
94 
100  using device_type =
101  Kokkos::Device<Kokkos::DefaultHostExecutionSpace,
102  Kokkos::HostSpace>;
103 
112  using dense_matrix_type =
113  Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>;
114 
120  using magnitude_type = double;
121 
122  private:
123  using matview_type = TSQR::MatView<ordinal_type, scalar_type>;
124  using node_tsqr_factory_type =
125  TSQR::NodeTsqrFactory<scalar_type, ordinal_type, device_type>;
126  // Don't need a "typename" here, because there are no template
127  // parameters involved in the type definition.
128  using node_tsqr_type = TSQR::NodeTsqr<ordinal_type, scalar_type>;
129  using dist_tsqr_type = TSQR::DistTsqr<ordinal_type, scalar_type>;
130  using tsqr_type = TSQR::Tsqr<ordinal_type, scalar_type>;
131 
132  public:
139  TsqrAdaptor(const Teuchos::RCP<Teuchos::ParameterList>& plist)
140  : nodeTsqr_(node_tsqr_factory_type::getNodeTsqr())
141  , distTsqr_(new dist_tsqr_type)
142  , tsqr_(new tsqr_type(nodeTsqr_, distTsqr_))
143  , ready_(false) {
144  setParameterList(plist);
145  }
146 
148  TsqrAdaptor()
149  : nodeTsqr_(node_tsqr_factory_type::getNodeTsqr())
150  , distTsqr_(new dist_tsqr_type)
151  , tsqr_(new tsqr_type(nodeTsqr_, distTsqr_))
152  , ready_(false) {
153  setParameterList(Teuchos::null);
154  }
155 
156  Teuchos::RCP<const Teuchos::ParameterList>
157  getValidParameters() const {
158  using Teuchos::ParameterList;
159  using Teuchos::parameterList;
160  using Teuchos::RCP;
161  using Teuchos::rcp;
162 
163  if (defaultParams_.is_null()) {
164  RCP<ParameterList> params = parameterList("TSQR implementation");
165  params->set("NodeTsqr", *(nodeTsqr_->getValidParameters()));
166  params->set("DistTsqr", *(distTsqr_->getValidParameters()));
167  defaultParams_ = params;
168  }
169  return defaultParams_;
170  }
171 
172  void
173  setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& plist) {
174  using Teuchos::ParameterList;
175  using Teuchos::parameterList;
176  using Teuchos::RCP;
177  using Teuchos::sublist;
178 
179  RCP<ParameterList> params = plist.is_null() ? parameterList(*getValidParameters()) : plist;
180  nodeTsqr_->setParameterList(sublist(params, "NodeTsqr"));
181  distTsqr_->setParameterList(sublist(params, "DistTsqr"));
182 
183  this->setMyParamList(params);
184  }
185 
207  void
208  factorExplicit(MV& A,
209  MV& Q,
210  dense_matrix_type& R,
211  const bool forceNonnegativeDiagonal = false) {
212  prepareTsqr(Q); // Finish initializing TSQR.
213 
214  scalar_type* const A_ptr = A.Values();
215  scalar_type* const Q_ptr = Q.Values();
216  scalar_type* const R_ptr = R.values();
217  const ordinal_type numRows = A.MyLength();
218  const ordinal_type numCols = A.NumVectors();
219  const ordinal_type lda = A.Stride();
220  const ordinal_type ldq = Q.Stride();
221  const ordinal_type ldr = R.stride();
222 
223  const bool contiguousCacheBlocks = false;
224  tsqr_->factorExplicitRaw(numRows, numCols, A_ptr, lda,
225  Q_ptr, ldq, R_ptr, ldr,
226  contiguousCacheBlocks,
227  forceNonnegativeDiagonal);
228  }
229 
260  int revealRank(MV& Q,
261  dense_matrix_type& R,
262  const magnitude_type& tol) {
263  TEUCHOS_TEST_FOR_EXCEPTION(!Q.ConstantStride(), std::invalid_argument,
264  "TsqrAdaptor::"
265  "revealRank: Input MultiVector Q must have constant stride.");
266  prepareTsqr(Q); // Finish initializing TSQR.
267  // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
268  // sure it is the same communicator as the one we are using in
269  // our dist_tsqr_type implementation.
270  return tsqr_->revealRankRaw(Q.MyLength(), Q.NumVectors(),
271  Q.Values(), Q.Stride(),
272  R.values(), R.stride(), tol, false);
273  }
274 
275  private:
277  Teuchos::RCP<node_tsqr_type> nodeTsqr_;
278 
280  Teuchos::RCP<dist_tsqr_type> distTsqr_;
281 
283  Teuchos::RCP<tsqr_type> tsqr_;
284 
286  mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
287 
289  bool ready_;
290 
309  void
310  prepareTsqr(const MV& mv) {
311  if (!ready_) {
312  prepareDistTsqr(mv);
313  ready_ = true;
314  }
315  }
316 
323  void
324  prepareDistTsqr(const MV& mv) {
325  using Teuchos::RCP;
326  using Teuchos::rcp;
327  using TSQR::Epetra::makeTsqrMessenger;
328  typedef TSQR::MessengerBase<scalar_type> base_mess_type;
329 
330  // If mv falls out of scope, its Epetra_Comm may become invalid.
331  // Thus, we clone the input Epetra_Comm, so that the messenger
332  // owns the object.
333  RCP<const Epetra_Comm> comm = rcp(mv.Comm().Clone());
334  RCP<base_mess_type> messBase = makeTsqrMessenger<scalar_type>(comm);
335  distTsqr_->init(messBase);
336  }
337 };
338 
339 } // namespace Epetra
340 
341 #endif // defined(TPETRA_ENABLE_DEPRECATED_CODE) && defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
342 
343 #endif // EPETRA_TSQRADAPTOR_HPP
Function for wrapping Epetra_Comm in a communicator wrapper that TSQR can use.