Dune Core Modules (2.6.0)

pardiso.hh
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
11extern "C" void pardisoinit(void *, int *, int *, int *, double *, int *);
12
13extern "C" void pardiso(void *, int *, int *, int *, int *, int *,
14 double *, int *, int *, int *, int *, int *,
15 int *, double *, double *, int *, double *);
16
17namespace 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
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:58
#define DUNE_DEPRECATED_MSG(text)
Mark some entity as deprecated.
Definition: deprecated.hh:169
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:10
Define general preconditioner interface.
Templates characterizing the type of a solver.
Category
Definition: solvercategory.hh:21
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:23
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)