dune-istl  2.6-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 
40  virtual SolverCategory::Category category() const
41  {
43  }
44 
50  DUNE_DEPRECATED_MSG("Pardiso isn't supported anymore. Please use SuperLU or UMFPack.")
51  SeqPardiso (const M& A)
52  : A_(A)
53  {
54  mtype_ = 11;
55  nrhs_ = 1;
56  num_procs_ = 1;
57  maxfct_ = 1;
58  mnum_ = 1;
59  msglvl_ = 0;
60  error_ = 0;
61 
62  n_ = A_.N();
63  int nnz = 0;
64  RowIterator endi = A_.end();
65  for (RowIterator i = A_.begin(); i != endi; ++i)
66  {
67  ColIterator endj = (*i).end();
68  for (ColIterator j = (*i).begin(); j != endj; ++j) {
69  if (j->size() != 1)
70  DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
71  nnz++;
72  }
73  }
74 
75  std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
76 
77  a_ = new double[nnz];
78  ia_ = new int[n_+1];
79  ja_ = new int[nnz];
80 
81  int count = 0;
82  for (RowIterator i = A_.begin(); i != endi; ++i)
83  {
84  ia_[i.index()] = count+1;
85  ColIterator endj = (*i).end();
86  for (ColIterator j = (*i).begin(); j != endj; ++j) {
87  a_[count] = *j;
88  ja_[count] = j.index()+1;
89 
90  count++;
91  }
92  }
93  ia_[n_] = count+1;
94 
95  pardisoinit(pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
96 
97  int phase = 11;
98  int idum;
99  double ddum;
100  iparm_[2] = num_procs_;
101 
102  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
103  &n_, a_, ia_, ja_, &idum, &nrhs_,
104  iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
105 
106  if (error_ != 0)
107  DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
108 
109  std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
110  }
111 
117  virtual void pre (X& x, Y& b) {}
118 
124  virtual void apply (X& v, const Y& d)
125  {
126  int phase = 33;
127 
128  iparm_[7] = 1; /* Max numbers of iterative refinement steps. */
129  int idum;
130 
131  double x[n_];
132  double b[n_];
133  for (int i = 0; i < n_; i++) {
134  x[i] = v[i];
135  b[i] = d[i];
136  }
137 
138  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
139  &n_, a_, ia_, ja_, &idum, &nrhs_,
140  iparm_, &msglvl_, b, x, &error_, dparm_);
141 
142  if (error_ != 0)
143  DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
144 
145  for (int i = 0; i < n_; i++)
146  v[i] = x[i];
147 
148  std::cout << "SeqPardiso: Backsolve completed." << std::endl;
149  }
150 
156  virtual void post (X& x) {}
157 
158  ~SeqPardiso()
159  {
160  int phase = -1; // Release internal memory.
161  int idum;
162  double ddum;
163 
164  pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
165  &n_, &ddum, ia_, ja_, &idum, &nrhs_,
166  iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
167  delete[] a_;
168  delete[] ia_;
169  delete[] ja_;
170  }
171 
172  private:
173  M A_;
174  int n_;
175  double *a_;
176  int *ia_;
177  int *ja_;
178  int mtype_;
179  int solver_;
180  int nrhs_;
181  void *pt_[64];
182  int iparm_[64];
183  double dparm_[64];
184  int maxfct_;
185  int mnum_;
186  int msglvl_;
187  int error_;
188  int num_procs_;
189  };
190 
191  template<class M, class X, class Y>
192  struct IsDirectSolver<SeqPardiso<M,X,Y> >
193  {
194  enum { value=true};
195  };
196 } // end namespace Dune
197 
198 #endif //HAVE_PARDISO
199 #endif
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:21
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Dune::IsDirectSolver::value
@ value
Whether this is a direct solver.
Definition: solvertype.hh:22
preconditioners.hh
Define general preconditioner interface.
Dune
Definition: allocator.hh:7
solvertype.hh
Templates characterizing the type of a solver.