PETScMatrix.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _PETSC_MATRIX_H
#define _PETSC_MATRIX_H
 
#include "DomainDecomposition.h"
#include "SystemMatrix.h"
#include "PETScSupport.h"
#include "PETScSolParams.h"
 
#include <array>
#include <memory>
 
using PetscIntVec = std::vector<PetscInt>;    
using PetscIntMat = std::vector<PetscIntVec>; 
using PetscRealVec = std::vector<PetscReal>;  
using ISVec = std::vector<IS>;                
using ISMat = std::vector<ISVec>;             
 
 
class PETScVector : public StdVector
{
public:
  explicit PETScVector(const ProcessAdm& padm);
  PETScVector(const ProcessAdm& padm, size_t n);
  PETScVector(const ProcessAdm& padm, const Real* values, size_t n);
  PETScVector(const PETScVector& vec);
  ~PETScVector() override;
 
  LinAlg::MatrixType getType() const override { return LinAlg::PETSC; }
 
  SystemVector* copy() const override { return new PETScVector(*this); }
 
  void init(Real value) override;
 
  void redim(size_t n) override;
 
  bool endAssembly() override;
 
  Real L1norm() const override;
 
  Real L2norm() const override;
 
  Real Linfnorm() const override;
 
  Vec& getVector() { return x; }
  const Vec& getVector() const { return x; }
 
  const ProcessAdm& getAdm() const { return adm; }
 
protected:
  Vec x;                  
  const ProcessAdm& adm;  
};
 
 
class PETScMatrix;
 
class PETScVectors : public SystemVector
{
public:
  PETScVectors(const PETScMatrix& A, int nvec);
 
  ~PETScVectors();
 
  LinAlg::MatrixType getType() const override { return LinAlg::PETSC; }
 
  SystemVector* copy() const override { return nullptr; }
 
  size_t dim() const override { return myDim; }
 
  void redim(size_t) override {}
 
  Real* getPtr() override { return nullptr; }
  const Real* getRef() const override { return nullptr; }
  const Vector& vec() const override { static Vector empty; return empty; }
 
  void init(Real) override {}
 
  void mult(Real) override {}
 
  void add(const SystemVector&, Real) override {}
 
  Real L1norm() const override { return 0.0; }
 
  Real L2norm() const override { return 0.0; }
 
  Real Linfnorm() const override { return 0.0; }
 
  void assemble(const Vectors& vecs, const IntVec& meqn, int = 0) override;
 
  Vec get(size_t idx) { return vectors[idx]; }
 
  size_t size() const { return vectors.size(); }
 
private:
  size_t myDim; 
  std::vector<Vec> vectors; 
  const PETScMatrix& myA; 
};
 
 
class PETScMatrix : public SystemMatrix
{
public:
  PETScMatrix(const ProcessAdm& padm, const LinSolParams& spar);
  ~PETScMatrix() override;
 
  LinAlg::MatrixType getType() const override { return LinAlg::PETSC; }
 
  SystemMatrix* copy() const override;
 
  size_t dim(int idim = 1) const override;
 
  bool init(int maxEq,
            const IntMat* elms = nullptr,
            const IntMat* neighs = nullptr,
            const IntVec* part = nullptr);
 
  void preAssemble(const std::vector<IntVec>& MMNPC, size_t nel) override;
 
  void initAssembly(const SAM& sam, char) override;
 
  void init() override;
 
  bool assemble(const Matrix& eM, const SAM& sam, int e) override;
  bool assemble(const Matrix& eM, const SAM& sam, SystemVector& B, int e) override;
 
  bool assemble(const Matrix& eM, const SAM& sam,
                SystemVector& B, const IntVec& meq) override;
 
  bool assemble(const Matrix& eM, const IntVec& meq) override;
 
  bool endAssembly() override;
 
  bool multiply(const SystemVector& B, SystemVector& C) const override;
 
  bool solve(SystemVector& B, Real*) override;
 
  bool solve(const SystemVector& B, SystemVector& x) override;
 
  bool solve(SystemVector& B, Matrix& sol) override;
 
  void mult(Real alpha) override;
 
  bool solveEig(PETScMatrix& B, RealArray& eigVal, Matrix& eigVec,
                int nev, Real shift = Real(0), int iop = 1);
 
  Real Linfnorm() const override;
 
  Mat& getMatrix() { return pA; }
  const Mat& getMatrix() const { return pA; }
 
  const std::vector<Mat>& getBlockMatrices() const { return matvec; }
 
  const std::vector<IS>& getIS() const { return isvec; }
 
  bool setParameters(bool setup);
 
  const ProcessAdm& getAdm() const { return adm; }
 
  const DomainDecomposition& getDD() const;
 
  bool add(Real sigma, int ieq) override;
 
protected:
  bool solve(const Vec& b, Vec& x, bool knoll);
 
  PETScMatrix(const PETScMatrix& A) = delete;
 
  PETScMatrix(const ProcessAdm& padm,
              const PETScSolParams& spar);
 
  void setupSparsity(const IntVec& elms, const SAM& sam);
 
  void setupBlockSparsity(const IntVec& elms, const SAM& sam);
 
  void setupGlb2Blk(const SAM& sam);
 
  Mat preAllocator(const int nrows, const int ncols = 0) const;
 
  std::vector<Mat> preAllocators() const;
 
  Mat                 pA;              
  PetscInt            nrow;            
  PetscInt            ncol;            
  KSP                 ksp;             
  MatNullSpace*       nsp;             
  const ProcessAdm&   adm;             
  PETScSolParams      solParams;       
  bool                setParams;       
  std::string         forcedKSPType;   
  PetscInt            ISsize;          
  PetscRealVec        coords;          
  ISMat               dirIndexSet;     
  int                 nLinSolves;      
  bool                assembled;       
  bool                factored;        
 
  IS glob2LocEq = nullptr; 
  std::vector<Mat> matvec; 
 
  std::vector<IS> isvec; 
 
  std::vector<std::array<int,2>> glb2Blk; 
  std::unique_ptr<DomainDecomposition> m_dd{}; 
};
 
 
PETScVector operator*(const SystemMatrix& A, const PETScVector& b);
 
PETScVector operator/(SystemMatrix& A, const PETScVector& b);
 
#endif