Tpetra_LinearProblem_def.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//          Tpetra: Templated Linear Algebra Services Package
//
// Copyright 2008 NTESS and the Tpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef TPETRA_LINEARPROBLEM_DEF_HPP
#define TPETRA_LINEARPROBLEM_DEF_HPP


#include "Teuchos_DataAccess.hpp"
#include "Teuchos_TestForException.hpp"
#include "Tpetra_Details_Behavior.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_MultiVector_decl.hpp"

namespace Tpetra {

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  LinearProblem ()
    : dist_object_type (Teuchos::rcp (new map_type ())),
      A_(Teuchos::null),
      X_(Teuchos::null),
      B_(Teuchos::null)
  {
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  LinearProblem (const Teuchos::RCP<row_matrix_type> & A,
                 const Teuchos::RCP<multivector_type>& X,
                 const Teuchos::RCP<multivector_type>& B)
    : dist_object_type (A->getDomainMap()),
      A_(A),
      X_(X),
      B_(B)
  {
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  LinearProblem (const LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Problem)
    : dist_object_type (Problem),
      A_(Problem.A_),
      X_(Problem.X_),
      B_(Problem.B_)
  {
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  leftScale(const Teuchos::RCP<const vector_type> & D, Teuchos::ETransp mode)
  {
    const Scalar ST0 = Teuchos::ScalarTraits<Scalar>::zero();
    const Scalar ST1 = Teuchos::ScalarTraits<Scalar>::one();
    if (mode == Teuchos::NO_TRANS) {
      A_->leftScale(*D);
      B_->elementWiseMultiply(ST1, *D, *B_, ST0);
    }
    else {
      A_->rightScale(*D);
      vector_type R(*D, Teuchos::DataAccess::Copy);
      R.reciprocal(*D);
      X_->elementWiseMultiply(ST1, R, *X_, ST0);
    }
  
    return;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  rightScale(const Teuchos::RCP<const vector_type> & D, Teuchos::ETransp mode)
  {
    const Scalar ST0 = Teuchos::ScalarTraits<Scalar>::zero();
    const Scalar ST1 = Teuchos::ScalarTraits<Scalar>::one();
    if (mode == Teuchos::NO_TRANS) {
      A_->rightScale(*D);
      vector_type R(*D, Teuchos::DataAccess::Copy);
      R.reciprocal(*D);
      X_->elementWiseMultiply(ST1, R, *X_, ST0);
    }
    else {
      A_->leftScale(*D);
      B_->elementWiseMultiply(ST1, *D, *B_, ST0);
    }
    return;
  }

  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  void LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  checkInput() const {

    const char tfecfFuncName[] = "checkInput: ";

    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(A_==Teuchos::null, std::runtime_error, 
      "Linear problem does not have a matrix (A_).");
  
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(X_==Teuchos::null, 
      std::logic_error, "Solution vector (X_) is unset.");

    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(B_==Teuchos::null, 
      std::logic_error, "RHS vector (B_) is unset.");
  
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!A_->getRowMap()->isSameAs(*(X_->getMap())),
      std::logic_error, "Domain map of matrix is not the 'same as' the solution map.");

    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!A_->getRowMap()->isSameAs(*(B_->getMap())),
      std::logic_error, "Range map of matrix is not the 'same as' the RHS map.");

    return;
  }

  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
  bool LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
  checkSizes (const SrcDistObject& sourceObj)
  {
    // Check whether the source object is a LinearProblem.  If not, then
    // we can't even compare sizes.
    typedef LinearProblem<Scalar, LocalOrdinal, GlobalOrdinal, Node> LP;
    const LP* src = dynamic_cast<const LP*> (&sourceObj);
    if (src == nullptr) {
      return false;
    }
    else {
      this->checkInput();
      src->checkInput();

      return ((this->A_->getDomainMap() == src->getMatrix()->getDomainMap()) and
              (this->A_->getRangeMap() == src->getMatrix()->getRangeMap()));
    }
  }

} // namespace Tpetra

//
// Explicit instantiation macro
//
// Must be expanded from within the Tpetra namespace!
//

#define TPETRA_LINEARPROBLEM_INSTANT(SCALAR,LO,GO,NODE) \
  template class LinearProblem< SCALAR , LO , GO , NODE >;


#endif // TPETRA_LINEARPROBLEM_DEF_HPP