Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_TpetraLinearOp_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
4 //
5 // Copyright 2004 NTESS and the Thyra contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef THYRA_TPETRA_LINEAR_OP_HPP
11 #define THYRA_TPETRA_LINEAR_OP_HPP
12 
15 #include "Teuchos_ScalarTraits.hpp"
17 
18 #include "Tpetra_CrsMatrix.hpp"
19 
20 #ifdef HAVE_THYRA_TPETRA_EPETRA
22 #endif
23 
24 namespace Thyra {
25 
26 
27 #ifdef HAVE_THYRA_TPETRA_EPETRA
28 
29 // Utilites
30 
31 
33  template<class Scalar, class LocalOrdinal, class GlobalOrdinal>
34 class GetTpetraEpetraRowMatrixWrapper {
35 public:
36  template<class TpetraMatrixType>
37  static
38  RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
39  get(const RCP<TpetraMatrixType> &tpetraMatrix)
40  {
41  return Teuchos::null;
42  }
43 };
44 
45 
46 // NOTE: We could support other ordinal types, but we have to
47 // specialize the EpetraRowMatrix
48 template<>
49 class GetTpetraEpetraRowMatrixWrapper<double, int, int> {
50 public:
51  template<class TpetraMatrixType>
52  static
53  RCP<Tpetra::EpetraRowMatrix<TpetraMatrixType> >
54  get(const RCP<TpetraMatrixType> &tpetraMatrix)
55  {
56  return Teuchos::rcp(
57  new Tpetra::EpetraRowMatrix<TpetraMatrixType>(tpetraMatrix,
59  *convertTpetraToThyraComm(tpetraMatrix->getRowMap()->getComm())
60  )
61  )
62  );
63  }
64 };
65 
66 
67 #endif // HAVE_THYRA_TPETRA_EPETRA
68 
69 
70 template <class Scalar>
71 inline
74 {
77  Exceptions::OpNotSupported,
78  "For complex scalars such as " + Teuchos::TypeNameTraits<Scalar>::name() +
79  ", Tpetra does not support conjugation without transposition."
80  );
81  return Teuchos::NO_TRANS; // For non-complex scalars, CONJ is equivalent to NOTRANS.
82 }
83 
84 
85 template <class Scalar>
86 inline
88 convertToTeuchosTransMode(const Thyra::EOpTransp transp)
89 {
90  switch (transp) {
91  case NOTRANS: return Teuchos::NO_TRANS;
92  case CONJ: return convertConjNoTransToTeuchosTransMode<Scalar>();
93  case TRANS: return Teuchos::TRANS;
94  case CONJTRANS: return Teuchos::CONJ_TRANS;
95  }
96 
97  // Should not escape the switch
99 }
100 
101 
102 // Constructors/initializers
103 
104 
105 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107 {}
108 
109 
110 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112  const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
113  const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
114  const RCP<Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
115  )
116 {
117  initializeImpl(rangeSpace, domainSpace, tpetraOperator);
118 }
119 
120 
121 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
123  const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
124  const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
125  const RCP<const Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tpetraOperator
126  )
127 {
128  initializeImpl(rangeSpace, domainSpace, tpetraOperator);
129 }
130 
131 
132 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
135 {
136  return tpetraOperator_.getNonconstObj();
137 }
138 
139 
140 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
143 {
144  return tpetraOperator_;
145 }
146 
147 
148 // Public Overridden functions from LinearOpBase
149 
150 
151 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
154 {
155  return rangeSpace_;
156 }
157 
158 
159 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
162 {
163  return domainSpace_;
164 }
165 
166 
167 // Overridden from EpetraLinearOpBase
168 
169 
170 #ifdef HAVE_THYRA_TPETRA_EPETRA
171 
172 
173 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175  const Ptr<RCP<Epetra_Operator> > &epetraOp,
176  const Ptr<EOpTransp> &epetraOpTransp,
177  const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
178  const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
179  )
180 {
182 }
183 
184 
185 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
186 void TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getEpetraOpView(
187  const Ptr<RCP<const Epetra_Operator> > &epetraOp,
188  const Ptr<EOpTransp> &epetraOpTransp,
189  const Ptr<EApplyEpetraOpAs> &epetraOpApplyAs,
190  const Ptr<EAdjointEpetraOp> &epetraOpAdjointSupport
191  ) const
192 {
193  using Teuchos::rcp_dynamic_cast;
194  typedef Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpetraRowMatrix_t;
195  if (nonnull(tpetraOperator_)) {
196  if (is_null(epetraOp_)) {
197  epetraOp_ = GetTpetraEpetraRowMatrixWrapper<Scalar,LocalOrdinal,GlobalOrdinal>::get(
198  rcp_dynamic_cast<const TpetraRowMatrix_t>(tpetraOperator_.getConstObj(), true));
199  }
200  *epetraOp = epetraOp_;
201  *epetraOpTransp = NOTRANS;
202  *epetraOpApplyAs = EPETRA_OP_APPLY_APPLY;
203  *epetraOpAdjointSupport = ( tpetraOperator_->hasTransposeApply()
205  }
206  else {
207  *epetraOp = Teuchos::null;
208  }
209 }
210 
211 
212 #endif // HAVE_THYRA_TPETRA_EPETRA
213 
214 
215 // Protected Overridden functions from LinearOpBase
216 
217 
218 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
220  Thyra::EOpTransp M_trans) const
221 {
222  if (is_null(tpetraOperator_))
223  return false;
224 
225  if (M_trans == NOTRANS)
226  return true;
227 
228  if (M_trans == CONJ) {
229  // For non-complex scalars, CONJ is always supported since it is equivalent to NO_TRANS.
230  // For complex scalars, Tpetra does not support conjugation without transposition.
232  }
233 
234  return tpetraOperator_->hasTransposeApply();
235 }
236 
237 
238 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
240  const Thyra::EOpTransp M_trans,
241  const Thyra::MultiVectorBase<Scalar> &X_in,
242  const Teuchos::Ptr<Thyra::MultiVectorBase<Scalar> > &Y_inout,
243  const Scalar alpha,
244  const Scalar beta
245  ) const
246 {
247  using Teuchos::rcpFromRef;
248  using Teuchos::rcpFromPtr;
250  ConverterT;
251  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
252  TpetraMultiVector_t;
253 
254  // Get Tpetra::MultiVector objects for X and Y
255 
257  ConverterT::getConstTpetraMultiVector(rcpFromRef(X_in));
258 
259  const RCP<TpetraMultiVector_t> tY =
260  ConverterT::getTpetraMultiVector(rcpFromPtr(Y_inout));
261 
262  const Teuchos::ETransp tTransp = convertToTeuchosTransMode<Scalar>(M_trans);
263 
264  // Apply the operator
265 
266  tpetraOperator_->apply(*tX, *tY, tTransp, alpha, beta);
267  // CAG: Commented out since the purpose seems unclear.
268  // Tpetra apply should do all the necessary fencing.
269  // Kokkos::fence();
270 }
271 
272 // Protected member functions overridden from ScaledLinearOpBase
273 
274 
275 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
277 {
278  return true;
279 }
280 
281 
282 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
284 {
285  return true;
286 }
287 
288 
289 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290 void
292 scaleLeftImpl(const VectorBase<Scalar> &row_scaling_in)
293 {
294  using Teuchos::rcpFromRef;
295 
298 
300  Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),true);
301 
302  rowMatrix->leftScale(*row_scaling);
303  Kokkos::fence();
304 }
305 
306 
307 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
308 void
310 scaleRightImpl(const VectorBase<Scalar> &col_scaling_in)
311 {
312  using Teuchos::rcpFromRef;
313 
316 
318  Teuchos::rcp_dynamic_cast<Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getNonconstObj(),true);
319 
320  rowMatrix->rightScale(*col_scaling);
321  Kokkos::fence();
322 }
323 
324 // Protected member functions overridden from RowStatLinearOpBase
325 
326 
327 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
330  const RowStatLinearOpBaseUtils::ERowStat rowStat) const
331 {
332  if (is_null(tpetraOperator_))
333  return false;
334 
335  switch (rowStat) {
336  case RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM:
337  case RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM:
338  return true;
339  default:
340  return false;
341  }
342 
343 }
344 
345 
346 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
348  const RowStatLinearOpBaseUtils::ERowStat rowStat,
349  const Ptr<VectorBase<Scalar> > &rowStatVec_in
350  ) const
351 {
352  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
353  TpetraVector_t;
354  typedef Teuchos::ScalarTraits<Scalar> STS;
355  typedef typename STS::magnitudeType MT;
356  typedef Teuchos::ScalarTraits<MT> STM;
357 
358  if ( (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) ||
359  (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_ROW_SUM) ) {
360 
361  TEUCHOS_ASSERT(nonnull(tpetraOperator_));
362  TEUCHOS_ASSERT(nonnull(rowStatVec_in));
363 
364  // Currently we only support the case of row sums for a concrete
365  // Tpetra::CrsMatrix where (1) the entire row is stored on a
366  // single process and (2) that the domain map, the range map and
367  // the row map are the SAME. These checks enforce that. Later on
368  // we hope to add complete support for any mapping to the concrete
369  // tpetra matrix types.
370 
371  const RCP<TpetraVector_t> tRowSumVec =
373 
375  Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(tpetraOperator_.getConstObj(),true);
376 
377  // EGP: The following assert fails when row sum scaling is applied to blocked Tpetra operators, but without the assert, the correct row sum scaling is obtained.
378  // Furthermore, no valgrind memory errors occur in this case when the assert is removed.
379  //TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getDomainMap()));
380  TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tCrsMatrix->getRangeMap()));
381  TEUCHOS_ASSERT(tCrsMatrix->getRowMap()->isSameAs(*tRowSumVec->getMap()));
382 
383  size_t numMyRows = tCrsMatrix->getLocalNumRows();
384 
385  using crs_t = Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
386  typename crs_t::local_inds_host_view_type indices;
387  typename crs_t::values_host_view_type values;
388 
389 
390  for (size_t row=0; row < numMyRows; ++row) {
391  MT sum = STM::zero ();
392  tCrsMatrix->getLocalRowView (row, indices, values);
393 
394  for (int col = 0; col < (int) values.size(); ++col) {
395  sum += STS::magnitude (values[col]);
396  }
397 
398  if (rowStat == RowStatLinearOpBaseUtils::ROW_STAT_INV_ROW_SUM) {
399  if (sum < STM::sfmin ()) {
400  TEUCHOS_TEST_FOR_EXCEPTION(sum == STM::zero (), std::runtime_error,
401  "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Inverse row sum "
402  << "requested for a matrix where one of the rows has a zero row sum!");
403  sum = STM::one () / STM::sfmin ();
404  }
405  else {
406  sum = STM::one () / sum;
407  }
408  }
409 
410  tRowSumVec->replaceLocalValue (row, Scalar (sum));
411  }
412 
413  }
414  else {
415  TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
416  "Error - Thyra::TpetraLinearOp::getRowStatImpl() - Column sum support not implemented!");
417  }
418  Kokkos::fence();
419 }
420 
421 
422 // private
423 
424 
425 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
426 template<class TpetraOperator_t>
428  const RCP<const VectorSpaceBase<Scalar> > &rangeSpace,
429  const RCP<const VectorSpaceBase<Scalar> > &domainSpace,
430  const RCP<TpetraOperator_t> &tpetraOperator
431  )
432 {
433 #ifdef THYRA_DEBUG
434  TEUCHOS_ASSERT(nonnull(rangeSpace));
435  TEUCHOS_ASSERT(nonnull(domainSpace));
436  TEUCHOS_ASSERT(nonnull(tpetraOperator));
437  // ToDo: Assert that spaces are comparible with tpetraOperator
438 #endif
439  rangeSpace_ = rangeSpace;
440  domainSpace_ = domainSpace;
441  tpetraOperator_ = tpetraOperator;
442 }
443 
444 
445 } // namespace Thyra
446 
447 
448 #endif // THYRA_TPETRA_LINEAR_OP_HPP
bool is_null(const boost::shared_ptr< T > &p)
void initializeImpl(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< TpetraOperator_t > &tpetraOperator)
void constInitialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< const Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraOperator)
Initialize.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
virtual void getRowStatImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat, const Ptr< VectorBase< Scalar > > &rowStatVec) const
bool opSupportedImpl(Thyra::EOpTransp M_trans) const
Traits class that enables the extraction of Tpetra operator/vector objects wrapped in Thyra operator/...
RCP< const Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraOperator() const
Get embedded const Tpetra::Operator.
void applyImpl(const Thyra::EOpTransp M_trans, const Thyra::MultiVectorBase< Scalar > &X_in, const Teuchos::Ptr< Thyra::MultiVectorBase< Scalar > > &Y_inout, const Scalar alpha, const Scalar beta) const
static RCP< const Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getConstTpetraVector(const RCP< const VectorBase< Scalar > > &v)
Get a const Tpetra::Vector from a const Thyra::VectorBase object.
virtual bool rowStatIsSupportedImpl(const RowStatLinearOpBaseUtils::ERowStat rowStat) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual bool supportsScaleLeftImpl() const
RCP< const Thyra::VectorSpaceBase< Scalar > > range() const
RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraOperator()
Get embedded non-const Tpetra::Operator.
static RCP< Tpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getTpetraVector(const RCP< VectorBase< Scalar > > &v)
Get a non-const Tpetra::Vector from a non-const Thyra::VectorBase object.
Teuchos::ETransp convertConjNoTransToTeuchosTransMode()
TpetraLinearOp()
Construct to uninitialized.
Teuchos::ETransp convertToTeuchosTransMode(const Thyra::EOpTransp transp)
RCP< const Epetra_Comm > get_Epetra_Comm(const Teuchos::Comm< Ordinal > &comm)
Get (or create) and Epetra_Comm given a Teuchos::Comm object.
virtual void scaleLeftImpl(const VectorBase< Scalar > &row_scaling)
bool nonnull(const boost::shared_ptr< T > &p)
virtual void scaleRightImpl(const VectorBase< Scalar > &col_scaling)
virtual bool supportsScaleRightImpl() const
void initialize(const RCP< const VectorSpaceBase< Scalar > > &rangeSpace, const RCP< const VectorSpaceBase< Scalar > > &domainSpace, const RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &tpetraOperator)
Initialize.
RCP< const Teuchos::Comm< Ordinal > > convertTpetraToThyraComm(const RCP< const Teuchos::Comm< int > > &tpetraComm)
Given an Tpetra Teuchos::Comm&lt;int&gt; object, return an equivalent Teuchos::Comm&lt;Ordinal&gt; object...
#define TEUCHOS_ASSERT(assertion_test)
RCP< const Thyra::VectorSpaceBase< Scalar > > domain() const
Concrete Thyra::LinearOpBase subclass for Tpetra::Operator.
Apply using Epetra_Operator::Apply(...)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)