pardiso.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_PARDISO_HH
#define DUNE_ISTL_PARDISO_HH
 
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvertype.hh>
 
#if HAVE_PARDISO
// PARDISO prototypes
extern "C" void pardisoinit(void *, int *, int *, int *, double *, int *);
 
extern "C" void pardiso(void *, int *, int *, int *, int *, int *,
                        double *, int *, int *, int *, int *, int *,
                        int *, double *, double *, int *, double *);
 
namespace Dune {
 
 
  template<class M, class X, class Y>
  class SeqPardiso : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
 
    typedef typename M::RowIterator RowIterator;
    typedef typename M::ColIterator ColIterator;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqPardiso (const M& A)
      : A_(A)
    {
      mtype_ = 11;
      nrhs_ = 1;
      num_procs_ = 1;
      maxfct_ = 1;
      mnum_   = 1;
      msglvl_ = 0;
      error_  = 0;
 
      n_ = A_.N();
      int nnz = 0;
      RowIterator endi = A_.end();
      for (RowIterator i = A_.begin(); i != endi; ++i)
      {
        ColIterator endj = (*i).end();
        for (ColIterator j = (*i).begin(); j != endj; ++j) {
          if (j->size() != 1)
            DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
          nnz++;
        }
      }
 
      std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
 
      a_ = new double[nnz];
      ia_ = new int[n_+1];
      ja_ = new int[nnz];
 
      int count = 0;
      for (RowIterator i = A_.begin(); i != endi; ++i)
      {
        ia_[i.index()] = count+1;
        ColIterator endj = (*i).end();
        for (ColIterator j = (*i).begin(); j != endj; ++j) {
          a_[count] = *j;
          ja_[count] = j.index()+1;
 
          count++;
        }
      }
      ia_[n_] = count+1;
 
      pardisoinit(pt_,  &mtype_, &solver_, iparm_, dparm_, &error_);
 
      int phase = 11;
      int idum;
      double ddum;
      iparm_[2]  = num_procs_;
 
      pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
              &n_, a_, ia_, ja_, &idum, &nrhs_,
              iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
 
      if (error_ != 0)
        DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
 
      std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
    }
 
    virtual void pre (X& x, Y& b) {}
 
    virtual void apply (X& v, const Y& d)
    {
      int phase = 33;
 
      iparm_[7] = 1;         /* Max numbers of iterative refinement steps. */
      int idum;
 
      double x[n_];
      double b[n_];
      for (int i = 0; i < n_; i++) {
        x[i] = v[i];
        b[i] = d[i];
      }
 
      pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
              &n_, a_, ia_, ja_, &idum, &nrhs_,
              iparm_, &msglvl_, b, x, &error_, dparm_);
 
      if (error_ != 0)
        DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
 
      for (int i = 0; i < n_; i++)
        v[i] = x[i];
 
      std::cout << "SeqPardiso: Backsolve completed." << std::endl;
    }
 
    virtual void post (X& x) {}
 
    ~SeqPardiso()
    {
      int phase = -1;                   // Release internal memory.
      int idum;
      double ddum;
 
      pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
              &n_, &ddum, ia_, ja_, &idum, &nrhs_,
              iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
      delete[] a_;
      delete[] ia_;
      delete[] ja_;
    }
 
  private:
    M A_; 
    int n_; 
    double *a_; 
    int *ia_; 
    int *ja_; 
    int mtype_; 
    int solver_; 
    int nrhs_; 
    void *pt_[64]; 
    int iparm_[64]; 
    double dparm_[64]; 
    int maxfct_;        
    int mnum_;  
    int msglvl_;    
    int error_;      
    int num_procs_; 
  };
 
  template<class M, class X, class Y>
  struct IsDirectSolver<SeqPardiso<M,X,Y> >
  {
    enum { value=true};
  };
} // end namespace Dune
 
#endif //HAVE_PARDISO
#endif