Xpetra_TripleMatrixMultiply_decl.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//             Xpetra: A linear algebra interface package
//
// Copyright 2012 NTESS and the Xpetra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_DECL_HPP_
#define PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_DECL_HPP_
 
#include "Xpetra_ConfigDefs.hpp"
 
// #include "Xpetra_BlockedCrsMatrix.hpp"
#include "Xpetra_CrsMatrixWrap.hpp"
#include "Xpetra_MapExtractor.hpp"
#include "Xpetra_Map.hpp"
#include "Xpetra_MatrixFactory.hpp"
#include "Xpetra_Matrix.hpp"
#include "Xpetra_StridedMapFactory.hpp"
#include "Xpetra_StridedMap.hpp"
#include "Xpetra_IO.hpp"
 
#ifdef HAVE_XPETRA_TPETRA
#include <TpetraExt_TripleMatrixMultiply.hpp>
#include <Xpetra_TpetraCrsMatrix.hpp>
#include <Tpetra_BlockCrsMatrix.hpp>
#include <Tpetra_BlockCrsMatrix_Helpers.hpp>
// #include <Xpetra_TpetraMultiVector.hpp>
// #include <Xpetra_TpetraVector.hpp>
#endif  // HAVE_XPETRA_TPETRA
 
namespace Xpetra {
 
template <class Scalar,
          class LocalOrdinal /*= int*/,
          class GlobalOrdinal /*= LocalOrdinal*/,
          class Node /*= Tpetra::KokkosClassic::DefaultNode::DefaultNodeType*/>
class TripleMatrixMultiply {
#undef XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
#include "Xpetra_UseShortNames.hpp"
 
 public:
  static void MultiplyRAP(const Matrix& R, bool transposeR,
                          const Matrix& A, bool transposeA,
                          const Matrix& P, bool transposeP,
                          Matrix& Ac,
                          bool call_FillComplete_on_result = true,
                          bool doOptimizeStorage           = true,
                          const std::string& label         = std::string(),
                          const RCP<ParameterList>& params = null);
 
};  // class TripleMatrixMultiply
 
#ifdef HAVE_XPETRA_EPETRA
// specialization TripleMatrixMultiply for SC=double, LO=GO=int
template <>
class TripleMatrixMultiply<double, int, int, EpetraNode> {
  typedef double Scalar;
  typedef int LocalOrdinal;
  typedef int GlobalOrdinal;
  typedef EpetraNode Node;
#include "Xpetra_UseShortNames.hpp"
 
 public:
  static void MultiplyRAP(const Matrix& R, bool transposeR,
                          const Matrix& A, bool transposeA,
                          const Matrix& P, bool transposeP,
                          Matrix& Ac,
                          bool call_FillComplete_on_result = true,
                          bool doOptimizeStorage           = true,
                          const std::string& label         = std::string(),
                          const RCP<ParameterList>& params = null) {
    TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
                               Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
    TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
                               Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
 
    TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
    TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
    TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
 
    bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
 
    if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
      throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
    } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
#ifdef HAVE_XPETRA_TPETRA
#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_INT))) || \
     (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_INT))))
      throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,int> ETI enabled."));
#else
      using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
      if (helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P) && helpers::isTpetraCrs(Ac)) {
        // All matrices are Crs
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(R);
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(P);
        Tpetra::CrsMatrix<SC, LO, GO, NO>& tpAc      = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(Ac);
 
        // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
        // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
        Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
      } else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
        // All matrices are BlockCrs (except maybe Ac)
        // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
        // we'll need to think about refactoring BlockCrs so we can do something smarter here.
        if (!A.getRowMap()->getComm()->getRank())
          std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
 
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(R);
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(P);
        //          Tpetra::BlockCrsMatrix<SC,LO,GO,NO> &       tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
 
        using CRS           = Tpetra::CrsMatrix<SC, LO, GO, NO>;
        RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
        RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
        RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
        //          RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
 
        // FIXME: The lines below only works because we're assuming Ac is Point
        RCP<CRS> Accrs              = Teuchos::rcp(new CRS(Rcrs->getRowMap(), 0));
        const bool do_fill_complete = true;
        Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
 
        // Temporary output matrix
        RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t       = Tpetra::convertToBlockCrsMatrix(*Accrs, A.GetStorageBlockSize());
        RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
        RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p            = Ac_x;
 
        // We can now cheat and replace the innards of Ac
        RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(Ac));
        Ac_w->replaceCrsMatrix(Ac_p);
 
      } else {
        // Mix and match (not supported)
        TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
      }
#endif
#else
      throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
#endif
      if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
        RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
        fillParams->set("Optimize Storage", doOptimizeStorage);
        Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
                        (transposeR) ? R.getDomainMap() : R.getRangeMap(),
                        fillParams);
      }
 
      // transfer striding information
      RCP<const Map> domainMap = Teuchos::null;
      RCP<const Map> rangeMap  = Teuchos::null;
 
      const std::string stridedViewLabel("stridedMaps");
      const size_t blkSize = 1;
      std::vector<size_t> stridingInfo(1, blkSize);
      LocalOrdinal stridedBlockId = -1;
 
      if (R.IsView(stridedViewLabel)) {
        rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
      } else {
        rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
        rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
      }
 
      if (P.IsView(stridedViewLabel)) {
        domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
      } else {
        domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
        domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
      }
      Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
    }
 
  }  // end Multiply
 
};  // end specialization on SC=double, GO=int and NO=EpetraNode
 
// specialization TripleMatrixMultiply for SC=double, GO=long long and NO=EpetraNode
template <>
class TripleMatrixMultiply<double, int, long long, EpetraNode> {
  typedef double Scalar;
  typedef int LocalOrdinal;
  typedef long long GlobalOrdinal;
  typedef EpetraNode Node;
#include "Xpetra_UseShortNames.hpp"
 
 public:
  static void MultiplyRAP(const Matrix& R, bool transposeR,
                          const Matrix& A, bool transposeA,
                          const Matrix& P, bool transposeP,
                          Matrix& Ac,
                          bool call_FillComplete_on_result = true,
                          bool doOptimizeStorage           = true,
                          const std::string& label         = std::string(),
                          const RCP<ParameterList>& params = null) {
    TEUCHOS_TEST_FOR_EXCEPTION(transposeR == false && Ac.getRowMap()->isSameAs(*R.getRowMap()) == false,
                               Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as row map of R");
    TEUCHOS_TEST_FOR_EXCEPTION(transposeR == true && Ac.getRowMap()->isSameAs(*R.getDomainMap()) == false,
                               Exceptions::RuntimeError, "XpetraExt::TripleMatrixMultiply::MultiplyRAP: row map of Ac is not same as domain map of R");
 
    TEUCHOS_TEST_FOR_EXCEPTION(!R.isFillComplete(), Exceptions::RuntimeError, "R is not fill-completed");
    TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Exceptions::RuntimeError, "A is not fill-completed");
    TEUCHOS_TEST_FOR_EXCEPTION(!P.isFillComplete(), Exceptions::RuntimeError, "P is not fill-completed");
 
    bool haveMultiplyDoFillComplete = call_FillComplete_on_result && doOptimizeStorage;
 
    if (Ac.getRowMap()->lib() == Xpetra::UseEpetra) {
      throw(Xpetra::Exceptions::RuntimeError("Xpetra::TripleMatrixMultiply::MultiplyRAP is only implemented for Tpetra"));
    } else if (Ac.getRowMap()->lib() == Xpetra::UseTpetra) {
#ifdef HAVE_XPETRA_TPETRA
#if ((defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_OPENMP) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))) || \
     (!defined(EPETRA_HAVE_OMP) && (!defined(HAVE_TPETRA_INST_SERIAL) || !defined(HAVE_TPETRA_INST_INT_LONG_LONG))))
      throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra <double,int,long long,EpetraNode> ETI enabled."));
#else
      using helpers = Xpetra::Helpers<SC, LO, GO, NO>;
      if (helpers::isTpetraCrs(R) && helpers::isTpetraCrs(A) && helpers::isTpetraCrs(P)) {
        // All matrices are Crs
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(R);
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(A);
        const Tpetra::CrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraCrs(P);
        Tpetra::CrsMatrix<SC, LO, GO, NO>& tpAc      = Xpetra::Helpers<SC, LO, GO, NO>::Op2NonConstTpetraCrs(Ac);
 
        // 18Feb2013 JJH I'm reenabling the code that allows the matrix matrix multiply to do the fillComplete.
        // Previously, Tpetra's matrix matrix multiply did not support fillComplete.
        Tpetra::TripleMatrixMultiply::MultiplyRAP(tpR, transposeR, tpA, transposeA, tpP, transposeP, tpAc, haveMultiplyDoFillComplete, label, params);
      } else if (helpers::isTpetraBlockCrs(R) && helpers::isTpetraBlockCrs(A) && helpers::isTpetraBlockCrs(P)) {
        // All matrices are BlockCrs (except maybe Ac)
        // FIXME: For the moment we're just going to clobber the innards of AC, so no reuse. Once we have a reuse kernel,
        // we'll need to think about refactoring BlockCrs so we can do something smarter here.
        if (!A.getRowMap()->getComm()->getRank())
          std::cout << "WARNING: Using inefficient BlockCrs Multiply Placeholder" << std::endl;
 
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpR = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(R);
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpA = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(A);
        const Tpetra::BlockCrsMatrix<SC, LO, GO, NO>& tpP = Xpetra::Helpers<SC, LO, GO, NO>::Op2TpetraBlockCrs(P);
        //          Tpetra::BlockCrsMatrix<SC,LO,GO,NO> &       tpAc = Xpetra::Helpers<SC,LO,GO,NO>::Op2NonConstTpetraBlockCrs(Ac);
 
        using CRS           = Tpetra::CrsMatrix<SC, LO, GO, NO>;
        RCP<const CRS> Rcrs = Tpetra::convertToCrsMatrix(tpR);
        RCP<const CRS> Acrs = Tpetra::convertToCrsMatrix(tpA);
        RCP<const CRS> Pcrs = Tpetra::convertToCrsMatrix(tpP);
        //          RCP<CRS> Accrs = Tpetra::convertToCrsMatrix(tpAc);
 
        // FIXME: The lines below only works because we're assuming Ac is Point
        RCP<CRS> Accrs              = Teuchos::rcp(new CRS(Rcrs->getRowMap(), 0));
        const bool do_fill_complete = true;
        Tpetra::TripleMatrixMultiply::MultiplyRAP(*Rcrs, transposeR, *Acrs, transposeA, *Pcrs, transposeP, *Accrs, do_fill_complete, label, params);
 
        // Temporary output matrix
        RCP<Tpetra::BlockCrsMatrix<SC, LO, GO, NO> > Ac_t       = Tpetra::convertToBlockCrsMatrix(*Accrs, A.GetStorageBlockSize());
        RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > Ac_x = Teuchos::rcp(new Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO>(Ac_t));
        RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > Ac_p            = Ac_x;
 
        // We can now cheat and replace the innards of Ac
        RCP<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> > Ac_w = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC, LO, GO, NO> >(Teuchos::rcpFromRef(Ac));
        Ac_w->replaceCrsMatrix(Ac_p);
      } else {
        // Mix and match
        TEUCHOS_TEST_FOR_EXCEPTION(1, Exceptions::RuntimeError, "Mix-and-match Crs/BlockCrs Multiply not currently supported");
      }
 
#endif
#else
      throw(Xpetra::Exceptions::RuntimeError("Xpetra must be compiled with Tpetra."));
#endif
      if (call_FillComplete_on_result && !haveMultiplyDoFillComplete) {
        RCP<Teuchos::ParameterList> fillParams = rcp(new Teuchos::ParameterList());
        fillParams->set("Optimize Storage", doOptimizeStorage);
        Ac.fillComplete((transposeP) ? P.getRangeMap() : P.getDomainMap(),
                        (transposeR) ? R.getDomainMap() : R.getRangeMap(),
                        fillParams);
      }
 
      // transfer striding information
      RCP<const Map> domainMap = Teuchos::null;
      RCP<const Map> rangeMap  = Teuchos::null;
 
      const std::string stridedViewLabel("stridedMaps");
      const size_t blkSize = 1;
      std::vector<size_t> stridingInfo(1, blkSize);
      LocalOrdinal stridedBlockId = -1;
 
      if (R.IsView(stridedViewLabel)) {
        rangeMap = transposeR ? R.getColMap(stridedViewLabel) : R.getRowMap(stridedViewLabel);
      } else {
        rangeMap = transposeR ? R.getDomainMap() : R.getRangeMap();
        rangeMap = StridedMapFactory::Build(rangeMap, stridingInfo, stridedBlockId);
      }
 
      if (P.IsView(stridedViewLabel)) {
        domainMap = transposeP ? P.getRowMap(stridedViewLabel) : P.getColMap(stridedViewLabel);
      } else {
        domainMap = transposeP ? P.getRangeMap() : P.getDomainMap();
        domainMap = StridedMapFactory::Build(domainMap, stridingInfo, stridedBlockId);
      }
      Ac.CreateView(stridedViewLabel, rangeMap, domainMap);
    }
 
  }  // end Multiply
 
};  // end specialization on GO=long long and NO=EpetraNode
#endif
 
}  // end namespace Xpetra
 
#define XPETRA_TRIPLEMATRIXMULTIPLY_SHORT
 
#endif /* PACKAGES_XPETRA_SUP_UTILS_XPETRA_TRIPLEMATRIXMULTIPLY_DECL_HPP_ */