Ifpack2 Templated Preconditioning Package Version 1.0
Ifpack2_Diagonal_def.hpp
00001 /*@HEADER
00002 // ***********************************************************************
00003 //
00004 //       Ifpack2: Tempated Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2009) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 */
00042 
00043 #ifndef IFPACK2_DIAGONAL_DEF_HPP
00044 #define IFPACK2_DIAGONAL_DEF_HPP
00045 
00046 #include "Ifpack2_Diagonal_decl.hpp"
00047 #include "Tpetra_CrsMatrix_def.hpp"
00048 #include "Ifpack2_Condest.hpp"
00049 
00050 namespace Ifpack2 {
00051 
00052 template<class MatrixType>
00053 Diagonal<MatrixType>::Diagonal (const Teuchos::RCP<const row_matrix_type>& A) :
00054   matrix_ (A),
00055   initializeTime_ (0.0),
00056   computeTime_ (0.0),
00057   applyTime_ (0.0),
00058   numInitialize_ (0),
00059   numCompute_ (0),
00060   numApply_ (0),
00061   condEst_ (-Teuchos::ScalarTraits<magnitude_type>::one ()),
00062   isInitialized_ (false),
00063   isComputed_ (false)
00064 {}
00065 
00066 template<class MatrixType>
00067 Diagonal<MatrixType>::Diagonal (const Teuchos::RCP<const crs_matrix_type>& A) :
00068   matrix_ (A),
00069   initializeTime_ (0.0),
00070   computeTime_ (0.0),
00071   applyTime_ (0.0),
00072   numInitialize_ (0),
00073   numCompute_ (0),
00074   numApply_ (0),
00075   condEst_ (-Teuchos::ScalarTraits<magnitude_type>::one ()),
00076   isInitialized_ (false),
00077   isComputed_ (false)
00078 {}
00079 
00080 template<class MatrixType>
00081 Diagonal<MatrixType>::
00082 Diagonal (const Teuchos::RCP<const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> >& diag) :
00083   userInverseDiag_ (diag),
00084   inverseDiag_ (diag),
00085   initializeTime_ (0.0),
00086   computeTime_ (0.0),
00087   applyTime_ (0.0),
00088   numInitialize_ (0),
00089   numCompute_ (0),
00090   numApply_ (0),
00091   condEst_ (-Teuchos::ScalarTraits<magnitude_type>::one ()),
00092   isInitialized_ (false),
00093   isComputed_ (false)
00094 {}
00095 
00096 
00097 template<class MatrixType>
00098 Diagonal<MatrixType>::~Diagonal()
00099 {}
00100 
00101 template<class MatrixType>
00102 Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
00103 Diagonal<MatrixType>::getDomainMap () const
00104 {
00105   if (matrix_.is_null ()) {
00106     if (userInverseDiag_.is_null ()) {
00107       TEUCHOS_TEST_FOR_EXCEPTION(
00108         true, std::runtime_error, "Ifpack2::Diagonal::getDomainMap: "
00109         "The input matrix A is null, and you did not provide a vector of "
00110         "inverse diagonal entries.  Please call setMatrix() with a nonnull "
00111         "input matrix before calling this method.");
00112     } else {
00113       return userInverseDiag_->getMap ();
00114     }
00115   } else {
00116     return matrix_->getDomainMap ();
00117   }
00118 }
00119 
00120 
00121 template<class MatrixType>
00122 Teuchos::RCP<const typename Diagonal<MatrixType>::map_type>
00123 Diagonal<MatrixType>::getRangeMap () const
00124 {
00125   if (matrix_.is_null ()) {
00126     if (userInverseDiag_.is_null ()) {
00127       TEUCHOS_TEST_FOR_EXCEPTION(
00128         true, std::runtime_error, "Ifpack2::Diagonal::getRangeMap: "
00129         "The input matrix A is null, and you did not provide a vector of "
00130         "inverse diagonal entries.  Please call setMatrix() with a nonnull "
00131         "input matrix before calling this method.");
00132     } else {
00133       return userInverseDiag_->getMap ();
00134     }
00135   } else {
00136     return matrix_->getRangeMap ();
00137   }
00138 }
00139 
00140 
00141 template<class MatrixType>
00142 void Diagonal<MatrixType>::setParameters (const Teuchos::ParameterList& /*params*/)
00143 {}
00144 
00145 
00146 template<class MatrixType>
00147 void Diagonal<MatrixType>::reset ()
00148 {
00149   inverseDiag_ = Teuchos::null;
00150   offsets_ = Teuchos::null;
00151   condEst_ = -Teuchos::ScalarTraits<magnitude_type>::one ();
00152   isInitialized_ = false;
00153   isComputed_ = false;
00154 }
00155 
00156 
00157 template<class MatrixType>
00158 void Diagonal<MatrixType>::
00159 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
00160 {
00161   if (A.getRawPtr () != matrix_.getRawPtr ()) { // it's a different matrix
00162     reset ();
00163     matrix_ = A;
00164   }
00165 }
00166 
00167 
00168 template<class MatrixType>
00169 void Diagonal<MatrixType>::initialize ()
00170 {
00171   // Either the matrix to precondition must be nonnull, or the user
00172   // must have provided a Vector of diagonal entries.
00173   TEUCHOS_TEST_FOR_EXCEPTION(
00174     matrix_.is_null () && userInverseDiag_.is_null (), std::runtime_error,
00175     "Ifpack2::Diagonal::initialize: The matrix to precondition is null, "
00176     "and you did not provide a Tpetra::Vector of diagonal entries.  "
00177     "Please call setMatrix() with a nonnull input before calling this method.");
00178 
00179   // If the user did provide an input matrix, then that takes
00180   // precedence over the vector of inverse diagonal entries, if they
00181   // provided one earlier.  This is only possible if they created this
00182   // Diagonal instance using the constructor that takes a
00183   // Tpetra::Vector pointer, and then called setMatrix() with a
00184   // nonnull input matrix.
00185   if (! matrix_.is_null ()) {
00186     // If you call initialize(), it means that you are asserting that
00187     // the structure of the input sparse matrix may have changed.
00188     // This means we should always recompute the diagonal offsets, if
00189     // the input matrix is a Tpetra::CrsMatrix.
00190     Teuchos::RCP<const crs_matrix_type> A_crs =
00191       Teuchos::rcp_dynamic_cast<const crs_matrix_type> (matrix_);
00192     if (A_crs.is_null ()) {
00193       offsets_ = Teuchos::null; // offsets are no longer valid
00194     }
00195     else {
00196       A_crs->getLocalDiagOffsets (offsets_);
00197     }
00198   }
00199 
00200   isInitialized_ = true;
00201   ++numInitialize_;
00202 }
00203 
00204 template<class MatrixType>
00205 void Diagonal<MatrixType>::compute ()
00206 {
00207   // Either the matrix to precondition must be nonnull, or the user
00208   // must have provided a Vector of diagonal entries.
00209   TEUCHOS_TEST_FOR_EXCEPTION(
00210     matrix_.is_null () && userInverseDiag_.is_null (), std::runtime_error,
00211     "Ifpack2::Diagonal::compute: The matrix to precondition is null, "
00212     "and you did not provide a Tpetra::Vector of diagonal entries.  "
00213     "Please call setMatrix() with a nonnull input before calling this method.");
00214 
00215   if (! isInitialized_) {
00216     initialize ();
00217   }
00218 
00219   // If the user did provide an input matrix, then that takes
00220   // precedence over the vector of inverse diagonal entries, if they
00221   // provided one earlier.  This is only possible if they created this
00222   // Diagonal instance using the constructor that takes a
00223   // Tpetra::Vector pointer, and then called setMatrix() with a
00224   // nonnull input matrix.
00225   if (matrix_.is_null ()) { // accept the user's diagonal
00226     inverseDiag_ = userInverseDiag_;
00227   }
00228   else {
00229     Teuchos::RCP<vector_type> tmpVec (new vector_type (matrix_->getRowMap ()));
00230     Teuchos::RCP<const crs_matrix_type> A_crs =
00231       Teuchos::rcp_dynamic_cast<const crs_matrix_type> (matrix_);
00232     if (A_crs.is_null ()) {
00233       // Get the diagonal entries from the Tpetra::RowMatrix.
00234       matrix_->getLocalDiagCopy (*tmpVec);
00235     }
00236     else {
00237       // Get the diagonal entries from the Tpetra::CrsMatrix using the
00238       // precomputed offsets.
00239       A_crs->getLocalDiagCopy (*tmpVec, offsets_ ());
00240     }
00241     tmpVec->reciprocal (*tmpVec); // invert the diagonal entries
00242     inverseDiag_ = tmpVec;
00243   }
00244 
00245   isComputed_ = true;
00246   ++numCompute_;
00247 }
00248 
00249 template<class MatrixType>
00250 void Diagonal<MatrixType>::
00251 apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
00252        Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
00253        Teuchos::ETransp /*mode*/,
00254        scalar_type alpha,
00255        scalar_type beta) const
00256 {
00257   TEUCHOS_TEST_FOR_EXCEPTION(
00258     ! isComputed (), std::runtime_error, "Ifpack2::Diagonal::apply: You "
00259     "must first call compute() before you may call apply().  Once you have "
00260     "called compute(), you need not call it again unless the values in the "
00261     "matrix have changed, or unless you have called setMatrix().");
00262 
00263   Y.elementWiseMultiply (alpha, *inverseDiag_, X, beta);
00264   ++numApply_;
00265 }
00266 
00267 template<class MatrixType>
00268 typename Teuchos::ScalarTraits<typename MatrixType::scalar_type>::magnitudeType
00269 Diagonal<MatrixType>::
00270 computeCondEst (CondestType CT,
00271                 local_ordinal_type MaxIters,
00272                 magnitude_type Tol,
00273                 const Teuchos::Ptr<const Tpetra::RowMatrix<scalar_type,local_ordinal_type,global_ordinal_type,node_type> > &matrix)
00274 {
00275   const magnitude_type minusOne = -Teuchos::ScalarTraits<magnitude_type>::one ();
00276 
00277   if (! isComputed ()) { // cannot compute right now
00278     return minusOne;
00279   }
00280   // NOTE: this is computing the *local* condest
00281   if (condEst_ == minusOne) {
00282     condEst_ = Ifpack2::Condest (*this, CT, MaxIters, Tol, matrix);
00283   }
00284   return condEst_;
00285 }
00286 
00287 template <class MatrixType>
00288 int Diagonal<MatrixType>::getNumInitialize() const {
00289   return(numInitialize_);
00290 }
00291 
00292 template <class MatrixType>
00293 int Diagonal<MatrixType>::getNumCompute() const {
00294   return(numCompute_);
00295 }
00296 
00297 template <class MatrixType>
00298 int Diagonal<MatrixType>::getNumApply() const {
00299   return(numApply_);
00300 }
00301 
00302 template <class MatrixType>
00303 double Diagonal<MatrixType>::getInitializeTime() const {
00304   return(initializeTime_);
00305 }
00306 
00307 template<class MatrixType>
00308 double Diagonal<MatrixType>::getComputeTime() const {
00309   return(computeTime_);
00310 }
00311 
00312 template<class MatrixType>
00313 double Diagonal<MatrixType>::getApplyTime() const {
00314   return(applyTime_);
00315 }
00316 
00317 template <class MatrixType>
00318 std::string Diagonal<MatrixType>::description () const
00319 {
00320   std::ostringstream out;
00321 
00322   // Output is a valid YAML dictionary in flow style.  If you don't
00323   // like everything on a single line, you should call describe()
00324   // instead.
00325   out << "\"Ifpack2::Diagonal\": "
00326       << "{";
00327   if (this->getObjectLabel () != "") {
00328     out << "Label: \"" << this->getObjectLabel () << "\", ";
00329   }
00330   if (matrix_.is_null ()) {
00331     out << "Matrix: null";
00332   }
00333   else {
00334     out << "Matrix: not null"
00335         << ", Global matrix dimensions: ["
00336         << matrix_->getGlobalNumRows () << ", "
00337         << matrix_->getGlobalNumCols () << "]";
00338   }
00339 
00340   out << "}";
00341   return out.str ();
00342 }
00343 
00344 template <class MatrixType>
00345 void Diagonal<MatrixType>::
00346 describe (Teuchos::FancyOStream &out,
00347           const Teuchos::EVerbosityLevel verbLevel) const
00348 {
00349   using std::endl;
00350 
00351   const Teuchos::EVerbosityLevel vl =
00352     (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
00353   if (vl != Teuchos::VERB_NONE) {
00354     Teuchos::OSTab tab0 (out);
00355     out << "\"Ifpack2::Diagonal\":";
00356     Teuchos::OSTab tab1 (out);
00357     out << "Template parameter: "
00358         << Teuchos::TypeNameTraits<MatrixType>::name () << endl;
00359     if (this->getObjectLabel () != "") {
00360       out << "Label: \"" << this->getObjectLabel () << "\", ";
00361     }
00362     out << "Number of initialize calls: " << numInitialize_ << endl
00363         << "Number of compute calls: " << numCompute_ << endl
00364         << "Number of apply calls: " << numApply_ << endl;
00365   }
00366 }
00367 
00368 } // namespace Ifpack2
00369 
00370 #define IFPACK2_DIAGONAL_INSTANT(S,LO,GO,N)                            \
00371   template class Ifpack2::Diagonal< Tpetra::CrsMatrix<S, LO, GO, N> >; \
00372   template class Ifpack2::Diagonal< Tpetra::RowMatrix<S, LO, GO, N> >;
00373 
00374 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends