operator/PETSc.hh

Go to the documentation of this file.

/* PETSc - use solvers from the PETSc library 
 * http://www-fp.mcs.anl.gov/petsc/
 */

#ifndef op_PETSC_hh
#define op_PETSC_hh

extern "C" {
  // forwards for PETSc
  typedef struct _p_SLES* SLES;
  typedef struct _p_Mat* Mat;
  typedef struct _p_Vec* Vec;
}
#include <string>

#include "compositions.hh"
#include "basics/typedefs.hh"
#include "matrix.hh"

namespace concepts {

  // forward declarations
  template<class F>
  class Function;

  template<class F>
  class Vector;

  template<class F>
  class Space;

  template<class F, class G>
  class BilinearForm;

  // **************************************************************** PETSc **

  class PETSc : public Operator<Real> {
  public:
    PETSc(Operator<Real>& A, const Real maxeps, const std::string ksptype, 
          const std::string pctype = std::string("none"));
    ~PETSc();
    virtual void operator()(const Function<Real>& fncY, Function<Real>& fncX);

    uint iterations() const { return iterations_; }
  protected:
    virtual std::ostream& info(std::ostream& os) const;
  private:
    typedef int (*multPointer)(Mat, Vec, Vec);
    typedef void (*PETScmultPointer)();

    Operator<Real>& A_;

    int iterations_;

    Mat PETScA_;
    bool destroyPETScA_;

    SLES solver_;

    Vec x_, y_;
  };

  // ************************************************************* PETScMat **

  class PETScMat : public Matrix<Real> {
  public:
    PETScMat(const Space<Real>& spcX, const Space<Real>& spcY);

    PETScMat(const Space<Real>& spc, BilinearForm<Real,Real>& bf);

    PETScMat(const SparseMatrix<Real>& matrix);

    ~PETScMat();

    virtual void operator()(const Function<double>& fncY,
                            Function<double>& fncX);
    virtual void operator()(const Function<std::complex<double> >& fncY,
                            Function<std::complex<double> >& fncX);
    void operator()(const Vector<double>& fncY,
                    Vector<double>& fncX);
    virtual void transpMult(const Vector<double>& fncY,
                            Vector<double>& fncX);
    virtual void transpMult(const Vector<std::complex<double> >& fncY,
                            Vector<std::complex<double> >& fncX);

    virtual Real operator()(const uint i, const uint j) const;
    virtual Real& operator()(const uint i, const uint j);

    virtual const Space<double>& spaceX() const { return spcX_; }
    virtual const Space<double>& spaceY() const { return spcY_; }

    operator Mat() { return A_; }

    void storeMatlab(const char* name) const;
  protected:
    std::ostream& info(std::ostream& os) const;
  private:
    const Space<Real>& spcX_;

    uint nX_;

    const Space<Real>& spcY_;

    uint nY_;

    Vec x_, y_;

    Mat A_;
  };

} // namespace concepts

#endif // op_PETSC_hh