/usr/include/dune/istl/pardiso.hh is in libdune-istl-dev 2.5.1-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 | // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_PARDISO_HH
#define DUNE_ISTL_PARDISO_HH
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvertype.hh>
#if HAVE_PARDISO
// PARDISO prototypes
extern "C" void pardisoinit(void *, int *, int *, int *, double *, int *);
extern "C" void pardiso(void *, int *, int *, int *, int *, int *,
double *, int *, int *, int *, int *, int *,
int *, double *, double *, int *, double *);
namespace Dune {
/*! \brief The sequential Pardiso preconditioner.
Put the Pardiso direct solver into the preconditioner framework.
*/
template<class M, class X, class Y>
class SeqPardiso : public Preconditioner<X,Y> {
public:
//! \brief The matrix type the preconditioner is for.
typedef M matrix_type;
//! \brief The domain type of the preconditioner.
typedef X domain_type;
//! \brief The range type of the preconditioner.
typedef Y range_type;
//! \brief The field type of the preconditioner.
typedef typename X::field_type field_type;
typedef typename M::RowIterator RowIterator;
typedef typename M::ColIterator ColIterator;
// define the category
enum {
//! \brief The category the preconditioner is part of
category=SolverCategory::sequential
};
/*! \brief Constructor.
Constructor gets all parameters to operate the prec.
\param A The matrix to operate on.
*/
SeqPardiso (const M& A)
: A_(A)
{
mtype_ = 11;
nrhs_ = 1;
num_procs_ = 1;
maxfct_ = 1;
mnum_ = 1;
msglvl_ = 0;
error_ = 0;
n_ = A_.N();
int nnz = 0;
RowIterator endi = A_.end();
for (RowIterator i = A_.begin(); i != endi; ++i)
{
ColIterator endj = (*i).end();
for (ColIterator j = (*i).begin(); j != endj; ++j) {
if (j->size() != 1)
DUNE_THROW(NotImplemented, "SeqPardiso: column blocksize != 1.");
nnz++;
}
}
std::cout << "dimension = " << n_ << ", number of nonzeros = " << nnz << std::endl;
a_ = new double[nnz];
ia_ = new int[n_+1];
ja_ = new int[nnz];
int count = 0;
for (RowIterator i = A_.begin(); i != endi; ++i)
{
ia_[i.index()] = count+1;
ColIterator endj = (*i).end();
for (ColIterator j = (*i).begin(); j != endj; ++j) {
a_[count] = *j;
ja_[count] = j.index()+1;
count++;
}
}
ia_[n_] = count+1;
pardisoinit(pt_, &mtype_, &solver_, iparm_, dparm_, &error_);
int phase = 11;
int idum;
double ddum;
iparm_[2] = num_procs_;
pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
&n_, a_, ia_, ja_, &idum, &nrhs_,
iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
if (error_ != 0)
DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);
std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;
}
/*!
\brief Prepare the preconditioner.
\copydoc Preconditioner::pre(X&,Y&)
*/
virtual void pre (X& x, Y& b) {}
/*!
\brief Apply the preconditioner.
\copydoc Preconditioner::apply(X&,const Y&)
*/
virtual void apply (X& v, const Y& d)
{
int phase = 33;
iparm_[7] = 1; /* Max numbers of iterative refinement steps. */
int idum;
double x[n_];
double b[n_];
for (int i = 0; i < n_; i++) {
x[i] = v[i];
b[i] = d[i];
}
pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
&n_, a_, ia_, ja_, &idum, &nrhs_,
iparm_, &msglvl_, b, x, &error_, dparm_);
if (error_ != 0)
DUNE_THROW(MathError, "SeqPardiso.apply: Backsolve failed. Error code " << error_);
for (int i = 0; i < n_; i++)
v[i] = x[i];
std::cout << "SeqPardiso: Backsolve completed." << std::endl;
}
/*!
\brief Clean up.
\copydoc Preconditioner::post(X&)
*/
virtual void post (X& x) {}
~SeqPardiso()
{
int phase = -1; // Release internal memory.
int idum;
double ddum;
pardiso(pt_, &maxfct_, &mnum_, &mtype_, &phase,
&n_, &ddum, ia_, ja_, &idum, &nrhs_,
iparm_, &msglvl_, &ddum, &ddum, &error_, dparm_);
delete[] a_;
delete[] ia_;
delete[] ja_;
}
private:
M A_; //!< The matrix we operate on.
int n_; //!< dimension of the system
double *a_; //!< matrix values
int *ia_; //!< indices to rows
int *ja_; //!< column indices
int mtype_; //!< matrix type, currently only 11 (real unsymmetric matrix) is supported
int solver_; //!< solver method
int nrhs_; //!< number of right hand sides
void *pt_[64]; //!< internal solver memory pointer
int iparm_[64]; //!< Pardiso integer control parameters
double dparm_[64]; //!< Pardiso double control parameters
int maxfct_; //!< Maximum number of numerical factorizations
int mnum_; //!< Which factorization to use
int msglvl_; //!< flag to print statistical information
int error_; //!< error flag
int num_procs_; //!< number of processors
};
template<class M, class X, class Y>
struct IsDirectSolver<SeqPardiso<M,X,Y> >
{
enum { value=true};
};
} // end namespace Dune
#endif //HAVE_PARDISO
#endif
|