SystemMatrix.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _SYSTEM_MATRIX_H
#define _SYSTEM_MATRIX_H
 
#include "MatVec.h"
#include "LinAlgenums.h"
 
class SAM;
class LinSolParams;
class ProcessAdm;
 
using IntVec = std::vector<int>; 
 
 
class SystemVector
{
public:
  static SystemVector* create(const ProcessAdm* adm, LinAlg::MatrixType vtype);
 
protected:
  SystemVector() {}
 
public:
  virtual ~SystemVector() {}
 
  virtual LinAlg::MatrixType getType() const = 0;
 
  virtual SystemVector* copy() const = 0;
 
  virtual size_t dim() const = 0;
 
  virtual void redim(size_t n) = 0;
 
  virtual void resize(size_t n, bool = false) { this->redim(n); }
 
  bool empty() const { return this->dim() == 0; }
 
  virtual Real* getPtr() = 0;
  virtual const Real* getRef() const = 0;
  virtual const Vector& vec() const = 0;
 
  virtual void init(Real value = Real(0)) = 0;
 
  SystemVector& copy(const SystemVector& x);
 
  virtual void assemble(const Vectors& vecs,
                        const IntVec& meqn, int neq) = 0;
 
  virtual bool endAssembly() { return true; }
 
  virtual void mult(Real alpha) = 0;
 
  virtual void add(const SystemVector& vec, Real scale = Real(1)) = 0;
 
  virtual Real L1norm() const = 0;
 
  virtual Real L2norm() const = 0;
 
  virtual Real Linfnorm() const = 0;
 
  virtual void dump(std::ostream&, LinAlg::StorageFormat,
                    const char* = nullptr) const {}
 
protected:
  virtual std::ostream& write(std::ostream& os) const { return os; }
 
  friend std::ostream& operator<<(std::ostream& os, const SystemVector& X)
  {
    return X.write(os);
  }
};
 
 
class StdVector : public SystemVector, public Vector
{
public:
  explicit StdVector(size_t n = 0) : Vector(n) {}
  StdVector(const Real* values, size_t n) : Vector(values,n) {}
  explicit StdVector(const std::vector<Real>& vec) : Vector(vec) {}
 
  virtual LinAlg::MatrixType getType() const { return LinAlg::DENSE; }
 
  virtual SystemVector* copy() const { return new StdVector(*this); }
 
  virtual void assemble(const Vectors& vecs,
                        const IntVec& meqn,
                        int neq);
 
  virtual size_t dim() const { return this->size(); }
 
  virtual void redim(size_t n)
  {
    this->Vector::resize(n,utl::RETAIN);
  }
 
  virtual void resize(size_t n, bool forceClear = false)
  {
    this->Vector::resize(n,forceClear);
  }
 
  virtual Real* getPtr() { return this->ptr(); }
  virtual const Real* getRef() const { return this->ptr(); }
  virtual const Vector& vec() const { return *this; }
 
  virtual void init(Real value = Real(0)) { this->fill(value); }
 
  virtual void mult(Real alpha) { this->operator*=(alpha); }
 
  virtual void add(const SystemVector& vec, Real scale)
  {
    this->Vector::add(static_cast<const StdVector&>(vec),scale);
  }
 
  virtual Real L1norm() const { return this->asum(); }
 
  virtual Real L2norm() const { return this->norm2(); }
 
  virtual Real Linfnorm() const { size_t off = 0; return this->normInf(off); }
 
  virtual void dump(std::ostream& os, LinAlg::StorageFormat format,
                    const char* label) const { dump(*this,label,format,os); }
 
  static void dump(const Vector& x, const char* label,
                   LinAlg::StorageFormat format, std::ostream& os);
 
protected:
  virtual std::ostream& write(std::ostream& os) const
  {
    return os << static_cast<const Vector&>(*this);
  }
};
 
 
class SystemMatrix
{
public:
  static SystemMatrix* create(const ProcessAdm* adm, LinAlg::MatrixType mType,
                              int num_thread_SLU = 1);
  static SystemMatrix* create(const ProcessAdm* adm, LinAlg::MatrixType mType,
                              const LinSolParams& spar);
 
protected:
  SystemMatrix() : nonZeroEqs({false}) {}
  SystemMatrix(const SystemMatrix& a) : nonZeroEqs(a.nonZeroEqs) {}
 
public:
  virtual ~SystemMatrix() {}
 
  virtual LinAlg::MatrixType getType() const = 0;
 
  virtual SystemMatrix* copy() const = 0;
 
  virtual bool lockPattern(bool) { return false; }
 
  virtual bool empty() const { return this->dim(0) == 0; }
  bool isZero() const;
 
  virtual size_t dim(int idim = 1) const = 0;
 
  virtual void initAssembly(const SAM& sam, char = 0) = 0;
 
  virtual void preAssemble(const std::vector<IntVec>&, size_t = 0) {}
 
  virtual void compressPattern() {}
 
  virtual void init() = 0;
 
  void initNonZeroEqs();
  bool flagNonZeroEq(int ieq = 0);
  bool flagNonZeroEqs(const SAM& sam, const IntVec& meq);
  bool flagNonZeroEqs(const SystemMatrix& B);
 
  virtual bool endAssembly();
 
  virtual bool assemble(const Matrix& eM, const SAM& sam, int e) = 0;
  virtual bool assemble(const Matrix& eM, const SAM& sam,
                        SystemVector& B, int e) = 0;
  virtual bool assemble(const Matrix& eM, const SAM& sam,
                        SystemVector& B, const IntVec& meq) = 0;
 
  virtual bool assemble(const Matrix& eM, const IntVec& meq) = 0;
 
  virtual bool augment(const SystemMatrix&, size_t, size_t) { return false; }
 
  virtual bool truncate(Real = Real(1.0e-16)) { return false; }
 
  virtual void mult(Real alpha) = 0;
 
  virtual bool add(const SystemMatrix&, Real = Real(1)) { return false; }
 
  virtual bool add(Real, int = 0) { return false; }
 
  virtual bool multiply(const SystemVector&, SystemVector&) const
  { return false; }
 
  virtual bool solve(SystemVector& b, Real* rc = nullptr) = 0;
 
  virtual bool solve(const SystemVector& b, SystemVector& x)
  {
    return this->solve(x.copy(b));
  }
 
  virtual bool solve(SystemVector& b, Matrix& x);
 
  virtual Real Linfnorm() const = 0;
 
  virtual void dump(std::ostream&, LinAlg::StorageFormat,
                    const char* = nullptr) {}
 
  virtual void dump(const char* fileName, std::streamsize precision = 0,
                    LinAlg::StorageFormat format = LinAlg::FLAT);
  virtual bool load(const char*, bool = false) { return false; }
 
  StdVector operator*(const SystemVector& b) const;
 
  StdVector operator/(const SystemVector& b);
 
protected:
  virtual std::ostream& write(std::ostream& os) const { return os; }
 
  friend std::ostream& operator<<(std::ostream& os, const SystemMatrix& A)
  {
    return A.write(os);
  }
 
private:
  std::vector<bool> nonZeroEqs; 
};
 
#endif