Tpetra parallel linear algebra  Version of the Day
Tpetra_CrsMatrixSolveOp.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_CRSMATRIXSOLVEOP_HPP
43 #define TPETRA_CRSMATRIXSOLVEOP_HPP
44 
49 
50 #include <Tpetra_CrsMatrix.hpp>
51 
52 namespace Tpetra {
53 
88  template <class Scalar,
89  class MatScalar = Scalar,
90  class LocalOrdinal = Details::DefaultTypes::local_ordinal_type,
91  class GlobalOrdinal = Details::DefaultTypes::global_ordinal_type,
94  public Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
95  public:
100 
102 
103 
105  CrsMatrixSolveOp (const Teuchos::RCP<const crs_matrix_type>& A) :
106  matrix_ (A)
107  {}
108 
110  virtual ~CrsMatrixSolveOp () {}
111 
113 
115 
118  void
121  Teuchos::ETransp mode = Teuchos::NO_TRANS,
122  Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
123  Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ()) const
124  {
125  typedef Teuchos::ScalarTraits<Scalar> STOS;
126  const char prefix[] = "Tpetra::CrsMatrixSolveOp::apply: ";
127 
128  TEUCHOS_TEST_FOR_EXCEPTION
129  (! matrix_->isFillComplete (), std::runtime_error,
130  prefix << "Underlying matrix is not fill complete.");
131  TEUCHOS_TEST_FOR_EXCEPTION
132  (! matrix_->isLowerTriangular () && ! matrix_->isUpperTriangular (),
133  std::runtime_error, prefix << "The matrix is neither lower nor upper "
134  "triangular. Remember that in Tpetra, triangular-ness is a local "
135  "(per MPI process) property.");
136  TEUCHOS_TEST_FOR_EXCEPTION
137  (X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
138  prefix << "X and Y must have the same number of columns (vectors). "
139  "X.getNumVectors() = " << X.getNumVectors ()
140  << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
141  TEUCHOS_TEST_FOR_EXCEPTION
142  (alpha != STOS::one () || beta != STOS::zero (), std::logic_error,
143  prefix << "The case alpha != 1 or beta != 0 has not yet been "
144  "implemented. Please speak with the Tpetra developers.");
145  if (mode == Teuchos::NO_TRANS) {
146  applyNonTranspose (X,Y);
147  } else if (mode == Teuchos::TRANS || mode == Teuchos::CONJ_TRANS) {
148  applyTranspose (X, Y, mode);
149  } else {
150  TEUCHOS_TEST_FOR_EXCEPTION
151  (true, std::invalid_argument, prefix << "The 'mode' argument has an "
152  "invalid value " << mode << ". Valid values are Teuchos::NO_TRANS="
153  << Teuchos::NO_TRANS << ", Teuchos::TRANS=" << Teuchos::TRANS << ", "
154  "and Teuchos::CONJ_TRANS=" << Teuchos::CONJ_TRANS << ".");
155  }
156  }
157 
159  bool hasTransposeApply () const {
160  return true;
161  }
162 
165  Teuchos::RCP<const map_type> getDomainMap () const
166  {
167  return matrix_->getRangeMap ();
168  }
169 
172  Teuchos::RCP<const map_type> getRangeMap () const {
173  return matrix_->getDomainMap ();
174  }
175 
177  protected:
179 
181  const Teuchos::RCP<const crs_matrix_type> matrix_;
182 
184  mutable Teuchos::RCP<MV> importMV_;
186  mutable Teuchos::RCP<MV> exportMV_;
187 
189  void applyNonTranspose (const MV& X_in, MV& Y_in) const
190  {
191  using Teuchos::NO_TRANS;
192  using Teuchos::null;
193  typedef Teuchos::ScalarTraits<Scalar> ST;
194 
195  // Solve U X = Y or L X = Y
196  // X belongs to domain map, while Y belongs to range map
197 
198  const size_t numVectors = X_in.getNumVectors();
199  Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer =
200  matrix_->getGraph ()->getImporter ();
201  Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
202  matrix_->getGraph ()->getExporter ();
203  Teuchos::RCP<const MV> X;
204 
205  // it is okay if X and Y reference the same data, because we can
206  // perform a triangular solve in-situ. however, we require that
207  // column access to each is strided.
208 
209  // set up import/export temporary multivectors
210  if (importer != null) {
211  if (importMV_ != null && importMV_->getNumVectors () != numVectors) {
212  importMV_ = null;
213  }
214  if (importMV_ == null) {
215  importMV_ = Teuchos::rcp (new MV (matrix_->getColMap (), numVectors));
216  }
217  }
218  if (exporter != null) {
219  if (exportMV_ != null && exportMV_->getNumVectors () != numVectors) {
220  exportMV_ = null;
221  }
222  if (exportMV_ == null) {
223  exportMV_ = Teuchos::rcp (new MV (matrix_->getRowMap (), numVectors));
224  }
225  }
226 
227  // solve(NO_TRANS): RangeMap -> DomainMap
228  // lclMatSolve_: RowMap -> ColMap
229  // importer: DomainMap -> ColMap
230  // exporter: RowMap -> RangeMap
231  //
232  // solve = reverse(exporter) o lclMatSolve_ o reverse(importer)
233  // RangeMap -> RowMap -> ColMap -> DomainMap
234  //
235  // If we have a non-trivial exporter, we must import elements that
236  // are permuted or are on other processors
237  if (exporter != null) {
238  exportMV_->doImport (X_in, *exporter, INSERT);
239  X = exportMV_;
240  }
241  else if (! X_in.isConstantStride ()) {
242  // cannot handle non-constant stride right now
243  // generate a copy of X_in
244  X = Teuchos::rcp (new MV (X_in));
245  }
246  else {
247  // just temporary, so this non-owning RCP is okay
248  X = Teuchos::rcpFromRef (X_in);
249  }
250 
251  // If we have a non-trivial importer, we must export elements that
252  // are permuted or belong to other processes. We will compute
253  // solution into the to-be-exported MV.
254  if (importer != null) {
255  matrix_->template localSolve<Scalar, Scalar> (*X, *importMV_, NO_TRANS);
256  // Make sure target is zero: necessary because we are adding.
257  Y_in.putScalar (ST::zero ());
258  Y_in.doExport (*importMV_, *importer, ADD);
259  }
260  // otherwise, solve into Y
261  else {
262  // can't solve into non-strided multivector
263  if (! Y_in.isConstantStride ()) {
264  // generate a strided copy of Y
265  MV Y (Y_in);
266  matrix_->template localSolve<Scalar, Scalar> (*X, Y, NO_TRANS);
267  Tpetra::deep_copy (Y_in, Y);
268  }
269  else {
270  matrix_->template localSolve<Scalar, Scalar> (*X, Y_in, NO_TRANS);
271  }
272  }
273  }
274 
276  void applyTranspose (const MV& X_in, MV& Y_in, const Teuchos::ETransp mode) const
277  {
278  typedef Teuchos::ScalarTraits<Scalar> ST;
279  using Teuchos::null;
280 
281  TEUCHOS_TEST_FOR_EXCEPTION
282  (mode != Teuchos::TRANS && mode != Teuchos::CONJ_TRANS, std::logic_error,
283  "Tpetra::CrsMatrixSolveOp::applyTranspose: mode is neither TRANS nor "
284  "CONJ_TRANS. Should never get here! Please report this bug to the "
285  "Tpetra developers.");
286 
287  const size_t numVectors = X_in.getNumVectors();
288  Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer =
289  matrix_->getGraph ()->getImporter ();
290  Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter =
291  matrix_->getGraph ()->getExporter ();
292  Teuchos::RCP<const MV> X;
293 
294  // it is okay if X and Y reference the same data, because we can
295  // perform a triangular solve in-situ. however, we require that
296  // column access to each is strided.
297 
298  // set up import/export temporary multivectors
299  if (importer != null) {
300  if (importMV_ != null && importMV_->getNumVectors() != numVectors) {
301  importMV_ = null;
302  }
303  if (importMV_ == null) {
304  importMV_ = Teuchos::rcp( new MV(matrix_->getColMap(),numVectors) );
305  }
306  }
307  if (exporter != null) {
308  if (exportMV_ != null && exportMV_->getNumVectors() != numVectors) {
309  exportMV_ = null;
310  }
311  if (exportMV_ == null) {
312  exportMV_ = Teuchos::rcp( new MV(matrix_->getRowMap(),numVectors) );
313  }
314  }
315 
316  // solve(TRANS): DomainMap -> RangeMap
317  // lclMatSolve_(TRANS): ColMap -> RowMap
318  // importer: DomainMap -> ColMap
319  // exporter: RowMap -> RangeMap
320  //
321  // solve = importer o lclMatSolve_ o exporter
322  // Domainmap -> ColMap -> RowMap -> RangeMap
323  //
324  // If we have a non-trivial importer, we must import elements that
325  // are permuted or are on other processes.
326  if (importer != null) {
327  importMV_->doImport(X_in,*importer,INSERT);
328  X = importMV_;
329  }
330  else if (X_in.isConstantStride() == false) {
331  // cannot handle non-constant stride right now
332  // generate a copy of X_in
333  X = Teuchos::rcp(new MV(X_in));
334  }
335  else {
336  // just temporary, so this non-owning RCP is okay
337  X = Teuchos::rcpFromRef (X_in);
338  }
339 
340 
341  // If we have a non-trivial exporter, we must export elements that
342  // are permuted or belong to other processes. We will compute
343  // solution into the to-be-exported MV; get a view.
344  if (exporter != null) {
345  matrix_->template localSolve<Scalar, Scalar> (*X, *exportMV_,
346  Teuchos::CONJ_TRANS);
347  // Make sure target is zero: necessary because we are adding
348  Y_in.putScalar(ST::zero());
349  Y_in.doExport(*importMV_, *importer, ADD);
350  }
351  // otherwise, solve into Y
352  else {
353  if (Y_in.isConstantStride() == false) {
354  // generate a strided copy of Y
355  MV Y(Y_in);
356  matrix_->template localSolve<Scalar, Scalar> (*X, Y, Teuchos::CONJ_TRANS);
357  Y_in = Y;
358  }
359  else {
360  matrix_->template localSolve<Scalar, Scalar> (*X, Y_in, Teuchos::CONJ_TRANS);
361  }
362  }
363  }
364  };
365 
373  template<class OpScalar,
374  class MatScalar,
375  class LocalOrdinal,
376  class GlobalOrdinal,
377  class Node>
378  Teuchos::RCP<CrsMatrixSolveOp<OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node> >
380  {
382  }
383 
384 } // namespace Tpetra
385 
386 #endif // TPETRA_CRSMATRIXSOLVEOP_HPP
Namespace Tpetra contains the class and methods constituting the Tpetra library.
Sparse matrix that presents a row-oriented interface that lets users read or modify entries...
KokkosClassic::DefaultNode::DefaultNodeType node_type
Default value of Node template parameter.
const Teuchos::RCP< const crs_matrix_type > matrix_
The underlying CrsMatrix.
One or more distributed dense vectors.
void deep_copy(MultiVector< DS, DL, DG, DN, dstClassic > &dst, const MultiVector< SS, SL, SG, SN, srcClassic > &src)
Copy the contents of the MultiVector src into dst.
CrsMatrix< MatScalar, LocalOrdinal, GlobalOrdinal, Node > crs_matrix_type
The specialization of CrsMatrix which this class wraps.
bool isConstantStride() const
Whether this multivector has constant stride between columns.
int local_ordinal_type
Default value of LocalOrdinal template parameter.
Teuchos::RCP< CrsMatrixSolveOp< OpScalar, MatScalar, LocalOrdinal, GlobalOrdinal, Node > > createCrsMatrixSolveOp(const Teuchos::RCP< const CrsMatrix< MatScalar, LocalOrdinal, GlobalOrdinal, Node > > &A)
Nonmember function that wraps a CrsMatrix in a CrsMatrixSolveOp.
Map< LocalOrdinal, GlobalOrdinal, Node > map_type
The specialization of Map which this class uses.
Teuchos::RCP< const map_type > getRangeMap() const
The range Map of this operator. This is the domain Map of the underlying CrsMatrix.
Insert new values that don&#39;t currently exist.
void applyNonTranspose(const MV &X_in, MV &Y_in) const
Do the non-transpose solve.
Abstract interface for operators (e.g., matrices and preconditioners).
Sum new values into existing values.
Teuchos::RCP< const map_type > getDomainMap() const
The domain Map of this operator. This is the range map of the underlying CrsMatrix.
Teuchos::RCP< MV > importMV_
Cached temporary destination of Import operation in apply().
void putScalar(const Scalar &value)
Set all values in the multivector with the given value.
bool hasTransposeApply() const
Whether apply() can solve with the (conjugate) transpose of the matrix.
Describes a parallel distribution of objects over processes.
Wrap a CrsMatrix instance&#39;s triangular solve in an Operator.
void applyTranspose(const MV &X_in, MV &Y_in, const Teuchos::ETransp mode) const
Do the transpose or conjugate transpose solve.
size_t getNumVectors() const
Number of columns in the multivector.
virtual ~CrsMatrixSolveOp()
Destructor (virtual for memory safety of derived classes).
void doExport(const SrcDistObject &source, const Export< LocalOrdinal, GlobalOrdinal, Node > &exporter, CombineMode CM)
Export data into this object using an Export object ("forward mode").
CrsMatrixSolveOp(const Teuchos::RCP< const crs_matrix_type > &A)
Constructor; takes a CrsMatrix to use for local triangular solves.
void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
Compute , where represents the result of the local triangular solve.
Teuchos::RCP< MV > exportMV_
Cached temporary source of Export operation in apply().