Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
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
51namespace Ifpack2 {
52
53template<class MatrixType>
55IdentitySolver (const Teuchos::RCP<const row_matrix_type>& A)
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
68template<class MatrixType>
72
73template<class MatrixType>
74void IdentitySolver<MatrixType>::setParameters (const Teuchos::ParameterList& /*params*/)
75{
76}
77
78template<class MatrixType>
80{
81 TEUCHOS_TEST_FOR_EXCEPTION(
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
86 TEUCHOS_TEST_FOR_EXCEPTION(
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
109template<class MatrixType>
111{
112 TEUCHOS_TEST_FOR_EXCEPTION(
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
125template<class MatrixType>
127apply (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;
134 typedef Teuchos::ScalarTraits<scalar_type> STS;
135 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
137
138 TEUCHOS_TEST_FOR_EXCEPTION(
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
166template <class MatrixType>
168 return(numInitialize_);
169}
170
171template <class MatrixType>
173 return(numCompute_);
174}
175
176template <class MatrixType>
178 return(numApply_);
179}
180
181template <class MatrixType>
183 return(initializeTime_);
184}
185
186template<class MatrixType>
188 return(computeTime_);
189}
190
191template<class MatrixType>
193 return(applyTime_);
194}
195
196template <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
225template <class MatrixType>
227describe (Teuchos::FancyOStream& out,
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
246template <class MatrixType>
247Teuchos::RCP<const typename IdentitySolver<MatrixType>::map_type> IdentitySolver<MatrixType>::getDomainMap() const
248{
249 TEUCHOS_TEST_FOR_EXCEPTION(
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
256template <class MatrixType>
257Teuchos::RCP<const typename IdentitySolver<MatrixType>::map_type> IdentitySolver<MatrixType>::getRangeMap() const
258{
259 TEUCHOS_TEST_FOR_EXCEPTION(
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
266template<class MatrixType>
268setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
269{
270 // Check in serial or one-process mode if the matrix is square.
271 TEUCHOS_TEST_FOR_EXCEPTION(
272 ! A.is_null () && A->getComm ()->getSize () == 1 &&
273 A->getLocalNumRows () != A->getLocalNumCols (),
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->getLocalNumRows () << " rows and "
277 << A->getLocalNumCols () << " 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
void initialize()
Initialize.
Definition Ifpack2_IdentitySolver_def.hpp:79
int getNumInitialize() const
Return the number of calls to initialize().
Definition Ifpack2_IdentitySolver_def.hpp:167
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
MatrixType::node_type node_type
Node type of the input matrix.
Definition Ifpack2_IdentitySolver_decl.hpp:78
void setParameters(const Teuchos::ParameterList &params)
Set this object's parameters.
Definition Ifpack2_IdentitySolver_def.hpp:74
double getInitializeTime() const
Return the time spent in initialize().
Definition Ifpack2_IdentitySolver_def.hpp:182
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
int getNumCompute() const
Return the number of calls to compute().
Definition Ifpack2_IdentitySolver_def.hpp:172
MatrixType::local_ordinal_type local_ordinal_type
Type of the local indices of the input matrix.
Definition Ifpack2_IdentitySolver_decl.hpp:74
IdentitySolver(const Teuchos::RCP< const row_matrix_type > &A)
Constructor: Takes the matrix to precondition.
Definition Ifpack2_IdentitySolver_def.hpp:55
double getApplyTime() const
Return the time spent in apply().
Definition Ifpack2_IdentitySolver_def.hpp:192
std::string description() const
Return a simple one-line description of this object.
Definition Ifpack2_IdentitySolver_def.hpp:197
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set this preconditioner's matrix.
Definition Ifpack2_IdentitySolver_def.hpp:268
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
double getComputeTime() const
Return the time spent in compute().
Definition Ifpack2_IdentitySolver_def.hpp:187
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
MatrixType::global_ordinal_type global_ordinal_type
Type of the global indices of the input matrix.
Definition Ifpack2_IdentitySolver_decl.hpp:76
MatrixType::scalar_type scalar_type
Type of the entries of the input matrix.
Definition Ifpack2_IdentitySolver_decl.hpp:72
int getNumApply() const
Return the number of calls to apply().
Definition Ifpack2_IdentitySolver_def.hpp:177
virtual ~IdentitySolver()
Destructor.
Definition Ifpack2_IdentitySolver_def.hpp:69
void compute()
Compute the preconditioner.
Definition Ifpack2_IdentitySolver_def.hpp:110
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74