This file is indexed.

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

#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>

/** \file
    \author Oliver Sander
    \brief Implementation of the BTDMatrix class
 */

namespace Dune {
  /**
   * @addtogroup ISTL_SPMV
   * @{
   */
  /** \brief A block-tridiagonal matrix

     \todo It would be safer and more efficient to have a real implementation of
     a block-tridiagonal matrix and not just subclassing from BCRSMatrix.  But that's
     quite a lot of work for that little advantage.*/
  template <class B, class A=std::allocator<B> >
  class BTDMatrix : public BCRSMatrix<B,A>
  {
  public:

    //===== type definitions and constants

    //! export the type representing the field
    typedef typename B::field_type field_type;

    //! export the type representing the components
    typedef B block_type;

    //! export the allocator type
    typedef A allocator_type;

    //! implement row_type with compressed vector
    //typedef BCRSMatrix<B,A>::row_type row_type;

    //! The type for the index access and the size
    typedef typename A::size_type size_type;

    //! increment block level counter
    enum {blocklevel = B::blocklevel+1};

    /** \brief Default constructor */
    BTDMatrix() : BCRSMatrix<B,A>() {}

    explicit BTDMatrix(int size)
      : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
    {
      // special handling for 1x1 matrices
      if (size==1) {

        this->BCRSMatrix<B,A>::setrowsize(0, 1);
        this->BCRSMatrix<B,A>::endrowsizes();

        this->BCRSMatrix<B,A>::addindex(0, 0);
        this->BCRSMatrix<B,A>::endindices();

        return;
      }

      // Set number of entries for each row
      this->BCRSMatrix<B,A>::setrowsize(0, 2);

      for (int i=1; i<size-1; i++)
        this->BCRSMatrix<B,A>::setrowsize(i, 3);

      this->BCRSMatrix<B,A>::setrowsize(size-1, 2);

      this->BCRSMatrix<B,A>::endrowsizes();

      // The actual entries for each row
      this->BCRSMatrix<B,A>::addindex(0, 0);
      this->BCRSMatrix<B,A>::addindex(0, 1);

      for (int i=1; i<size-1; i++) {
        this->BCRSMatrix<B,A>::addindex(i, i-1);
        this->BCRSMatrix<B,A>::addindex(i, i  );
        this->BCRSMatrix<B,A>::addindex(i, i+1);
      }

      this->BCRSMatrix<B,A>::addindex(size-1, size-2);
      this->BCRSMatrix<B,A>::addindex(size-1, size-1);

      this->BCRSMatrix<B,A>::endindices();

    }

    //! assignment
    BTDMatrix& operator= (const BTDMatrix& other) {
      this->BCRSMatrix<B,A>::operator=(other);
      return *this;
    }

    //! assignment from scalar
    BTDMatrix& operator= (const field_type& k) {
      this->BCRSMatrix<B,A>::operator=(k);
      return *this;
    }

    /** \brief Use the Thomas algorithm to solve the system Ax=b in O(n) time
     *
     * \exception ISTLError if the matrix is singular
     *
     */
    template <class V>
    void solve (V& x, const V& rhs) const {

      // special handling for 1x1 matrices.  The generic algorithm doesn't work for them
      if (this->N()==1) {
        (*this)[0][0].solve(x[0],rhs[0]);
        return;
      }

      // Make copies of the rhs and the right matrix band
      V d = rhs;
      std::vector<block_type> c(this->N()-1);
      for (size_t i=0; i<this->N()-1; i++)
        c[i] = (*this)[i][i+1];

      /* Modify the coefficients. */
      block_type a_00_inv = (*this)[0][0];
      a_00_inv.invert();

      //c[0] /= (*this)[0][0];	/* Division by zero risk. */
      block_type c_0_tmp = c[0];
      FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]);

      // d = a^{-1} d        /* Division by zero would imply a singular matrix. */
      typename V::block_type d_0_tmp = d[0];
      a_00_inv.mv(d_0_tmp,d[0]);

      for (unsigned int i = 1; i < this->N(); i++) {

        // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
        block_type tmp;
        FMatrixHelp::multMatrix((*this)[i][i-1],c[i-1], tmp);
        block_type id = (*this)[i][i];
        id -= tmp;
        id.invert();         /* Division by zero risk. */

        if (i<c.size()) {
          // c[i] *= id
          tmp = c[i];
          FMatrixHelp::multMatrix(id,tmp, c[i]);                      /* Last value calculated is redundant. */
        }

        // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
        (*this)[i][i-1].mmv(d[i-1], d[i]);
        typename V::block_type tmpVec = d[i];
        id.mv(tmpVec, d[i]);
        //d[i] *= id;

      }

      /* Now back substitute. */
      x[this->N() - 1] = d[this->N() - 1];
      for (int i = this->N() - 2; i >= 0; i--) {
        //x[i] = d[i] - c[i] * x[i + 1];
        x[i] = d[i];
        c[i].mmv(x[i+1], x[i]);
      }

    }

  private:

    // ////////////////////////////////////////////////////////////////////////////
    //   The following methods from the base class should now actually be called
    // ////////////////////////////////////////////////////////////////////////////

    // createbegin and createend should be in there, too, but I can't get it to compile
    //     BCRSMatrix<B,A>::CreateIterator createbegin () {}
    //     BCRSMatrix<B,A>::CreateIterator createend () {}
    void setrowsize (size_type i, size_type s) {}
    void incrementrowsize (size_type i) {}
    void endrowsizes () {}
    void addindex (size_type row, size_type col) {}
    void endindices () {}
  };
  /** @}*/

}  // end namespace Dune

#endif