Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Diagonal_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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 
49 namespace Ifpack2 {
50 
51 template<class MatrixType>
52 Diagonal<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 
64 template<class MatrixType>
65 Diagonal<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 
77 template<class MatrixType>
78 Diagonal<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 
91 template<class MatrixType>
92 Diagonal<MatrixType>::~Diagonal ()
93 {}
94 
95 template<class MatrixType>
96 Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
97 Diagonal<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 
114 template<class MatrixType>
115 Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
116 Diagonal<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 
133 template<class MatrixType>
134 void Diagonal<MatrixType>::
135 setParameters (const Teuchos::ParameterList& /*params*/)
136 {}
137 
138 template<class MatrixType>
139 void Diagonal<MatrixType>::reset ()
140 {
141  inverseDiag_ = Teuchos::null;
142  offsets_ = Teuchos::null;
143  isInitialized_ = false;
144  isComputed_ = false;
145 }
146 
147 template<class MatrixType>
148 void Diagonal<MatrixType>::
149 setMatrix (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 
157 template<class MatrixType>
158 void 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  if (A_crs.is_null ()) {
182  offsets_ = Teuchos::null; // offsets are no longer valid
183  }
184  else {
185  A_crs->getLocalDiagOffsets (offsets_);
186  }
187  }
188 
189  isInitialized_ = true;
190  ++numInitialize_;
191 }
192 
193 template<class MatrixType>
194 void Diagonal<MatrixType>::compute ()
195 {
196  // Either the matrix to precondition must be nonnull, or the user
197  // must have provided a Vector of diagonal entries.
198  TEUCHOS_TEST_FOR_EXCEPTION(
199  matrix_.is_null () && userInverseDiag_.is_null (), std::runtime_error,
200  "Ifpack2::Diagonal::compute: The matrix to precondition is null, "
201  "and you did not provide a Tpetra::Vector of diagonal entries. "
202  "Please call setMatrix() with a nonnull input before calling this method.");
203 
204  if (! isInitialized_) {
205  initialize ();
206  }
207 
208  // If the user did provide an input matrix, then that takes
209  // precedence over the vector of inverse diagonal entries, if they
210  // provided one earlier. This is only possible if they created this
211  // Diagonal instance using the constructor that takes a
212  // Tpetra::Vector pointer, and then called setMatrix() with a
213  // nonnull input matrix.
214  if (matrix_.is_null ()) { // accept the user's diagonal
215  inverseDiag_ = userInverseDiag_;
216  }
217  else {
218  Teuchos::RCP<vector_type> tmpVec (new vector_type (matrix_->getRowMap ()));
219  Teuchos::RCP<const crs_matrix_type> A_crs =
220  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (matrix_);
221  if (A_crs.is_null ()) {
222  // Get the diagonal entries from the Tpetra::RowMatrix.
223  matrix_->getLocalDiagCopy (*tmpVec);
224  }
225  else {
226  // Get the diagonal entries from the Tpetra::CrsMatrix using the
227  // precomputed offsets.
228  A_crs->getLocalDiagCopy (*tmpVec, offsets_ ());
229  }
230  tmpVec->reciprocal (*tmpVec); // invert the diagonal entries
231  inverseDiag_ = tmpVec;
232  }
233 
234  isComputed_ = true;
235  ++numCompute_;
236 }
237 
238 template<class MatrixType>
239 void Diagonal<MatrixType>::
240 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
241  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
242  Teuchos::ETransp /*mode*/,
243  scalar_type alpha,
244  scalar_type beta) const
245 {
246  TEUCHOS_TEST_FOR_EXCEPTION(
247  ! isComputed (), std::runtime_error, "Ifpack2::Diagonal::apply: You "
248  "must first call compute() before you may call apply(). Once you have "
249  "called compute(), you need not call it again unless the values in the "
250  "matrix have changed, or unless you have called setMatrix().");
251 
252  // FIXME (mfh 12 Sep 2014) This assumes that row Map == range Map ==
253  // domain Map. If the preconditioner has a matrix, we should ask
254  // the matrix whether we need to do an Import before and/or an
255  // Export after.
256 
257  Y.elementWiseMultiply (alpha, *inverseDiag_, X, beta);
258  ++numApply_;
259 }
260 
261 template <class MatrixType>
262 int Diagonal<MatrixType>::getNumInitialize() const {
263  return numInitialize_;
264 }
265 
266 template <class MatrixType>
267 int Diagonal<MatrixType>::getNumCompute() const {
268  return numCompute_;
269 }
270 
271 template <class MatrixType>
272 int Diagonal<MatrixType>::getNumApply() const {
273  return numApply_;
274 }
275 
276 template <class MatrixType>
277 double Diagonal<MatrixType>::getInitializeTime() const {
278  return initializeTime_;
279 }
280 
281 template<class MatrixType>
282 double Diagonal<MatrixType>::getComputeTime() const {
283  return computeTime_;
284 }
285 
286 template<class MatrixType>
287 double Diagonal<MatrixType>::getApplyTime() const {
288  return applyTime_;
289 }
290 
291 template <class MatrixType>
292 std::string Diagonal<MatrixType>::description () const
293 {
294  std::ostringstream out;
295 
296  // Output is a valid YAML dictionary in flow style. If you don't
297  // like everything on a single line, you should call describe()
298  // instead.
299  out << "\"Ifpack2::Diagonal\": "
300  << "{";
301  if (this->getObjectLabel () != "") {
302  out << "Label: \"" << this->getObjectLabel () << "\", ";
303  }
304  if (matrix_.is_null ()) {
305  out << "Matrix: null";
306  }
307  else {
308  out << "Matrix: not null"
309  << ", Global matrix dimensions: ["
310  << matrix_->getGlobalNumRows () << ", "
311  << matrix_->getGlobalNumCols () << "]";
312  }
313 
314  out << "}";
315  return out.str ();
316 }
317 
318 template <class MatrixType>
319 void Diagonal<MatrixType>::
320 describe (Teuchos::FancyOStream &out,
321  const Teuchos::EVerbosityLevel verbLevel) const
322 {
323  using std::endl;
324 
325  const Teuchos::EVerbosityLevel vl =
326  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
327  if (vl != Teuchos::VERB_NONE) {
328  Teuchos::OSTab tab0 (out);
329  out << "\"Ifpack2::Diagonal\":";
330  Teuchos::OSTab tab1 (out);
331  out << "Template parameter: "
332  << Teuchos::TypeNameTraits<MatrixType>::name () << endl;
333  if (this->getObjectLabel () != "") {
334  out << "Label: \"" << this->getObjectLabel () << "\", ";
335  }
336  out << "Number of initialize calls: " << numInitialize_ << endl
337  << "Number of compute calls: " << numCompute_ << endl
338  << "Number of apply calls: " << numApply_ << endl;
339  }
340 }
341 
342 } // namespace Ifpack2
343 
344 #define IFPACK2_DIAGONAL_INSTANT(S,LO,GO,N) \
345  template class Ifpack2::Diagonal< Tpetra::RowMatrix<S, LO, GO, N> >;
346 
347 #endif
virtual bool isComputed() const
Returns true if partitions have been computed successfully.
Definition: Ifpack2_OverlappingPartitioner_def.hpp:357
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72