dune-istl  2.2.0
pardiso.hh
Go to the documentation of this file.
00001 #ifndef DUNE_PARDISO_HH
00002 #define DUNE_PARDISO_HH
00003 
00004 #include <dune/istl/preconditioners.hh>
00005 #include <dune/istl/solvertype.hh>
00006 /* Change this, if your Fortran compiler does not append underscores. */
00007 /* e.g. the AIX compiler:  #define F77_FUNC(func) func                */
00008 
00009 #ifdef AIX
00010 #define F77_FUNC(func)  func 
00011 #else
00012 #define F77_FUNC(func)  func ## _
00013 #endif
00014 
00015 
00016 #ifdef HAVE_PARDISO 
00017 /* PARDISO prototype. */
00018 extern "C" int F77_FUNC(pardisoinit)
00019     (void *, int *, int *);
00020 
00021 extern "C" int F77_FUNC(pardiso)
00022     (void *, int *, int *, int *, int *, int *, 
00023      double *, int *, int *, int *, int *, int *, 
00024      int *, double *, double *, int *);
00025 #endif 
00026 
00027 namespace Dune {
00028 
00029 
00034   template<class M, class X, class Y>
00035   class SeqPardiso : public Preconditioner<X,Y> {
00036   public:
00038     typedef M matrix_type;
00040     typedef X domain_type;
00042     typedef Y range_type;
00044     typedef typename X::field_type field_type;
00045     
00046     typedef typename M::RowIterator RowIterator;
00047     typedef typename M::ColIterator ColIterator;
00048 
00049     // define the category
00050     enum {
00052       category=SolverCategory::sequential
00053     };
00054 
00060     SeqPardiso (const M& A)
00061       : A_(A)
00062     {
00063 #ifdef HAVE_PARDISO
00064         
00065         mtype_ = 11;
00066         nrhs_ = 1;
00067         num_procs_ = 1;
00068         maxfct_ = 1;    
00069         mnum_   = 1;         
00070         msglvl_ = 0;        
00071         error_  = 0;        
00072 
00073         n_ = A_.rowdim();
00074         int nnz = 0;
00075         RowIterator endi = A_.end();
00076         for (RowIterator i = A_.begin(); i != endi; ++i)
00077         {
00078                 if (A_.rowdim(i.index()) != 1)
00079                         DUNE_THROW(NotImplemented, "SeqPardiso: row blocksize != 1.");
00080                 ColIterator endj = (*i).end();
00081                 for (ColIterator j = (*i).begin(); j != endj; ++j) {
00082                         if (A_.coldim(j.index()) != 1)
00083                                 DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
00084                         nnz++;
00085                 }
00086         }
00087                   
00088         std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
00089         
00090         a_ = new double[nnz];
00091         ia_ = new int[n_+1];
00092         ja_ = new int[nnz];
00093         
00094         int count = 0;
00095         for (RowIterator i = A_.begin(); i != endi; ++i)
00096         {
00097                 ia_[i.index()] = count+1;
00098                 ColIterator endj = (*i).end();
00099                 for (ColIterator j = (*i).begin(); j != endj; ++j) {
00100                         a_[count] = *j;
00101                         ja_[count] = j.index()+1;
00102                         
00103                         count++;
00104                 }
00105         }
00106         ia_[n_] = count+1;
00107         
00108         F77_FUNC(pardisoinit) (pt_,  &mtype_, iparm_); 
00109 
00110         int phase = 11;
00111         int idum;
00112         double ddum;
00113         iparm_[2]  = num_procs_;
00114         
00115         F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
00116                        &n_, a_, ia_, ja_, &idum, &nrhs_,
00117                        iparm_, &msglvl_, &ddum, &ddum, &error_);
00118       
00119         if (error_ != 0) 
00120                 DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
00121         
00122         std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
00123         
00124 #else
00125         DUNE_THROW(NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
00126 #endif
00127     }
00128 
00134     virtual void pre (X& x, Y& b) {}
00135 
00141     virtual void apply (X& v, const Y& d)
00142     {
00143 #ifdef HAVE_PARDISO
00144         int phase = 33;
00145 
00146         iparm_[7] = 1;       /* Max numbers of iterative refinement steps. */
00147         int idum;
00148 
00149         double x[n_];
00150         double b[n_];
00151         for (int i = 0; i < n_; i++) {
00152                 x[i] = v[i];
00153                 b[i] = d[i];
00154         }
00155         
00156         F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
00157                        &n_, a_, ia_, ja_, &idum, &nrhs_,
00158                        iparm_, &msglvl_, b, x, &error_);
00159         
00160         if (error_ != 0) 
00161                 DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
00162         
00163         for (int i = 0; i < n_; i++) 
00164                 v[i] = x[i];
00165         
00166         std::cout << "SeqPardiso: Backsolve completed." << std::endl;
00167 #endif
00168     }
00169 
00175     virtual void post (X& x) {}
00176 
00177     ~SeqPardiso() 
00178     {
00179 #ifdef HAVE_PARDISO
00180         int phase = -1;                 // Release internal memory. 
00181         int idum;
00182         double ddum;
00183 
00184         F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
00185                        &n_, &ddum, ia_, ja_, &idum, &nrhs_,
00186                        iparm_, &msglvl_, &ddum, &ddum, &error_);
00187         delete[] a_;
00188         delete[] ia_;
00189         delete[] ja_;
00190 #endif
00191     }
00192   
00193   private:
00194     M A_; 
00195     int n_; 
00196     double *a_; 
00197     int *ia_; 
00198     int *ja_; 
00199     int mtype_; 
00200     int nrhs_; 
00201     void *pt_[64]; 
00202     int iparm_[64]; 
00203     int maxfct_;        
00204     int mnum_;  
00205     int msglvl_;    
00206     int error_;      
00207     int num_procs_; 
00208   };
00209 
00210   template<class M, class X, class Y>
00211   struct IsDirectSolver<SeqPardiso<M,X,Y> >
00212   {
00213     enum{ value=true};
00214   };
00215   
00216   
00217 } // end namespace Dune
00218 #endif
00219