dune-istl  3.0-git
pardiso.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_PARDISO_HH
4 #define DUNE_ISTL_PARDISO_HH
5 
8 
9 #if HAVE_PARDISO
10 // PARDISO prototypes
11 extern "C" void pardisoinit(void *, int *, int *, int *, double *, int *);
12 
13 extern "C" void pardiso(void *, int *, int *, int *, int *, int *,
14  double *, int *, int *, int *, int *, int *,
15  int *, double *, double *, int *, double *);
16 
17 namespace Dune {
18 
19 
24  template<class M, class X, class Y>
25  class SeqPardiso : public Preconditioner<X,Y> {
26  public:
28  typedef M matrix_type;
30  typedef X domain_type;
32  typedef Y range_type;
34  typedef typename X::field_type field_type;
35 
36  typedef typename M::RowIterator RowIterator;
37  typedef typename M::ColIterator ColIterator;
38 
39  // define the category
40  enum {
43  };
44 
50  SeqPardiso (const M& A)
51  : A_(A)
52  {
53  mtype_ = 11;
54  nrhs_ = 1;
55  num_procs_ = 1;
56  maxfct_ = 1;
57  mnum_ = 1;
58  msglvl_ = 0;
59  error_ = 0;
60 
61  n_ = A_.N();
62  int nnz = 0;
63  RowIterator endi = A_.end();
64  for (RowIterator i = A_.begin(); i != endi; ++i)
65  {
66  ColIterator endj = (*i).end();
67  for (ColIterator j = (*i).begin(); j != endj; ++j) {
68  if (j->size() != 1)
69  DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
70  nnz++;
71  }
72  }
73 
74  std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
75 
76  a_ = new double[nnz];
77  ia_ = new int[n_+1];
78  ja_ = new int[nnz];
79 
80  int count = 0;
81  for (RowIterator i = A_.begin(); i != endi; ++i)
82  {
83  ia_[i.index()] = count+1;
84  ColIterator endj = (*i).end();
85  for (ColIterator j = (*i).begin(); j != endj; ++j) {
86  a_[count] = *j;
87  ja_[count] = j.index()+1;
88 
89  count++;
90  }
91  }
92  ia_[n_] = count+1;
93 
94  pardisoinit(pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
95 
96  int phase = 11;
97  int idum;
98  double ddum;
99  iparm_[2] = num_procs_;
100 
101  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
102  &n_, a_, ia_, ja_, &idum, &nrhs_,
103  iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
104 
105  if (error_ != 0)
106  DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
107 
108  std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
109  }
110 
116  virtual void pre (X& x, Y& b) {}
117 
123  virtual void apply (X& v, const Y& d)
124  {
125  int phase = 33;
126 
127  iparm_[7] = 1; /* Max numbers of iterative refinement steps. */
128  int idum;
129 
130  double x[n_];
131  double b[n_];
132  for (int i = 0; i < n_; i++) {
133  x[i] = v[i];
134  b[i] = d[i];
135  }
136 
137  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
138  &n_, a_, ia_, ja_, &idum, &nrhs_,
139  iparm_, &msglvl_, b, x, &error_, dparm_);
140 
141  if (error_ != 0)
142  DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
143 
144  for (int i = 0; i < n_; i++)
145  v[i] = x[i];
146 
147  std::cout << "SeqPardiso: Backsolve completed." << std::endl;
148  }
149 
155  virtual void post (X& x) {}
156 
157  ~SeqPardiso()
158  {
159  int phase = -1; // Release internal memory.
160  int idum;
161  double ddum;
162 
163  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
164  &n_, &ddum, ia_, ja_, &idum, &nrhs_,
165  iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
166  delete[] a_;
167  delete[] ia_;
168  delete[] ja_;
169  }
170 
171  private:
172  M A_;
173  int n_;
174  double *a_;
175  int *ia_;
176  int *ja_;
177  int mtype_;
178  int solver_;
179  int nrhs_;
180  void *pt_[64];
181  int iparm_[64];
182  double dparm_[64];
183  int maxfct_;
184  int mnum_;
185  int msglvl_;
186  int error_;
187  int num_procs_;
188  };
189 
190  template<class M, class X, class Y>
191  struct IsDirectSolver<SeqPardiso<M,X,Y> >
192  {
193  enum { value=true};
194  };
195 } // end namespace Dune
196 
197 #endif //HAVE_PARDISO
198 #endif
Category for sequential solvers.
Definition: solvercategory.hh:21
Templates characterizing the type of a solver.
Define general preconditioner interface.
Definition: basearray.hh:19