MueLu_TpetraOperatorAsRowMatrix.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//        MueLu: A package for multigrid based preconditioning
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// 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
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER
#include <Tpetra_RowMatrix.hpp>

namespace MueLu {

template <class Scalar        = Tpetra::Operator<>::scalar_type,
          class LocalOrdinal  = typename Tpetra::Operator<Scalar>::local_ordinal_type,
          class GlobalOrdinal = typename Tpetra::Operator<Scalar, LocalOrdinal>::global_ordinal_type,
          class Node          = typename Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal>::node_type>
class TpetraOperatorAsRowMatrix : public Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
 public:
  using op_type  = Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
  using vec_type = Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>;

  using row_matrix_type = Tpetra::RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>;

  using impl_scalar_type = typename row_matrix_type::impl_scalar_type;
  using mag_type         = typename Kokkos::ArithTraits<impl_scalar_type>::mag_type;

  using local_inds_device_view_type =
      typename row_matrix_type::local_inds_device_view_type;
  using local_inds_host_view_type =
      typename row_matrix_type::local_inds_host_view_type;
  using nonconst_local_inds_host_view_type =
      typename row_matrix_type::nonconst_local_inds_host_view_type;

  using global_inds_device_view_type =
      typename row_matrix_type::global_inds_device_view_type;
  using global_inds_host_view_type =
      typename row_matrix_type::global_inds_host_view_type;
  using nonconst_global_inds_host_view_type =
      typename row_matrix_type::nonconst_global_inds_host_view_type;

  using values_device_view_type =
      typename row_matrix_type::values_device_view_type;
  using values_host_view_type =
      typename row_matrix_type::values_host_view_type;
  using nonconst_values_host_view_type =
      typename row_matrix_type::nonconst_values_host_view_type;



  TpetraOperatorAsRowMatrix(const RCP<op_type>& op)
    : op_(op)
    , diag_(Teuchos::null) {}

  TpetraOperatorAsRowMatrix(const RCP<op_type>& op,
                            const RCP<vec_type>& diag)
    : op_(op)
    , diag_(diag) {}

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getDomainMap() const {
    return op_->getDomainMap();
  }

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getRangeMap() const {
    return op_->getRangeMap();
  }


  void apply(const Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& X,
             Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& Y,
             Teuchos::ETransp mode = Teuchos::NO_TRANS,
             Scalar alpha          = Teuchos::ScalarTraits<Scalar>::one(),
             Scalar beta           = Teuchos::ScalarTraits<Scalar>::zero()) const {
    op_->apply(X, Y, mode, alpha, beta);
  }

  // Fake RowMatrix interface
  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getRowMap() const {
    return op_->getRangeMap();
  }

  Teuchos::RCP<const Tpetra::Map<LocalOrdinal, GlobalOrdinal, Node> > getColMap() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  typename row_matrix_type::local_ordinal_type getBlockSize() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  Teuchos::RCP<const Teuchos::Comm<int> > getComm() const {
    return op_->getDomainMap()->getComm();
  }

  Teuchos::RCP<const Tpetra::RowGraph<LocalOrdinal, GlobalOrdinal, Node> > getGraph() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  Tpetra::global_size_t getGlobalNumRows() const {
    return getRowMap()->getGlobalNumElements();
  }

  Tpetra::global_size_t getGlobalNumCols() const {
    return getDomainMap()->getGlobalNumElements();
  }

  size_t getLocalNumRows() const {
    return getRowMap()->getLocalNumElements();
  }

  size_t getLocalNumCols() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  GlobalOrdinal getIndexBase() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  Tpetra::global_size_t getGlobalNumEntries() const {
    return 0;
  }

  size_t getLocalNumEntries() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  size_t getGlobalMaxNumRowEntries() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  size_t getLocalMaxNumRowEntries() const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  bool hasColMap() const {
    return false;
  }

  bool isLocallyIndexed() const {
    return true;
  }

  bool isGloballyIndexed() const {
    return true;
  }

  bool isFillComplete() const {
    return true;
  }

  bool supportsRowViews() const {
    return false;
  }

  void
  getGlobalRowCopy(GlobalOrdinal GlobalRow,
                   nonconst_global_inds_host_view_type& Indices,
                   nonconst_values_host_view_type& Values,
                   size_t& NumEntries) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  void
  getLocalRowCopy(LocalOrdinal LocalRow,
                  nonconst_local_inds_host_view_type& Indices,
                  nonconst_values_host_view_type& Values,
                  size_t& NumEntries) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  void
  getGlobalRowView(GlobalOrdinal GlobalRow,
                   global_inds_host_view_type& indices,
                   values_host_view_type& values) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  void
  getLocalRowView(LocalOrdinal LocalRow,
                  local_inds_host_view_type& indices,
                  values_host_view_type& values) const {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  void getLocalDiagCopy(Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& diag) const {
    if (diag_.is_null())
      throw MueLu::Exceptions::RuntimeError("No diagonal available.");
    else
      diag = *diag_;
  }

  void leftScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  void rightScale(const Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x) {
    throw MueLu::Exceptions::RuntimeError("Not implemented.");
  }

  mag_type getFrobeniusNorm() const {
    return 0.;
  }

  // void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel) const {
  //   using std::setw;
  //   using std::endl;
  //   const size_t numRows = nearField_->getRowMap()->getGlobalNumElements();
  //   const size_t nnzNearField = nearField_->getGlobalNumEntries();
  //   const double nnzNearPerRow = Teuchos::as<double>(nnzNearField)/numRows;
  //   const size_t nnzKernelApprox = kernelApproximations_->pointA_->getGlobalNumEntries();
  //   const size_t numClusterPairs = kernelApproximations_->blockA_->getGlobalNumEntries();
  //   const size_t nnzBasis = basisMatrix_->getGlobalNumEntries();
  //   size_t nnzTransfer = 0;
  //   for (size_t i = 0; i<transferMatrices_.size(); i++)
  //     nnzTransfer += transferMatrices_[i]->pointA_->getGlobalNumEntries();
  //   const size_t nnzTotal = nnzNearField+nnzKernelApprox+nnzBasis+nnzTransfer;
  //   const double nnzTotalPerRow = Teuchos::as<double>(nnzTotal)/numRows;
  //   std::ostringstream oss;
  //   oss << std::left;
  //   oss << setw(9) << "rows"  << setw(12) << "nnz(near)"  << setw(14) << "nnz(near)/row" << setw(12) << "nnz(basis)" << setw(15) << "#cluster pairs" << setw(12)<< "nnz(kernel)"    << setw(14) << "nnz(transfer)" << setw(12) << "nnz(total)" << setw(14) << "nnz(total)/row" << endl;
  //   oss << setw(9) << numRows << setw(12) << nnzNearField << setw(14) << nnzNearPerRow   << setw(12) << nnzBasis     << setw(15) << numClusterPairs  << setw(12) << nnzKernelApprox << setw(14) << nnzTransfer     << setw(12) << nnzTotal     << setw(14) << nnzTotalPerRow   << endl;
  //   out << oss.str();
  // }

 private:
  RCP<op_type> op_;
  RCP<vec_type> diag_;
};

}  // namespace MueLu