This file is indexed.

/usr/include/dune/istl/ilusubdomainsolver.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
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
#define DUNE_ISTL_ILUSUBDOMAINSOLVER_HH

#include <map>
#include <dune/common/typetraits.hh>
#include "matrix.hh"
#include <cmath>
#include <cstdlib>

namespace Dune {

  /**
   * @file
   * @brief Various local subdomain solvers based on ILU
   * for SeqOverlappingSchwarz.
   * @author Markus Blatt
   */
  /**
   * @addtogroup ISTL
   * @{
   */
  /**
   * @brief base class encapsulating common algorithms of ILU0SubdomainSolver
   * and ILUNSubdomainSolver.
   * @tparam M The type of the matrix.
   * @tparam X The type of the vector for the domain.
   * @tparam X The type of the vector for the range.
   *
   */
  template<class M, class X, class Y>
  class ILUSubdomainSolver  {
  public:
    //! \brief The matrix type the preconditioner is for.
    typedef typename Dune::remove_const<M>::type 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 Apply the subdomain solver.
     *
     * On entry v=? and d=b-A(x) (although this might not be
     * computed in that way. On exit v contains the update
     */
    virtual void apply (X& v, const Y& d) =0;

    virtual ~ILUSubdomainSolver()
    {}

  protected:
    /**
     * @brief Copy the local part of the global matrix to ILU.
     * @param A The global matrix.
     * @param rowset The global indices of the local problem.
     */
    template<class S>
    std::size_t copyToLocalMatrix(const M& A, S& rowset);

    //! \brief The ILU0 decomposition of the matrix, or the local matrix
    // for ILUN
    matrix_type ILU;
  };

  /**
   * @brief Exact subdomain solver using ILU(p) with appropriate p.
   * @tparam M The type of the matrix.
   * @tparam X The type of the vector for the domain.
   * @tparam X The type of the vector for the range.
   */
  template<class M, class X, class Y>
  class ILU0SubdomainSolver
    : public ILUSubdomainSolver<M,X,Y>{
  public:
    //! \brief The matrix type the preconditioner is for.
    typedef typename Dune::remove_const<M>::type matrix_type;
    typedef typename Dune::remove_const<M>::type rilu_type;
    //! \brief The domain type of the preconditioner.
    typedef X domain_type;
    //! \brief The range type of the preconditioner.
    typedef Y range_type;


    /**
     * @brief Apply the subdomain solver.
     * @copydoc ILUSubdomainSolver::apply
     */
    void apply (X& v, const Y& d)
    {
      bilu_backsolve(this->ILU,v,d);
    }
    /**
     * @brief Set the data of the local problem.
     *
     * @param A The global matrix.
     * @param rowset The global indices of the local problem.
     * @tparam S The type of the set with the indices.
     */
    template<class S>
    void setSubMatrix(const M& A, S& rowset);

  };

  template<class M, class X, class Y>
  class ILUNSubdomainSolver
    : public ILUSubdomainSolver<M,X,Y>{
  public:
    //! \brief The matrix type the preconditioner is for.
    typedef typename Dune::remove_const<M>::type matrix_type;
    typedef typename Dune::remove_const<M>::type rilu_type;
    //! \brief The domain type of the preconditioner.
    typedef X domain_type;
    //! \brief The range type of the preconditioner.
    typedef Y range_type;

    /**
     * @brief Apply the subdomain solver.
     * @copydoc ILUSubdomainSolver::apply
     */
    void apply (X& v, const Y& d)
    {
      bilu_backsolve(RILU,v,d);
    }

    /**
     * @brief Set the data of the local problem.
     *
     * @param A The global matrix.
     * @param rowset The global indices of the local problem.
     * @tparam S The type of the set with the indices.
     */
    template<class S>
    void setSubMatrix(const M& A, S& rowset);

  private:
    /**
     * @brief Storage for the ILUN decomposition.
     */
    rilu_type RILU;
  };



  template<class M, class X, class Y>
  template<class S>
  std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
  {
    // Calculate consecutive indices for local problem
    // while perserving the ordering
    typedef typename M::size_type size_type;
    typedef std::map<typename S::value_type,size_type> IndexMap;
    typedef typename IndexMap::iterator IMIter;
    IndexMap indexMap;
    IMIter guess = indexMap.begin();
    size_type localIndex=0;

    typedef typename S::const_iterator SIter;
    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
        rowIdx!= rowEnd; ++rowIdx, ++localIndex)
      guess = indexMap.insert(guess,
                              std::make_pair(*rowIdx,localIndex));


    // Build Matrix for local subproblem
    ILU.setSize(rowSet.size(),rowSet.size());
    ILU.setBuildMode(matrix_type::row_wise);

    // Create sparsity pattern
    typedef typename matrix_type::CreateIterator CIter;
    CIter rowCreator = ILU.createbegin();
    std::size_t offset=0;
    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
        rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
      // See wich row entries are in our subset and add them to
      // the sparsity pattern
      guess = indexMap.begin();

      for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
          endcol=A[*rowIdx].end(); col != endcol; ++col) {
        // search for the entry in the row set
        guess = indexMap.find(col.index());
        if(guess!=indexMap.end()) {
          // add local index to row
          rowCreator.insert(guess->second);
          offset=std::max(offset,(std::size_t)std::abs((int)(guess->second-rowCreator.index())));
        }
      }

    }

    // Insert the matrix values for the local problem
    typename matrix_type::iterator iluRow=ILU.begin();

    for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
        rowIdx!= rowEnd; ++rowIdx, ++iluRow) {
      // See wich row entries are in our subset and add them to
      // the sparsity pattern
      typename matrix_type::ColIterator localCol=iluRow->begin();
      for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
          endcol=A[*rowIdx].end(); col != endcol; ++col) {
        // search for the entry in the row set
        guess = indexMap.find(col.index());
        if(guess!=indexMap.end()) {
          // set local value
          (*localCol)=(*col);
          ++localCol;
        }
      }
    }
    return offset;
  }


  template<class M, class X, class Y>
  template<class S>
  void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
  {
    this->copyToLocalMatrix(A,rowSet);
    bilu0_decomposition(this->ILU);
  }

  template<class M, class X, class Y>
  template<class S>
  void ILUNSubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
  {
    std::size_t offset=copyToLocalMatrix(A,rowSet);
    RILU.setSize(rowSet.size(),rowSet.size(), (1+2*offset)*rowSet.size());
    RILU.setBuildMode(matrix_type::row_wise);
    bilu_decomposition(this->ILU, (offset+1)/2, RILU);
  }

  /** @} */
} // end name space DUNE


#endif