Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_IdentitySolver_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
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 
43 #ifndef IFPACK2_IDENTITY_SOLVER_DEF_HPP
44 #define IFPACK2_IDENTITY_SOLVER_DEF_HPP
45 
46 #include "Ifpack2_IdentitySolver_decl.hpp"
47 #include "Tpetra_Map.hpp"
48 #include "Tpetra_MultiVector.hpp"
49 #include "Tpetra_Export.hpp"
50 
51 namespace Ifpack2 {
52 
53 template<class MatrixType>
56  : matrix_ (A),
57  isInitialized_ (false),
58  isComputed_ (false),
59  numInitialize_ (0),
60  numCompute_ (0),
61  numApply_ (0),
62  initializeTime_(0.0),
63  computeTime_(0.0),
64  applyTime_(0.0)
65 {
66 }
67 
68 template<class MatrixType>
70 {
71 }
72 
73 template<class MatrixType>
75 {
76 }
77 
78 template<class MatrixType>
80 {
82  matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver: "
83  "You must call setMatrix() with a nonnull input matrix "
84  "before you may call initialize() or compute().");
85 
87  ! matrix_->getDomainMap ()->isCompatible (* (matrix_->getRangeMap ())),
88  std::invalid_argument,
89  "Ifpack2::IdentitySolver: The domain and range Maps "
90  "of the input matrix must be compatible.");
91 
92  // If the domain and range Maps are not the same, then we need to
93  // construct an Export from the domain Map to the range Map, so that
94  // this operator is really the identity and not a permutation.
95  if (! matrix_->getDomainMap ()->isSameAs (* (matrix_->getRangeMap ()))) {
96  export_ = Teuchos::rcp (new export_type (matrix_->getDomainMap (),
97  matrix_->getRangeMap ()));
98  }
99  else {
100  // If the Export is null, we won't do the Export in apply().
101  // Thus, we need to set it to null here as a flag.
102  export_ = Teuchos::null;
103  }
104 
105  isInitialized_ = true;
106  ++numInitialize_;
107 }
108 
109 template<class MatrixType>
111 {
113  matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver: "
114  "You must call setMatrix() with a nonnull input matrix "
115  "before you may call initialize() or compute().");
116 
117  if (! isInitialized_) {
118  initialize ();
119  }
120 
121  isComputed_ = true;
122  ++numCompute_;
123 }
124 
125 template<class MatrixType>
127 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
128  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
129  Teuchos::ETransp /*mode*/,
130  scalar_type alpha,
131  scalar_type beta) const
132 {
133  using Teuchos::RCP;
135  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
137 
139  ! isComputed (), std::runtime_error,
140  "Ifpack2::IdentitySolver::apply: If compute() has not yet been called, "
141  "or if you have changed the matrix via setMatrix(), "
142  "you must call compute() before you may call this method.");
143 
144  // "Identity solver" does what it says: it's the identity operator.
145  // We have to Export if the domain and range Maps are not the same.
146  // Otherwise, this operator would be a permutation, not the identity.
147  if (export_.is_null ()) {
148  Y.update (alpha, X, beta);
149  }
150  else {
151  if (alpha == STS::one () && beta == STS::zero ()) { // the common case
152  Y.doExport (X, *export_, Tpetra::REPLACE);
153  }
154  else {
155  // We know that the domain and range Maps are compatible. First
156  // bring X into the range Map via Export. Then compute in place
157  // in Y.
158  MV X_tmp (Y.getMap (), Y.getNumVectors ());
159  X_tmp.doExport (X, *export_, Tpetra::REPLACE);
160  Y.update (alpha, X_tmp, beta);
161  }
162  }
163  ++numApply_;
164 }
165 
166 template <class MatrixType>
168  return(numInitialize_);
169 }
170 
171 template <class MatrixType>
173  return(numCompute_);
174 }
175 
176 template <class MatrixType>
178  return(numApply_);
179 }
180 
181 template <class MatrixType>
183  return(initializeTime_);
184 }
185 
186 template<class MatrixType>
188  return(computeTime_);
189 }
190 
191 template<class MatrixType>
193  return(applyTime_);
194 }
195 
196 template <class MatrixType>
198 {
199  std::ostringstream os;
200 
201  // Output is a valid YAML dictionary in flow style. If you don't
202  // like everything on a single line, you should call describe()
203  // instead.
204  os << "\"Ifpack2::IdentitySolver\": {";
205  if (this->getObjectLabel () != "") {
206  os << "Label: \"" << this->getObjectLabel () << "\", ";
207  }
208  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
209  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
210 
211  if (matrix_.is_null ()) {
212  os << "Matrix: null";
213  }
214  else {
215  os << "Matrix: not null"
216  << ", Global matrix dimensions: ["
217  << matrix_->getGlobalNumRows () << ", "
218  << matrix_->getGlobalNumCols () << "]";
219  }
220 
221  os << "}";
222  return os.str ();
223 }
224 
225 template <class MatrixType>
228  const Teuchos::EVerbosityLevel verbLevel) const
229 {
230  using std::endl;
231  const Teuchos::EVerbosityLevel vl
232  = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
233 
234  if (vl != Teuchos::VERB_NONE) {
235  // By convention, describe() should always begin with a tab.
236  Teuchos::OSTab tab0 (out);
237  out << "\"Ifpack2::IdentitySolver\":" << endl;
238  Teuchos::OSTab tab1 (out);
239  out << "MatrixType: " << Teuchos::TypeNameTraits<MatrixType>::name () << endl;
240  out << "numInitialize: " << numInitialize_ << endl;
241  out << "numCompute: " << numCompute_ << endl;
242  out << "numApply: " << numApply_ << endl;
243  }
244 }
245 
246 template <class MatrixType>
248 {
250  matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver::getDomainMap: "
251  "The matrix is null. Please call setMatrix() with a nonnull input "
252  "before calling this method.");
253  return matrix_->getDomainMap ();
254 }
255 
256 template <class MatrixType>
258 {
260  matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver::getRangeMap: "
261  "The matrix is null. Please call setMatrix() with a nonnull input "
262  "before calling this method.");
263  return matrix_->getRangeMap ();
264 }
265 
266 template<class MatrixType>
269 {
270  // Check in serial or one-process mode if the matrix is square.
272  ! A.is_null () && A->getComm ()->getSize () == 1 &&
273  A->getNodeNumRows () != A->getNodeNumCols (),
274  std::runtime_error, "Ifpack2::IdentitySolver::setMatrix: If A's communicator only "
275  "contains one process, then A must be square. Instead, you provided a "
276  "matrix A with " << A->getNodeNumRows () << " rows and "
277  << A->getNodeNumCols () << " columns.");
278 
279  // It's legal for A to be null; in that case, you may not call
280  // initialize() until calling setMatrix() with a nonnull input.
281  // Regardless, setting the matrix invalidates the preconditioner.
282  isInitialized_ = false;
283  isComputed_ = false;
284  export_ = Teuchos::null;
285 
286  matrix_ = A;
287 }
288 
289 } // namespace Ifpack2
290 
291 #define IFPACK2_IDENTITYSOLVER_INSTANT(S,LO,GO,N) \
292  template class Ifpack2::IdentitySolver< Tpetra::RowMatrix<S, LO, GO, N> >;
293 
294 #endif // IFPACK2_IDENTITY_SOLVER_DEF_HPP
int getNumApply() const
Return the number of calls to apply().
Definition: Ifpack2_IdentitySolver_def.hpp:177
MatrixType::node_type node_type
Node type of the input matrix.
Definition: Ifpack2_IdentitySolver_decl.hpp:78
double getComputeTime() const
Return the time spent in compute().
Definition: Ifpack2_IdentitySolver_def.hpp:187
int getNumCompute() const
Return the number of calls to compute().
Definition: Ifpack2_IdentitySolver_def.hpp:172
virtual ~IdentitySolver()
Destructor.
Definition: Ifpack2_IdentitySolver_def.hpp:69
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, and put the result in Y.
Definition: Ifpack2_IdentitySolver_def.hpp:127
Teuchos::RCP< const map_type > getRangeMap() const
Return the Tpetra::Map object associated with the range of this operator.
Definition: Ifpack2_IdentitySolver_def.hpp:257
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setParameters(const Teuchos::ParameterList &params)
Set this object&#39;s parameters.
Definition: Ifpack2_IdentitySolver_def.hpp:74
void initialize()
Initialize.
Definition: Ifpack2_IdentitySolver_def.hpp:79
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const map_type > getDomainMap() const
Return the Tpetra::Map object associated with the domain of this operator.
Definition: Ifpack2_IdentitySolver_def.hpp:247
IdentitySolver(const Teuchos::RCP< const row_matrix_type > &A)
Constructor: Takes the matrix to precondition.
Definition: Ifpack2_IdentitySolver_def.hpp:55
MatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the input matrix.
Definition: Ifpack2_IdentitySolver_decl.hpp:74
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set this preconditioner&#39;s matrix.
Definition: Ifpack2_IdentitySolver_def.hpp:268
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_IdentitySolver_def.hpp:197
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_IdentitySolver_def.hpp:227
void compute()
Compute the preconditioner.
Definition: Ifpack2_IdentitySolver_def.hpp:110
int getNumInitialize() const
Return the number of calls to initialize().
Definition: Ifpack2_IdentitySolver_def.hpp:167
MatrixType::global_ordinal_type global_ordinal_type
Type of the global indices of the input matrix.
Definition: Ifpack2_IdentitySolver_decl.hpp:76
double getInitializeTime() const
Return the time spent in initialize().
Definition: Ifpack2_IdentitySolver_def.hpp:182
MatrixType::scalar_type scalar_type
Type of the entries of the input matrix.
Definition: Ifpack2_IdentitySolver_decl.hpp:72
double getApplyTime() const
Return the time spent in apply().
Definition: Ifpack2_IdentitySolver_def.hpp:192
static std::string name()
bool is_null() const