Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Diagonal_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_DIAGONAL_DEF_HPP
44#define IFPACK2_DIAGONAL_DEF_HPP
45
46#include "Ifpack2_Diagonal_decl.hpp"
47#include "Tpetra_CrsMatrix.hpp"
48
49namespace Ifpack2 {
50
51template<class MatrixType>
52Diagonal<MatrixType>::Diagonal (const Teuchos::RCP<const row_matrix_type>& A) :
53 matrix_ (A),
54 initializeTime_ (0.0),
55 computeTime_ (0.0),
56 applyTime_ (0.0),
57 numInitialize_ (0),
58 numCompute_ (0),
59 numApply_ (0),
60 isInitialized_ (false),
61 isComputed_ (false)
62{}
63
64template<class MatrixType>
65Diagonal<MatrixType>::Diagonal (const Teuchos::RCP<const crs_matrix_type>& A) :
66 matrix_ (A),
67 initializeTime_ (0.0),
68 computeTime_ (0.0),
69 applyTime_ (0.0),
70 numInitialize_ (0),
71 numCompute_ (0),
72 numApply_ (0),
73 isInitialized_ (false),
74 isComputed_ (false)
75{}
76
77template<class MatrixType>
78Diagonal<MatrixType>::Diagonal (const Teuchos::RCP<const vector_type>& diag) :
79 userInverseDiag_ (diag),
80 inverseDiag_ (diag),
81 initializeTime_ (0.0),
82 computeTime_ (0.0),
83 applyTime_ (0.0),
84 numInitialize_ (0),
85 numCompute_ (0),
86 numApply_ (0),
87 isInitialized_ (false),
88 isComputed_ (false)
89{}
90
91template<class MatrixType>
92Diagonal<MatrixType>::~Diagonal ()
93{}
94
95template<class MatrixType>
96Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
97Diagonal<MatrixType>::getDomainMap () const
98{
99 if (matrix_.is_null ()) {
100 if (userInverseDiag_.is_null ()) {
101 TEUCHOS_TEST_FOR_EXCEPTION(
102 true, std::runtime_error, "Ifpack2::Diagonal::getDomainMap: "
103 "The input matrix A is null, and you did not provide a vector of "
104 "inverse diagonal entries. Please call setMatrix() with a nonnull "
105 "input matrix before calling this method.");
106 } else {
107 return userInverseDiag_->getMap ();
108 }
109 } else {
110 return matrix_->getDomainMap ();
111 }
112}
113
114template<class MatrixType>
115Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
116Diagonal<MatrixType>::getRangeMap () const
117{
118 if (matrix_.is_null ()) {
119 if (userInverseDiag_.is_null ()) {
120 TEUCHOS_TEST_FOR_EXCEPTION(
121 true, std::runtime_error, "Ifpack2::Diagonal::getRangeMap: "
122 "The input matrix A is null, and you did not provide a vector of "
123 "inverse diagonal entries. Please call setMatrix() with a nonnull "
124 "input matrix before calling this method.");
125 } else {
126 return userInverseDiag_->getMap ();
127 }
128 } else {
129 return matrix_->getRangeMap ();
130 }
131}
132
133template<class MatrixType>
134void Diagonal<MatrixType>::
135setParameters (const Teuchos::ParameterList& /*params*/)
136{}
137
138template<class MatrixType>
139void Diagonal<MatrixType>::reset ()
140{
141 inverseDiag_ = Teuchos::null;
142 offsets_ = offsets_type ();
143 isInitialized_ = false;
144 isComputed_ = false;
145}
146
147template<class MatrixType>
148void Diagonal<MatrixType>::
149setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
150{
151 if (A.getRawPtr () != matrix_.getRawPtr ()) { // it's a different matrix
152 reset ();
153 matrix_ = A;
154 }
155}
156
157template<class MatrixType>
158void Diagonal<MatrixType>::initialize ()
159{
160 // Either the matrix to precondition must be nonnull, or the user
161 // must have provided a Vector of diagonal entries.
162 TEUCHOS_TEST_FOR_EXCEPTION(
163 matrix_.is_null () && userInverseDiag_.is_null (), std::runtime_error,
164 "Ifpack2::Diagonal::initialize: The matrix to precondition is null, "
165 "and you did not provide a Tpetra::Vector of diagonal entries. "
166 "Please call setMatrix() with a nonnull input before calling this method.");
167
168 // If the user did provide an input matrix, then that takes
169 // precedence over the vector of inverse diagonal entries, if they
170 // provided one earlier. This is only possible if they created this
171 // Diagonal instance using the constructor that takes a
172 // Tpetra::Vector pointer, and then called setMatrix() with a
173 // nonnull input matrix.
174 if (! matrix_.is_null ()) {
175 // If you call initialize(), it means that you are asserting that
176 // the structure of the input sparse matrix may have changed.
177 // This means we should always recompute the diagonal offsets, if
178 // the input matrix is a Tpetra::CrsMatrix.
179 Teuchos::RCP<const crs_matrix_type> A_crs =
180 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (matrix_);
181
182 if (A_crs.is_null ()) {
183 offsets_ = offsets_type (); // offsets are no longer valid
184 }
185 else {
186 const size_t lclNumRows = A_crs->getLocalNumRows ();
187 if (offsets_.extent (0) < lclNumRows) {
188 offsets_ = offsets_type (); // clear first to save memory
189 offsets_ = offsets_type ("offsets", lclNumRows);
190 }
191 A_crs->getCrsGraph ()->getLocalDiagOffsets (offsets_);
192 }
193 }
194
195 isInitialized_ = true;
196 ++numInitialize_;
197}
198
199template<class MatrixType>
200void Diagonal<MatrixType>::compute ()
201{
202 // Either the matrix to precondition must be nonnull, or the user
203 // must have provided a Vector of diagonal entries.
204 TEUCHOS_TEST_FOR_EXCEPTION(
205 matrix_.is_null () && userInverseDiag_.is_null (), std::runtime_error,
206 "Ifpack2::Diagonal::compute: The matrix to precondition is null, "
207 "and you did not provide a Tpetra::Vector of diagonal entries. "
208 "Please call setMatrix() with a nonnull input before calling this method.");
209
210 if (! isInitialized_) {
211 initialize ();
212 }
213
214 // If the user did provide an input matrix, then that takes
215 // precedence over the vector of inverse diagonal entries, if they
216 // provided one earlier. This is only possible if they created this
217 // Diagonal instance using the constructor that takes a
218 // Tpetra::Vector pointer, and then called setMatrix() with a
219 // nonnull input matrix.
220 if (matrix_.is_null ()) { // accept the user's diagonal
221 inverseDiag_ = userInverseDiag_;
222 }
223 else {
224 Teuchos::RCP<vector_type> tmpVec (new vector_type (matrix_->getRowMap ()));
225 Teuchos::RCP<const crs_matrix_type> A_crs =
226 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (matrix_);
227 if (A_crs.is_null ()) {
228 // Get the diagonal entries from the Tpetra::RowMatrix.
229 matrix_->getLocalDiagCopy (*tmpVec);
230 }
231 else {
232 // Get the diagonal entries from the Tpetra::CrsMatrix using the
233 // precomputed offsets.
234 A_crs->getLocalDiagCopy (*tmpVec, offsets_);
235 }
236 tmpVec->reciprocal (*tmpVec); // invert the diagonal entries
237 inverseDiag_ = tmpVec;
238 }
239
240 isComputed_ = true;
241 ++numCompute_;
242}
243
244template<class MatrixType>
245void Diagonal<MatrixType>::
246apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
247 Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
248 Teuchos::ETransp /*mode*/,
249 scalar_type alpha,
250 scalar_type beta) const
251{
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 ! isComputed (), std::runtime_error, "Ifpack2::Diagonal::apply: You "
254 "must first call compute() before you may call apply(). Once you have "
255 "called compute(), you need not call it again unless the values in the "
256 "matrix have changed, or unless you have called setMatrix().");
257
258 // FIXME (mfh 12 Sep 2014) This assumes that row Map == range Map ==
259 // domain Map. If the preconditioner has a matrix, we should ask
260 // the matrix whether we need to do an Import before and/or an
261 // Export after.
262
263 Y.elementWiseMultiply (alpha, *inverseDiag_, X, beta);
264 ++numApply_;
265}
266
267template <class MatrixType>
268int Diagonal<MatrixType>::getNumInitialize() const {
269 return numInitialize_;
270}
271
272template <class MatrixType>
273int Diagonal<MatrixType>::getNumCompute() const {
274 return numCompute_;
275}
276
277template <class MatrixType>
278int Diagonal<MatrixType>::getNumApply() const {
279 return numApply_;
280}
281
282template <class MatrixType>
283double Diagonal<MatrixType>::getInitializeTime() const {
284 return initializeTime_;
285}
286
287template<class MatrixType>
288double Diagonal<MatrixType>::getComputeTime() const {
289 return computeTime_;
290}
291
292template<class MatrixType>
293double Diagonal<MatrixType>::getApplyTime() const {
294 return applyTime_;
295}
296
297template <class MatrixType>
298std::string Diagonal<MatrixType>::description () const
299{
300 std::ostringstream out;
301
302 // Output is a valid YAML dictionary in flow style. If you don't
303 // like everything on a single line, you should call describe()
304 // instead.
305 out << "\"Ifpack2::Diagonal\": "
306 << "{";
307 if (this->getObjectLabel () != "") {
308 out << "Label: \"" << this->getObjectLabel () << "\", ";
309 }
310 if (matrix_.is_null ()) {
311 out << "Matrix: null";
312 }
313 else {
314 out << "Matrix: not null"
315 << ", Global matrix dimensions: ["
316 << matrix_->getGlobalNumRows () << ", "
317 << matrix_->getGlobalNumCols () << "]";
318 }
319
320 out << "}";
321 return out.str ();
322}
323
324template <class MatrixType>
325void Diagonal<MatrixType>::
326describe (Teuchos::FancyOStream &out,
327 const Teuchos::EVerbosityLevel verbLevel) const
328{
329 using std::endl;
330
331 const Teuchos::EVerbosityLevel vl =
332 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
333 if (vl != Teuchos::VERB_NONE) {
334 Teuchos::OSTab tab0 (out);
335 out << "\"Ifpack2::Diagonal\":";
336 Teuchos::OSTab tab1 (out);
337 out << "Template parameter: "
338 << Teuchos::TypeNameTraits<MatrixType>::name () << endl;
339 if (this->getObjectLabel () != "") {
340 out << "Label: \"" << this->getObjectLabel () << "\", ";
341 }
342 out << "Number of initialize calls: " << numInitialize_ << endl
343 << "Number of compute calls: " << numCompute_ << endl
344 << "Number of apply calls: " << numApply_ << endl;
345 }
346}
347
348} // namespace Ifpack2
349
350#define IFPACK2_DIAGONAL_INSTANT(S,LO,GO,N) \
351 template class Ifpack2::Diagonal< Tpetra::RowMatrix<S, LO, GO, N> >;
352
353#endif
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74