This file is indexed.

/usr/include/dune/istl/pardiso.hh is in libdune-istl-dev 2.4.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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
// -*- 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>
/* Change this, if your Fortran compiler does not append underscores. */
/* e.g. the AIX compiler:  #define F77_FUNC(func) func                */

#ifdef AIX
#define F77_FUNC(func)  func
#else
#define F77_FUNC(func)  func ## _
#endif


#ifdef HAVE_PARDISO
/* PARDISO prototype. */
extern "C" int F77_FUNC(pardisoinit)
  (void *, int *, int *);

extern "C" int F77_FUNC(pardiso)
  (void *, int *, int *, int *, int *, int *,
  double *, int *, int *, int *, int *, int *,
  int *, double *, double *, int *);
#endif

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)
    {
#ifdef HAVE_PARDISO

      mtype_ = 11;
      nrhs_ = 1;
      num_procs_ = 1;
      maxfct_ = 1;
      mnum_   = 1;
      msglvl_ = 0;
      error_  = 0;

      n_ = A_.rowdim();
      int nnz = 0;
      RowIterator endi = A_.end();
      for (RowIterator i = A_.begin(); i != endi; ++i)
      {
        if (A_.rowdim(i.index()) != 1)
          DUNE_THROW(NotImplemented, "SeqPardiso: row blocksize != 1.");
        ColIterator endj = (*i).end();
        for (ColIterator j = (*i).begin(); j != endj; ++j) {
          if (A_.coldim(j.index()) != 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;

      F77_FUNC(pardisoinit) (pt_,  &mtype_, iparm_);

      int phase = 11;
      int idum;
      double ddum;
      iparm_[2]  = num_procs_;

      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
                         &n_, a_, ia_, ja_, &idum, &nrhs_,
                         iparm_, &msglvl_, &ddum, &ddum, &error_);

      if (error_ != 0)
        DUNE_THROW(MathError, "Constructor SeqPardiso: Factorization failed. Error code " << error_);

      std::cout << "Constructor SeqPardiso: Factorization completed." << std::endl;

#else
      DUNE_THROW(NotImplemented, "no Pardiso library available, reconfigure with correct --with-pardiso options");
#endif
    }

    /*!
       \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)
    {
#ifdef HAVE_PARDISO
      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];
      }

      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
                         &n_, a_, ia_, ja_, &idum, &nrhs_,
                         iparm_, &msglvl_, b, x, &error_);

      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;
#endif
    }

    /*!
       \brief Clean up.

       \copydoc Preconditioner::post(X&)
     */
    virtual void post (X& x) {}

    ~SeqPardiso()
    {
#ifdef HAVE_PARDISO
      int phase = -1;                   // Release internal memory.
      int idum;
      double ddum;

      F77_FUNC(pardiso) (pt_, &maxfct_, &mnum_, &mtype_, &phase,
                         &n_, &ddum, ia_, ja_, &idum, &nrhs_,
                         iparm_, &msglvl_, &ddum, &ddum, &error_);
      delete[] a_;
      delete[] ia_;
      delete[] ja_;
#endif
    }

  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 nrhs_; //!< number of right hand sides
    void *pt_[64]; //!< internal solver memory pointer
    int iparm_[64]; //!< Pardiso 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