3#ifndef DUNE_ISTL_PARDISO_HH
4#define DUNE_ISTL_PARDISO_HH
11extern "C" void pardisoinit(
void *,
int *,
int *,
int *,
double *,
int *);
13extern "C" void pardiso(
void *,
int *,
int *,
int *,
int *,
int *,
14 double *,
int *,
int *,
int *,
int *,
int *,
15 int *,
double *,
double *,
int *,
double *);
24 template<
class M,
class X,
class Y>
25 class SeqPardiso :
public Preconditioner<X,Y> {
28 typedef M matrix_type;
30 typedef X domain_type;
34 typedef typename X::field_type field_type;
36 typedef typename M::RowIterator RowIterator;
37 typedef typename M::ColIterator ColIterator;
51 SeqPardiso (const M& A)
64 RowIterator endi = A_.end();
65 for (RowIterator i = A_.begin(); i != endi; ++i)
67 ColIterator endj = (*i).end();
68 for (ColIterator j = (*i).begin(); j != endj; ++j) {
70 DUNE_THROW(NotImplemented,
"SeqPardiso: column blocksize != 1.");
75 std::cout <<
"dimension = " << n_ <<
", number of nonzeros = " << nnz << std::endl;
82 for (RowIterator i = A_.begin(); i != endi; ++i)
84 ia_[i.index()] = count+1;
85 ColIterator endj = (*i).end();
86 for (ColIterator j = (*i).begin(); j != endj; ++j) {
88 ja_[count] = j.index()+1;
95 pardisoinit(pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
100 iparm_[2] = num_procs_;
102 pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
103 &n_, a_, ia_, ja_, &idum, &nrhs_,
104 iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
107 DUNE_THROW(MathError,
"Constructor SeqPardiso: Factorization failed. Error code " << error_);
109 std::cout <<
"Constructor SeqPardiso: Factorization completed." << std::endl;
117 virtual void pre (X& x, Y& b) {}
124 virtual void apply (X& v,
const Y& d)
133 for (
int i = 0; i < n_; i++) {
138 pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
139 &n_, a_, ia_, ja_, &idum, &nrhs_,
140 iparm_, &msglvl_, b, x, &error_, dparm_);
143 DUNE_THROW(MathError,
"SeqPardiso.apply: Backsolve failed. Error code " << error_);
145 for (
int i = 0; i < n_; i++)
148 std::cout <<
"SeqPardiso: Backsolve completed." << std::endl;
156 virtual void post (X& x) {}
164 pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
165 &n_, &ddum, ia_, ja_, &idum, &nrhs_,
166 iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
191 template<
class M,
class X,
class Y>
192 struct IsDirectSolver<SeqPardiso<M,X,Y> >
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