This file is indexed.

/usr/include/dune/istl/ilu.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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_ILU_HH
#define DUNE_ISTL_ILU_HH

#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip>
#include <string>
#include <set>
#include <map>

#include <dune/common/fmatrix.hh>
#include "istlexception.hh"

/** \file
 * \brief  ???
 */

namespace Dune {

  /** @addtogroup ISTL_Kernel
          @{
   */

  class MatrixBlockError : public virtual Dune::FMatrixError {
  public:
    int r, c;
  };


  //! compute ILU decomposition of A. A is overwritten by its decomposition
  template<class M>
  void bilu0_decomposition (M& A)
  {
    // iterator types
    typedef typename M::RowIterator rowiterator;
    typedef typename M::ColIterator coliterator;
    typedef typename M::block_type block;

    // implement left looking variant with stored inverse
    rowiterator endi=A.end();
    for (rowiterator i=A.begin(); i!=endi; ++i)
    {
      // coliterator is diagonal after the following loop
      coliterator endij=(*i).end();           // end of row i
      coliterator ij;

      // eliminate entries left of diagonal; store L factor
      for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
      {
        // find A_jj which eliminates A_ij
        coliterator jj = A[ij.index()].find(ij.index());

        // compute L_ij = A_jj^-1 * A_ij
        (*ij).rightmultiply(*jj);

        // modify row
        coliterator endjk=A[ij.index()].end();                 // end of row j
        coliterator jk=jj; ++jk;
        coliterator ik=ij; ++ik;
        while (ik!=endij && jk!=endjk)
          if (ik.index()==jk.index())
          {
            block B(*jk);
            B.leftmultiply(*ij);
            *ik -= B;
            ++ik; ++jk;
          }
          else
          {
            if (ik.index()<jk.index())
              ++ik;
            else
              ++jk;
          }
      }

      // invert pivot and store it in A
      if (ij.index()!=i.index())
        DUNE_THROW(ISTLError,"diagonal entry missing");
      try {
        (*ij).invert();   // compute inverse of diagonal block
      }
      catch (Dune::FMatrixError & e) {
        DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
                   << i.index() << "][" << ij.index() << "]" << e.what();
                   th__ex.r=i.index(); th__ex.c=ij.index(););
      }
    }
  }

  //! LU backsolve with stored inverse
  template<class M, class X, class Y>
  void bilu_backsolve (const M& A, X& v, const Y& d)
  {
    // iterator types
    typedef typename M::ConstRowIterator rowiterator;
    typedef typename M::ConstColIterator coliterator;
    typedef typename Y::block_type dblock;
    typedef typename X::block_type vblock;

    // lower triangular solve
    rowiterator endi=A.end();
    for (rowiterator i=A.begin(); i!=endi; ++i)
    {
      dblock rhs(d[i.index()]);
      for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
        (*j).mmv(v[j.index()],rhs);
      v[i.index()] = rhs;           // Lii = I
    }

    // upper triangular solve
    rowiterator rendi=A.beforeBegin();
    for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
    {
      vblock rhs(v[i.index()]);
      coliterator j;
      for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
        (*j).mmv(v[j.index()],rhs);
      v[i.index()] = 0;
      (*j).umv(rhs,v[i.index()]);           // diagonal stores inverse!
    }
  }



  // recursive function template to access first entry of a matrix
  template<class M>
  typename M::field_type& firstmatrixelement (M& A)
  {
    return firstmatrixelement(*(A.begin()->begin()));
  }

  template<class K, int n, int m>
  K& firstmatrixelement (FieldMatrix<K,n,m>& A)
  {
    return A[0][0];
  }

  template<class K>
  K& firstmatrixelement (FieldMatrix<K,1,1>& A)
  {
    return A[0][0];
  }


  /*! ILU decomposition of order n
          Computes ILU decomposition of order n. The matrix ILU should
      be an empty matrix in row_wise creation mode. This allows the user
      to either specify the number of nonzero elements or to
          determine it automatically at run-time.
   */
  template<class M>
  void bilu_decomposition (const M& A, int n, M& ILU)
  {
    // iterator types
    typedef typename M::ColIterator coliterator;
    typedef typename M::ConstRowIterator crowiterator;
    typedef typename M::ConstColIterator ccoliterator;
    typedef typename M::CreateIterator createiterator;
    typedef typename M::field_type K;
    typedef std::map<size_t, int> map;
    typedef typename map::iterator mapiterator;

    // symbolic factorization phase, store generation number in first matrix element
    crowiterator endi=A.end();
    createiterator ci=ILU.createbegin();
    for (crowiterator i=A.begin(); i!=endi; ++i)
    {
      //		std::cout << "in row " << i.index() << std::endl;
      map rowpattern;           // maps column index to generation

      // initialize pattern with row of A
      for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
        rowpattern[j.index()] = 0;

      // eliminate entries in row which are to the left of the diagonal
      for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
      {
        if ((*ik).second<n)
        {
          //                            std::cout << "  eliminating " << i.index() << "," << (*ik).first
          //                                              << " level " << (*ik).second << std::endl;

          coliterator endk = ILU[(*ik).first].end();                       // end of row k
          coliterator kj = ILU[(*ik).first].find((*ik).first);                       // diagonal in k
          for (++kj; kj!=endk; ++kj)                       // row k eliminates in row i
          {
            // we misuse the storage to store an int. If the field_type is std::complex, we have to access the real part
            // starting from C++11, we can use std::real to always return a real value, even if it is double/float
            int generation = (int) std::real( firstmatrixelement(*kj) );
            if (generation<n)
            {
              mapiterator ij = rowpattern.find(kj.index());
              if (ij==rowpattern.end())
              {
                //std::cout << "    new entry " << i.index() << "," << kj.index()
                //                                                << " level " << (*ik).second+1 << std::endl;

                rowpattern[kj.index()] = generation+1;
              }
            }
          }
        }
      }

      // create row
      for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
        ci.insert((*ik).first);
      ++ci;           // now row i exist

      // write generation index into entries
      coliterator endILUij = ILU[i.index()].end();;
      for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
        firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
    }

    //	printmatrix(std::cout,ILU,"ilu pattern","row",10,2);

    // copy entries of A
    for (crowiterator i=A.begin(); i!=endi; ++i)
    {
      coliterator ILUij;
      coliterator endILUij = ILU[i.index()].end();;
      for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
        (*ILUij) = 0;           // clear row
      ccoliterator Aij = (*i).begin();
      ccoliterator endAij = (*i).end();
      ILUij = ILU[i.index()].begin();
      while (Aij!=endAij && ILUij!=endILUij)
      {
        if (Aij.index()==ILUij.index())
        {
          *ILUij = *Aij;
          ++Aij; ++ILUij;
        }
        else
        {
          if (Aij.index()<ILUij.index())
            ++Aij;
          else
            ++ILUij;
        }
      }
    }

    // call decomposition on pattern
    bilu0_decomposition(ILU);
  }


  /** @} end documentation */

} // end namespace

#endif