MultiVector_decl.hpp

#ifndef MULTIVECTOR_DECL_hpp
#define MULTIVECTOR_DECL_hpp
 
/*
 We need to remove some functions from the MultiVector class that are
 only supported for a scalar type (SC) of double. To achieve this is
 simple using concepts in C++20. But for backwards compatibility, we
 offer a C++17 alternative. This alternative is based on SFINAE. A new
 function template parameter is introduced, which is necessary to make
 use of SFINAE. We cannot use substitution failure directly on the
 class template parameters, as this would lead to a compilation error.
 On the function level, if substitution fails, the corresponding
 function is removed. Thus, we can keep SC = double and remove all
 other possibilities. Since the function now has a new template
 parameter, it must be explicitly instantiated if this is done for
 the class as well. The problematic functions, for which only
 SC = double is supported, are functions that return a Thyra vector.
 A Thyra vector does not exist for SC = int, for example.
 TODO: The C++17 compatible code is quite difficult to read.
       I suggest to remove it rather sooner than later when C++17
       compatibility can be dropped.
*/
#if __cplusplus >= 202002L // >= C++20
    #include <concepts>
#else // <= C++17
    #include <type_traits>
#endif
 
#include <MatrixMarket_Tpetra.hpp>
#include <Teuchos_VerboseObject.hpp>
#include <Tpetra_MultiVector.hpp>
#include <Tpetra_Map.hpp>
#include <Tpetra_MultiVector.hpp>
#include <Tpetra_Vector.hpp>
#include <Tpetra_Export.hpp>
#include <Tpetra_Import.hpp>
#include <Thyra_DefaultProductVectorSpace.hpp>
#include <Thyra_TpetraThyraWrappers.hpp>
 
#include "feddlib/core/FEDDCore.hpp"
#include "Map.hpp"
#include "BlockMap.hpp"
#include "BlockMultiVector.hpp"

 
namespace FEDD {
// Forward Declaration of BlockMultiVector to avoid circular includes (BlockMultiVector includes MultiVector) and speed up compilation time.
template <class SC, class LO, class GO, class NO>
class BlockMultiVector;
template <class SC = default_sc, class LO = default_lo, class GO = default_go, class NO = default_no>
// template <class SC = default_sc, class LO = default_lo, class GO = default_go, class NO = default_no>
 
class MultiVector {
 
public:
 
    typedef MultiVector<SC,LO,GO,NO> MultiVector_Type;
    typedef Teuchos::RCP<MultiVector_Type> MultiVectorPtr_Type;
    typedef Teuchos::RCP<const MultiVector_Type> MultiVectorConstPtr_Type;
 
    typedef BlockMultiVector<SC,LO,GO,NO> BlockMultiVector_Type;
    typedef Teuchos::RCP<BlockMultiVector_Type> BlockMultiVectorPtr_Type;
    typedef Teuchos::RCP<const BlockMultiVector_Type> BlockMultiVectorConstPtr_Type;
    
    
    typedef Teuchos::Comm<int> Comm_Type;
    typedef Teuchos::RCP<Comm_Type> CommPtr_Type;    
    typedef Teuchos::RCP<const Comm_Type> CommConstPtr_Type;
 
    // -------------
    typedef Map<LO,GO,NO> Map_Type;
    typedef Teuchos::RCP<Map_Type> MapPtr_Type;
    typedef Teuchos::RCP<const Map_Type> MapConstPtr_Type;
 
    typedef Tpetra::Map<LO,GO,NO> TpetraMap_Type;
    typedef Teuchos::RCP<TpetraMap_Type> TpetraMapPtr_Type;
    typedef Teuchos::RCP<const TpetraMap_Type> TpetraMapConstPtr_Type;
    typedef const TpetraMapConstPtr_Type TpetraMapConstPtrConst_Type;
 
    typedef Tpetra::MultiVector<SC,LO,GO,NO> TpetraMultiVector_Type;
    typedef Teuchos::RCP<TpetraMultiVector_Type> TpetraMultiVectorPtr_Type;
    typedef Teuchos::RCP<const TpetraMultiVector_Type> TpetraMultiVectorConstPtr_Type;
    typedef const TpetraMultiVectorConstPtr_Type TpetraMultiVectorConstPtrConst_Type;
 
    typedef Tpetra::Import<LO,GO,NO> TpetraImport_Type;
    typedef Teuchos::RCP<TpetraImport_Type> TpetraImportPtr_Type;
 
    typedef Tpetra::Export<LO,GO,NO> TpetraExport_Type;
    typedef Teuchos::RCP<TpetraExport_Type> TpetraExportPtr_Type;
 
 
 
    MultiVector();
 
    
    MultiVector( MapConstPtr_Type map, UN nmbVectors=1 );

    MultiVector( TpetraMultiVectorPtr_Type& TpetraMVPtrIn );

    MultiVector( MultiVectorConstPtr_Type mvIn );

    ~MultiVector();

    MultiVector_Type& operator= (const MultiVector_Type& rhs) {
        Tpetra::deep_copy<SC,LO,GO,NO>(*multiVector_,*rhs.getTpetraMultiVector()); // (destination, source)
        return *this;
    }

    bool is_null() const;

    MapConstPtr_Type getMap() const;
    
    MapPtr_Type getMapNonConst();

    TpetraMapConstPtr_Type getMapTpetra() const;

    void replaceGlobalValue (GO globalRow, UN vectorIndex, const SC &value);

    void replaceLocalValue (LO localRow, UN vectorIndex, const SC &value);

    void sumIntoGlobalValue (GO globalRow, UN vectorIndex, const SC &value);
 
    LO getLocalLength() const;

    Teuchos::ArrayRCP< const SC >  getData(UN i) const;

    Teuchos::ArrayRCP< SC > getDataNonConst(UN i) const;

    UN getNumVectors() const;

    void print(Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_EXTREME) const;
 
    TpetraMultiVectorConstPtr_Type getTpetraMultiVector() const;
 
    TpetraMultiVectorPtr_Type getTpetraMultiVectorNonConst();
 
#if __cplusplus >= 202002L // >= C++20
    Teuchos::RCP< Thyra::MultiVectorBase<SC> > getThyraMultiVector( ) requires std::same_as<SC,double>;
 
    Teuchos::RCP<const Thyra::MultiVectorBase<SC> > getThyraMultiVectorConst( ) const requires std::same_as<SC,double>;
#else // <= C++17
    template <typename SCALAR = SC, typename = std::enable_if_t<std::is_same_v<SCALAR,double>>>
    Teuchos::RCP< Thyra::MultiVectorBase<SCALAR> > getThyraMultiVector( );
 
    template <typename SCALAR = SC, typename = std::enable_if_t<std::is_same_v<SCALAR,double>>>
    Teuchos::RCP<const Thyra::MultiVectorBase<SCALAR> > getThyraMultiVectorConst( ) const;
#endif
 
    void fromThyraMultiVector( Teuchos::RCP< Thyra::MultiVectorBase<SC> > thyraMV); 
 
    void norm2(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &norms) const;
 
    void normInf(const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &norms) const;
 
    void dot(MultiVectorConstPtr_Type a, const Teuchos::ArrayView<typename Teuchos::ScalarTraits<SC>::magnitudeType> &dots) const;
 
    // Calculate absolute value of Multivector
    void abs(MultiVectorConstPtr_Type a);
    //this = alpha*A + beta*this
    void update( const SC& alpha, const MultiVector_Type& A, const SC& beta );
 
    //this = alpha*A + beta*B + gamma*this
    void update( const SC& alpha, const MultiVector_Type& A, const SC& beta , const MultiVector_Type& B, const SC& gamma);
 
    // Matrix-matrix multiplication: this = beta*this + alpha*op(A)*op(B).
    void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const SC &alpha, MultiVectorConstPtr_Type &A, MultiVectorConstPtr_Type &B, const SC &beta);
 
#if __cplusplus >= 202002L // >= C++20
    void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const SC &alpha, BlockMultiVectorConstPtr_Type &A, BlockMultiVectorConstPtr_Type &B, const SC &beta) requires std::same_as<SC,double>;
#else // <= C++17
    template <typename SCALAR = SC, typename = std::enable_if_t<std::is_same_v<SCALAR,double>>>
    void multiply(Teuchos::ETransp transA, Teuchos::ETransp transB, const SCALAR &alpha, BlockMultiVectorConstPtr_Type &A, BlockMultiVectorConstPtr_Type &B, const SCALAR &beta);
#endif
    
    void putScalar( const SC& alpha );
 
    void scale( const SC& alpha );
 
    void importFromVector( MultiVectorConstPtr_Type mvIn, bool reuseImport = false, std::string combineMode = "Insert", std::string type="Forward" );
    
    void exportFromVector( MultiVectorConstPtr_Type mvIn, bool reuseExport = false, std::string combineMode = "Insert", std::string type="Forward" );
    
    void writeMM(std::string fileName="mv.mm") const;
    
    void readMM(std::string fileName) const;
    
    MultiVectorConstPtr_Type getVector( int i ) const;
    
    MultiVectorPtr_Type sumColumns() const;
    
    SC getMax() const;
    
private:
 
    TpetraMultiVectorPtr_Type multiVector_;
    MapConstPtr_Type map_;
    TpetraImportPtr_Type importer_;
    TpetraExportPtr_Type exporter_;
};
}
 
#endif