Ifpack2_IdentitySolver_def.hpp

/*@HEADER
// ***********************************************************************
//
//       Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2009) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
//@HEADER
*/

#ifndef IFPACK2_IDENTITY_SOLVER_DEF_HPP
#define IFPACK2_IDENTITY_SOLVER_DEF_HPP

#include "Ifpack2_IdentitySolver_decl.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Export.hpp"

namespace Ifpack2 {

template<class MatrixType>
IdentitySolver<MatrixType>::
IdentitySolver (const Teuchos::RCP<const row_matrix_type>& A)
  : matrix_ (A),
    isInitialized_ (false),
    isComputed_ (false),
    numInitialize_ (0),
    numCompute_ (0),
    numApply_ (0),
    initializeTime_(0.0),
    computeTime_(0.0),
    applyTime_(0.0)
{
}

template<class MatrixType>
IdentitySolver<MatrixType>::~IdentitySolver ()
{
}

template<class MatrixType>
void IdentitySolver<MatrixType>::setParameters (const Teuchos::ParameterList& /*params*/)
{
}

template<class MatrixType>
void IdentitySolver<MatrixType>::initialize ()
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver: "
    "You must call setMatrix() with a nonnull input matrix "
    "before you may call initialize() or compute().");

  TEUCHOS_TEST_FOR_EXCEPTION(
    ! matrix_->getDomainMap ()->isCompatible (* (matrix_->getRangeMap ())),
    std::invalid_argument,
    "Ifpack2::IdentitySolver: The domain and range Maps "
    "of the input matrix must be compatible.");

  // If the domain and range Maps are not the same, then we need to
  // construct an Export from the domain Map to the range Map, so that
  // this operator is really the identity and not a permutation.
  if (! matrix_->getDomainMap ()->isSameAs (* (matrix_->getRangeMap ()))) {
    export_ = Teuchos::rcp (new export_type (matrix_->getDomainMap (),
                                             matrix_->getRangeMap ()));
  }
  else {
    // If the Export is null, we won't do the Export in apply().
    // Thus, we need to set it to null here as a flag.
    export_ = Teuchos::null;
  }

  isInitialized_ = true;
  ++numInitialize_;
}

template<class MatrixType>
void IdentitySolver<MatrixType>::compute ()
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver: "
    "You must call setMatrix() with a nonnull input matrix "
    "before you may call initialize() or compute().");

  if (! isInitialized_) {
    initialize ();
  }

  isComputed_ = true;
  ++numCompute_;
}

template<class MatrixType>
void IdentitySolver<MatrixType>::
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*/,
       scalar_type alpha,
       scalar_type beta) const
{
  using Teuchos::RCP;
  typedef Teuchos::ScalarTraits<scalar_type> STS;
  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
                              global_ordinal_type, node_type> MV;

  TEUCHOS_TEST_FOR_EXCEPTION(
    ! isComputed (), std::runtime_error,
    "Ifpack2::IdentitySolver::apply: If compute() has not yet been called, "
    "or if you have changed the matrix via setMatrix(), "
    "you must call compute() before you may call this method.");

  // "Identity solver" does what it says: it's the identity operator.
  // We have to Export if the domain and range Maps are not the same.
  // Otherwise, this operator would be a permutation, not the identity.
  if (export_.is_null ()) {
    Y.update (alpha, X, beta);
  }
  else {
    if (alpha == STS::one () && beta == STS::zero ()) { // the common case
      Y.doExport (X, *export_, Tpetra::REPLACE);
    }
    else {
      // We know that the domain and range Maps are compatible.  First
      // bring X into the range Map via Export.  Then compute in place
      // in Y.
      MV X_tmp (Y.getMap (), Y.getNumVectors ());
      X_tmp.doExport (X, *export_, Tpetra::REPLACE);
      Y.update (alpha, X_tmp, beta);
    }
  }
  ++numApply_;
}

template <class MatrixType>
int IdentitySolver<MatrixType>::getNumInitialize() const {
  return(numInitialize_);
}

template <class MatrixType>
int IdentitySolver<MatrixType>::getNumCompute() const {
  return(numCompute_);
}

template <class MatrixType>
int IdentitySolver<MatrixType>::getNumApply() const {
  return(numApply_);
}

template <class MatrixType>
double IdentitySolver<MatrixType>::getInitializeTime() const {
  return(initializeTime_);
}

template<class MatrixType>
double IdentitySolver<MatrixType>::getComputeTime() const {
  return(computeTime_);
}

template<class MatrixType>
double IdentitySolver<MatrixType>::getApplyTime() const {
  return(applyTime_);
}

template <class MatrixType>
std::string IdentitySolver<MatrixType>::description () const
{
  std::ostringstream os;

  // Output is a valid YAML dictionary in flow style.  If you don't
  // like everything on a single line, you should call describe()
  // instead.
  os << "\"Ifpack2::IdentitySolver\": {";
  if (this->getObjectLabel () != "") {
    os << "Label: \"" << this->getObjectLabel () << "\", ";
  }
  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
     << "Computed: " << (isComputed () ? "true" : "false") << ", ";

  if (matrix_.is_null ()) {
    os << "Matrix: null";
  }
  else {
    os << "Matrix: not null"
       << ", Global matrix dimensions: ["
       << matrix_->getGlobalNumRows () << ", "
       << matrix_->getGlobalNumCols () << "]";
  }

  os << "}";
  return os.str ();
}

template <class MatrixType>
void IdentitySolver<MatrixType>::
describe (Teuchos::FancyOStream& out,
          const Teuchos::EVerbosityLevel verbLevel) const
{
  using std::endl;
  const Teuchos::EVerbosityLevel vl
    = (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;

  if (vl != Teuchos::VERB_NONE) {
    // By convention, describe() should always begin with a tab.
    Teuchos::OSTab tab0 (out);
    out << "\"Ifpack2::IdentitySolver\":" << endl;
    Teuchos::OSTab tab1 (out);
    out << "MatrixType: " << Teuchos::TypeNameTraits<MatrixType>::name () << endl;
    out << "numInitialize: " << numInitialize_ << endl;
    out << "numCompute: " << numCompute_ << endl;
    out << "numApply: " << numApply_ << endl;
  }
}

template <class MatrixType>
Teuchos::RCP<const typename IdentitySolver<MatrixType>::map_type> IdentitySolver<MatrixType>::getDomainMap() const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver::getDomainMap: "
    "The matrix is null.  Please call setMatrix() with a nonnull input "
    "before calling this method.");
  return matrix_->getDomainMap ();
}

template <class MatrixType>
Teuchos::RCP<const typename IdentitySolver<MatrixType>::map_type> IdentitySolver<MatrixType>::getRangeMap() const
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    matrix_.is_null (), std::runtime_error, "Ifpack2::IdentitySolver::getRangeMap: "
    "The matrix is null.  Please call setMatrix() with a nonnull input "
    "before calling this method.");
  return matrix_->getRangeMap ();
}

template<class MatrixType>
void IdentitySolver<MatrixType>::
setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
{
  // Check in serial or one-process mode if the matrix is square.
  TEUCHOS_TEST_FOR_EXCEPTION(
    ! A.is_null () && A->getComm ()->getSize () == 1 &&
    A->getNodeNumRows () != A->getNodeNumCols (),
    std::runtime_error, "Ifpack2::IdentitySolver::setMatrix: If A's communicator only "
    "contains one process, then A must be square.  Instead, you provided a "
    "matrix A with " << A->getNodeNumRows () << " rows and "
    << A->getNodeNumCols () << " columns.");

  // It's legal for A to be null; in that case, you may not call
  // initialize() until calling setMatrix() with a nonnull input.
  // Regardless, setting the matrix invalidates the preconditioner.
  isInitialized_ = false;
  isComputed_ = false;
  export_ = Teuchos::null;

  matrix_ = A;
}

} // namespace Ifpack2

#define IFPACK2_IDENTITYSOLVER_INSTANT(S,LO,GO,N)                            \
  template class Ifpack2::IdentitySolver< Tpetra::RowMatrix<S, LO, GO, N> >;

#endif // IFPACK2_IDENTITY_SOLVER_DEF_HPP