Tpetra_FECrsMatrix_def.hpp

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

#ifndef TPETRA_FECRSMATRIX_DEF_HPP
#define TPETRA_FECRSMATRIX_DEF_HPP
#include <sstream>
#include "Tpetra_CrsMatrix.hpp"

namespace Tpetra {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
    FECrsMatrix(const Teuchos::RCP<const fe_crs_graph_type>& graph,
                const Teuchos::RCP<Teuchos::ParameterList>& params)
  :  // We want the OWNED_PLUS_SHARED graph here
     // NOTE: The casts below are terrible, but necesssary
  crs_matrix_type(graph->inactiveCrsGraph_.is_null() ? Teuchos::rcp_const_cast<crs_graph_type>(Teuchos::rcp_dynamic_cast<const crs_graph_type>(graph)) : graph->inactiveCrsGraph_, params)
  , feGraph_(graph)

{
  const char tfecfFuncName[] = "FECrsMatrix(RCP<const FECrsGraph>[, RCP<ParameterList>]): ";

  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(graph.is_null(), std::runtime_error, "Input graph is null.");
  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(!graph->isFillComplete(), std::runtime_error,
                                        "Input graph is not "
                                        "fill complete. You must call fillComplete on the graph before using "
                                        "it to construct a FECrsMatrix.  Note that calling resumeFill on the "
                                        "graph makes it not fill complete, even if you had previously called "
                                        "fillComplete.  In that case, you must call fillComplete on the graph "
                                        "again.");
  TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(*graph->activeCrsGraph_ != FE::ACTIVE_OWNED, std::runtime_error,
                                        "Input graph must be in FE::ACTIVE_OWNED mode when this constructor is called.");

  bool start_owned = false;
  if (!params.is_null()) {
    if (params->isParameter("start owned")) {
      start_owned = params->get<bool>("start owned", start_owned);
    }
  }
  if (start_owned) {
    activeCrsMatrix_ = Teuchos::rcp(new FE::WhichActive(FE::ACTIVE_OWNED));
  } else {
    activeCrsMatrix_ = Teuchos::rcp(new FE::WhichActive(FE::ACTIVE_OWNED_PLUS_SHARED));
  }

  // Make an "inactive" matrix, if we need to
  if (!graph->inactiveCrsGraph_.is_null()) {
    // We are *requiring* memory aliasing here, so we'll grab the first chunk of the Owned+Shared matrix's values array to make the
    // guy for the Owned matrix.
    inactiveCrsMatrix_ = Teuchos::rcp(new crs_matrix_type(*this, graph));
  }

  fillState_ = Teuchos::rcp(new FE::FillState(FE::FillState::closed));
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doOwnedPlusSharedToOwned(const CombineMode CM) {
  if (!inactiveCrsMatrix_.is_null() && *activeCrsMatrix_ == FE::ACTIVE_OWNED_PLUS_SHARED) {
    // Do a self-export in "restricted mode"
    this->doExport(*this, *feGraph_->ownedRowsImporter_, CM, true);
    inactiveCrsMatrix_->fillComplete();
  }
  crs_matrix_type::fillComplete();
}  // end doOverlapToLocal

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::doOwnedToOwnedPlusShared(const CombineMode /* CM */) {
  // This should be a no-op for all of our purposes
}  // end doLocalToOverlap

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::switchActiveCrsMatrix() {
  if (*activeCrsMatrix_ == FE::ACTIVE_OWNED_PLUS_SHARED)
    *activeCrsMatrix_ = FE::ACTIVE_OWNED;
  else
    *activeCrsMatrix_ = FE::ACTIVE_OWNED_PLUS_SHARED;

  if (inactiveCrsMatrix_.is_null()) return;

  this->swap(*inactiveCrsMatrix_);

}  // end switchActiveCrsMatrix

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::endFill() {
  if (*activeCrsMatrix_ == FE::ACTIVE_OWNED_PLUS_SHARED) {
    doOwnedPlusSharedToOwned(Tpetra::ADD);
    switchActiveCrsMatrix();
  } else
    throw std::runtime_error("FECrsMatrix: Local CrsMatrix already active.  Cannot endFill()");
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::beginFill() {
  // Note: This does not throw an error since the on construction, the FECRS is in overlap mode.  Ergo, calling beginFill(),
  // like one should expect to do in a rational universe, should not cause an error.
  if (*activeCrsMatrix_ == FE::ACTIVE_OWNED) {
    this->resumeFill();
    switchActiveCrsMatrix();
  }
  this->resumeFill();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::beginAssembly() {
  const char tfecfFuncName[] = "FECrsMatrix::beginAssembly: ";
  if (*fillState_ != FE::FillState::closed) {
    std::ostringstream errmsg;
    errmsg << "Cannot begin assembly, matrix is not in a closed state "
           << "but is currently open for "
           << (*fillState_ == FE::FillState::open ? "assembly" : "modification");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  *fillState_ = FE::FillState::open;
  this->beginFill();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::endAssembly() {
  const char tfecfFuncName[] = "FECrsMatrix::endAssembly: ";
  if (*fillState_ != FE::FillState::open) {
    std::ostringstream errmsg;
    errmsg << "Cannot end assembly, matrix is not open for assembly "
           << "but is currently "
           << (*fillState_ == FE::FillState::closed ? "closed" : "open for modification");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  *fillState_ = FE::FillState::closed;
  this->endFill();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::beginModify() {
  const char tfecfFuncName[] = "FECrsMatrix::beginModify: ";
  if (*fillState_ != FE::FillState::closed) {
    std::ostringstream errmsg;
    errmsg << "Cannot begin modifying, matrix is not in a closed state "
           << "but is currently open for "
           << (*fillState_ == FE::FillState::open ? "assembly" : "modification");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  *fillState_ = FE::FillState::modify;
  this->resumeFill();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::endModify() {
  const char tfecfFuncName[] = "FECrsMatrix::endModify: ";
  if (*fillState_ != FE::FillState::modify) {
    std::ostringstream errmsg;
    errmsg << "Cannot end modifying, matrix is not open to modify but is currently "
           << (*fillState_ == FE::FillState::open ? "open for assembly" : "closed");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  *fillState_ = FE::FillState::closed;
  this->fillComplete();
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceGlobalValuesImpl(
    impl_scalar_type rowVals[],
    const crs_graph_type& graph,
    const RowInfo& rowInfo,
    const GlobalOrdinal inds[],
    const impl_scalar_type newVals[],
    const LocalOrdinal numElts) {
  const char tfecfFuncName[] = "FECrsMatrix::replaceGlobalValues: ";
  if (*fillState_ != FE::FillState::open) {
    std::ostringstream errmsg;
    errmsg << "Cannot replace global values, matrix is not open for assembly "
           << "but is currently "
           << (*fillState_ == FE::FillState::modify ? "open for modification" : "closed");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  return crs_matrix_type::replaceGlobalValuesImpl(rowVals, graph, rowInfo, inds, newVals, numElts);
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceLocalValuesImpl(
    impl_scalar_type rowVals[],
    const crs_graph_type& graph,
    const RowInfo& rowInfo,
    const LocalOrdinal inds[],
    const impl_scalar_type newVals[],
    const LocalOrdinal numElts) {
  const char tfecfFuncName[] = "FECrsMatrix::replaceLocalValues: ";
  if (*fillState_ != FE::FillState::open && *fillState_ != FE::FillState::modify) {
    std::ostringstream errmsg;
    errmsg << "Cannot replace local values, matrix is not open to fill/modify. "
           << "The matrix is currently closed";
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  return crs_matrix_type::replaceLocalValuesImpl(rowVals, graph, rowInfo, inds, newVals, numElts);
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sumIntoGlobalValuesImpl(
    impl_scalar_type rowVals[],
    const crs_graph_type& graph,
    const RowInfo& rowInfo,
    const GlobalOrdinal inds[],
    const impl_scalar_type newVals[],
    const LocalOrdinal numElts,
    const bool atomic) {
  const char tfecfFuncName[] = "FECrsMatrix::sumIntoGlobalValues: ";
  if (*fillState_ != FE::FillState::open) {
    std::ostringstream errmsg;
    errmsg << "Cannot sum in to global values, matrix is not open for assembly. "
           << "The matrix is currently "
           << (*fillState_ == FE::FillState::modify ? "open for modification" : "closed");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  return crs_matrix_type::sumIntoGlobalValuesImpl(
      rowVals, graph, rowInfo, inds, newVals, numElts, atomic);
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
LocalOrdinal
FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sumIntoLocalValuesImpl(
    impl_scalar_type rowVals[],
    const crs_graph_type& graph,
    const RowInfo& rowInfo,
    const LocalOrdinal inds[],
    const impl_scalar_type newVals[],
    const LocalOrdinal numElts,
    const bool atomic) {
  const char tfecfFuncName[] = "FECrsMatrix::sumIntoLocalValues: ";
  if (*fillState_ != FE::FillState::open) {
    std::ostringstream errmsg;
    errmsg << "Cannot sum in to local values, matrix is not open for assembly. "
           << "The matrix is currently "
           << (*fillState_ == FE::FillState::modify ? "open for modification" : "closed");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  return crs_matrix_type::sumIntoLocalValuesImpl(
      rowVals, graph, rowInfo, inds, newVals, numElts, atomic);
}

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void FECrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::insertGlobalValuesImpl(
    crs_graph_type& graph,
    RowInfo& rowInfo,
    const GlobalOrdinal gblColInds[],
    const impl_scalar_type vals[],
    const size_t numInputEnt) {
  const char tfecfFuncName[] = "FECrsMatrix::insertGlobalValues: ";
  if (*fillState_ != FE::FillState::open) {
    std::ostringstream errmsg;
    errmsg << "Cannot insert global values, matrix is not open for assembly. "
           << "The matrix is currently "
           << (*fillState_ == FE::FillState::modify ? "open for modification" : "closed");
    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(true, std::logic_error, errmsg.str());
  }
  return crs_matrix_type::insertGlobalValuesImpl(graph, rowInfo, gblColInds, vals, numInputEnt);
}

}  // end namespace Tpetra

//
// Explicit instantiation macro
//
// Must be expanded from within the Tpetra namespace!
//
#define TPETRA_FECRSMATRIX_INSTANT(SCALAR, LO, GO, NODE) \
  template class FECrsMatrix<SCALAR, LO, GO, NODE>;

#endif  // TPETRA_FECRSMATRIX_DEF