BelosTeuchosDenseAdapter.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                 Belos: Block Linear Solvers Package
//
// Copyright 2004-2016 NTESS and the Belos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
//
#ifndef BELOS_TEUCHOS_DENSE_MAT_TRAITS_HPP
#define BELOS_TEUCHOS_DENSE_MAT_TRAITS_HPP
 
#include "Teuchos_RCP.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_SerialDenseSolver.hpp"
#include "Teuchos_SerialSpdDenseSolver.hpp"
 
#include "BelosDenseMatTraits.hpp" 
#include "BelosDenseSolver.hpp"
 
namespace Belos {
 
  template<class ScalarType>
  class TeuchosDenseSolver : public DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType> >
  {
  public:
 
 
 
    TeuchosDenseSolver() {}
 
    virtual ~TeuchosDenseSolver() {}
 
 
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::setMatrix;
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::setVectors;
 
 
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::setSPD;
 
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::factorWithEquilibration;
 
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::solveWithTransposeFlag;
 
 
 
 
    int factor() 
    {
       //std::cout << "Calling DenseMatrix<Teuchos> factor!" << std::endl;
       // Create solver object if one doesn't exist.
       if ( (dSolver_==Teuchos::null) && (spdSolver_==Teuchos::null) )
       {
         if (spd_)
           spdSolver_ = Teuchos::rcp( new Teuchos::SerialSpdDenseSolver<int,ScalarType>() );
         else
           dSolver_ = Teuchos::rcp( new Teuchos::SerialDenseSolver<int,ScalarType>() );
      }
 
      int info = 0;
 
      // Set the matrix and factor
      if (spd_)
      {
        if (newMatrix_)
        {
          spdMatrix_ = Teuchos::rcp( new Teuchos::SerialSymDenseMatrix<int,ScalarType>(Teuchos::View, true, A_->values(),
                                                                                       A_->numRows(), A_->numCols()) );
        } 
        spdSolver_->setMatrix( spdMatrix_ );
        spdSolver_->factorWithEquilibration( equilibrate_ );
        info = spdSolver_->factor();
      }
      else
      {
        dSolver_->setMatrix( A_ );
        dSolver_->solveWithTransposeFlag( TRANS_ );
        dSolver_->factorWithEquilibration( equilibrate_ );
        info = dSolver_->factor();
      }     
 
      if (!info) 
        newMatrix_ = false;
      return info;
    }
 
 
    int solve() 
    {
      //std::cout << "Calling DenseMatrix<Teuchos> solve!" << std::endl;
      // Check if this is a new matrix that has not been factored.
      if (newMatrix_)
        factor();
 
      int info = 0;
 
      if (spd_)
      {
        spdSolver_->setVectors( X_, B_ );
        info = spdSolver_->solve();
      }
      else
      {
        dSolver_->setVectors( X_, B_ );
        info = dSolver_->solve();
      }
 
      return info;
    }
 
 
  private:
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::newMatrix_;
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::equilibrate_;
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::TRANS_;
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::spd_;
 
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::A_;
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::X_;
    using DenseSolver<ScalarType,Teuchos::SerialDenseMatrix<int,ScalarType>>::B_;
 
    // Pointer to dense solver for general linear systems
    Teuchos::RCP<Teuchos::SerialDenseSolver<int,ScalarType>> dSolver_;
 
    // Pointer to dense solver for SPD linear systems
    Teuchos::RCP<Teuchos::SerialSymDenseMatrix<int,ScalarType>> spdMatrix_;
    Teuchos::RCP<Teuchos::SerialSpdDenseSolver<int,ScalarType>> spdSolver_;
 
  };
 
 
  template<class ScalarType>
  class DenseMatTraits<ScalarType, Teuchos::SerialDenseMatrix<int,ScalarType>>{
  public:
    
 
    static Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType>> Create() { 
      return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>());
    }     
 
    static Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType>> Create( const int numrows, const int numcols, bool initZero = true) { 
      return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(numrows,numcols,initZero));
    }     
 
    static Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType>> CreateCopy(const Teuchos::SerialDenseMatrix<int,ScalarType> & dm, bool transpose=false) {
      if (transpose)    
        return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(dm, Teuchos::CONJ_TRANS));
      else 
        return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, dm));
    }
 
    static ScalarType* GetRawHostPtr( Teuchos::SerialDenseMatrix<int,ScalarType> & dm )
    { return dm.values(); }     
 
    static ScalarType const * GetConstRawHostPtr(const Teuchos::SerialDenseMatrix<int,ScalarType> & dm )  
    { return const_cast<ScalarType const *>(dm.values()); }     
 
    //static void RawPtrDataModified(Teuchos::SerialDenseMatrix<int,ScalarType> & dm ) { }
 
    //        Row and column indexing is zero-based.
    //        Source  - Reference to another dense matrix from which values are to be copied.
    //        numRows - The number of rows in this matrix.
    //        numCols - The number of columns in this matrix.
    //        startRow  - The row of Source from which the submatrix copy should start.
    //        startCol  - The column of Source from which the submatrix copy should start.
    //  
    //        Should ints be const? Should they be ints or some other ordinal type?
    static Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType>> Subview(Teuchos::SerialDenseMatrix<int,ScalarType> & source, int numRows, int numCols, int startRow=0, int startCol=0)
    { return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View, source, numRows, numCols, startRow, startCol)); }    
 
    static Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType>> SubviewConst( 
                              const Teuchos::SerialDenseMatrix<int,ScalarType> & source, int numRows, int numCols, int startRow=0, int startCol=0)
    { return Teuchos::rcp(new const Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View, source, numRows, numCols, startRow, startCol)); }    
 
    static Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType>> SubviewCopy( const Teuchos::SerialDenseMatrix<int,ScalarType>& source, int numRows, int numCols, int startRow=0, int startCol=0)
    { return Teuchos::rcp(new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::Copy, source, numRows, numCols, startRow, startCol)); }    
 
 
    static int GetNumRows( const Teuchos::SerialDenseMatrix<int,ScalarType>& dm )
    { return dm.numRows(); }     
 
    static int GetNumCols( const Teuchos::SerialDenseMatrix<int,ScalarType>& dm )
    { return dm.numCols(); }     
 
    static int GetStride( const Teuchos::SerialDenseMatrix<int,ScalarType>& dm )
    { return dm.stride(); }    
 
 
 
    /* \brief Reshaping method for changing the size of \c dm to have \c numrows rows and \c numcols columns.
     *        All values will be initialized to zero if the final argument is true. 
     *        If the final argument is fale, the previous entries in 
     *        the matrix will be maintained. For new entries that did not exist in the previous matrix, values will
     *        contain noise from memory. 
    */
    static void Reshape( Teuchos::SerialDenseMatrix<int,ScalarType>& dm, const int numrows, const int numcols, bool initZero = false) { 
      if(initZero){
        int err =  dm.shape(numrows,numcols); 
        if(err != 0){throw std::runtime_error ("Error in DenseMatrixTraits::shape. Teuchos::SerialDenseMatrix.shape failed.");}
      }
      else{
        int err =  dm.reshape(numrows,numcols); 
        if(err != 0){throw std::runtime_error ("Error in DenseMatrixTraits::reshape. Teuchos::SerialDenseMatrix.reshape failed.");}
      }
    }     
 
 
 
    static ScalarType & Value( Teuchos::SerialDenseMatrix<int,ScalarType>& dm, const int i, const int j )
    { 
      return dm(i,j);
    }
 
    static const ScalarType & ValueConst( const Teuchos::SerialDenseMatrix<int,ScalarType>& dm, const int i, const int j ) 
    { 
      return dm(i,j);
    }
 
    static void SyncDeviceToHost(Teuchos::SerialDenseMatrix<int,ScalarType> &){ }
 
    static void SyncHostToDevice(Teuchos::SerialDenseMatrix<int,ScalarType> &){ }
 
    
    static void Add( Teuchos::SerialDenseMatrix<int,ScalarType>& thisDM, const Teuchos::SerialDenseMatrix<int,ScalarType>& sourceDM){ 
      thisDM += sourceDM; 
    }
 
    static void PutScalar( Teuchos::SerialDenseMatrix<int,ScalarType> & dm, ScalarType value = Teuchos::ScalarTraits<ScalarType>::zero()){ 
      dm.putScalar(value);
    }
 
    static void Scale( Teuchos::SerialDenseMatrix<int,ScalarType>& dm, ScalarType value){
      dm.scale(value);
    }
 
    //TODO What to do here? Kinda needs random synced version??
    static void Randomize( Teuchos::SerialDenseMatrix<int,ScalarType>& dm){ 
      dm.random();
    }
 
    static void Assign( Teuchos::SerialDenseMatrix<int,ScalarType>& dest, const Teuchos::SerialDenseMatrix<int,ScalarType>& source){ 
      dest.assign(source);
    }
 
    static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType NormFrobenius(const Teuchos::SerialDenseMatrix<int,ScalarType>& dm) {
      return dm.normFrobenius(); 
    }
 
    static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType NormOne(const Teuchos::SerialDenseMatrix<int,ScalarType>& dm) {
      return dm.normOne();
    }
 
    
    static Teuchos::RCP<DenseSolver<ScalarType, Teuchos::SerialDenseMatrix<int,ScalarType>>> 
      createDenseSolver() { 
   
      Teuchos::RCP<DenseSolver<ScalarType, Teuchos::SerialDenseMatrix<int,ScalarType>>> newSolver
        = Teuchos::rcp( new TeuchosDenseSolver<ScalarType>() );
      return newSolver;
    }
 
  };
 
} // namespace Belos
 
#endif // end file BELOS_TEUCHOS_DENSE_MAT_TRAITS_HPP