dune-istl  2.8.0
btdmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_BTDMATRIX_HH
4 #define DUNE_ISTL_BTDMATRIX_HH
5 
6 #include <dune/common/fmatrix.hh>
7 #include <dune/common/scalarvectorview.hh>
8 #include <dune/common/scalarmatrixview.hh>
10 #include <dune/istl/blocklevel.hh>
11 
17 namespace Dune {
27  template <class B, class A=std::allocator<B> >
28  class BTDMatrix : public BCRSMatrix<B,A>
29  {
30  public:
31 
32  //===== type definitions and constants
33 
35  using field_type = typename Imp::BlockTraits<B>::field_type;
36 
38  typedef B block_type;
39 
41  typedef A allocator_type;
42 
44  //typedef BCRSMatrix<B,A>::row_type row_type;
45 
47  typedef typename A::size_type size_type;
48 
50  [[deprecated("Use free blockLevel function. Will be removed after 2.8.")]]
51  static constexpr auto blocklevel = blockLevel<B>()+1;
52 
54  BTDMatrix() : BCRSMatrix<B,A>() {}
55 
56  explicit BTDMatrix(size_type size)
57  : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
58  {
59  // Set number of entries for each row
60  // All rows get three entries, except for the first and the last one
61  for (size_t i=0; i<size; i++)
62  this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
63 
65 
66  // The actual entries for each row
67  for (size_t i=0; i<size; i++) {
68  if (i>0)
69  this->BCRSMatrix<B,A>::addindex(i, i-1);
70  this->BCRSMatrix<B,A>::addindex(i, i );
71  if (i<size-1)
72  this->BCRSMatrix<B,A>::addindex(i, i+1);
73  }
74 
76  }
77 
79  void setSize(size_type size)
80  {
81  auto nonZeros = (size==0) ? 0 : size + 2*(size-1);
82  this->BCRSMatrix<B,A>::setSize(size, // rows
83  size, // columns
84  nonZeros);
85 
86  // Set number of entries for each row
87  // All rows get three entries, except for the first and the last one
88  for (size_t i=0; i<size; i++)
89  this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
90 
92 
93  // The actual entries for each row
94  for (size_t i=0; i<size; i++) {
95  if (i>0)
96  this->BCRSMatrix<B,A>::addindex(i, i-1);
97  this->BCRSMatrix<B,A>::addindex(i, i );
98  if (i<size-1)
99  this->BCRSMatrix<B,A>::addindex(i, i+1);
100  }
101 
103  }
104 
106  BTDMatrix& operator= (const BTDMatrix& other) {
107  this->BCRSMatrix<B,A>::operator=(other);
108  return *this;
109  }
110 
114  return *this;
115  }
116 
122  template <class V>
123  void solve (V& x, const V& rhs) const {
124 
125  // special handling for 1x1 matrices. The generic algorithm doesn't work for them
126  if (this->N()==1) {
127  auto&& x0 = Impl::asVector(x[0]);
128  auto&& rhs0 = Impl::asVector(rhs[0]);
129  Impl::asMatrix((*this)[0][0]).solve(x0, rhs0);
130  return;
131  }
132 
133  // Make copies of the rhs and the right matrix band
134  V d = rhs;
135  std::vector<block_type> c(this->N()-1);
136  for (size_t i=0; i<this->N()-1; i++)
137  c[i] = (*this)[i][i+1];
138 
139  /* Modify the coefficients. */
140  block_type a_00_inv = (*this)[0][0];
141  Impl::asMatrix(a_00_inv).invert();
142 
143  //c[0] /= (*this)[0][0]; /* Division by zero risk. */
144  block_type tmp = a_00_inv;
145  Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[0]));
146  c[0] = tmp;
147 
148  // d = a^{-1} d /* Division by zero would imply a singular matrix. */
149  auto d_0_tmp = d[0];
150  auto&& d_0 = Impl::asVector(d[0]);
151  Impl::asMatrix(a_00_inv).mv(Impl::asVector(d_0_tmp),d_0);
152 
153  for (unsigned int i = 1; i < this->N(); i++) {
154 
155  // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
156  block_type tmp;
157  tmp = (*this)[i][i-1];
158  Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[i-1]));
159 
160  block_type id = (*this)[i][i];
161  id -= tmp;
162  Impl::asMatrix(id).invert(); /* Division by zero risk. */
163 
164  if (i<c.size()) {
165  Impl::asMatrix(c[i]).leftmultiply(Impl::asMatrix(id)); /* Last value calculated is redundant. */
166  }
167 
168  // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
169  auto&& d_i = Impl::asVector(d[i]);
170  Impl::asMatrix((*this)[i][i-1]).mmv(Impl::asVector(d[i-1]), d_i);
171  auto tmpVec = d[i];
172  Impl::asMatrix(id).mv(Impl::asVector(tmpVec), d_i);
173  }
174 
175  /* Now back substitute. */
176  x[this->N() - 1] = d[this->N() - 1];
177  for (int i = this->N() - 2; i >= 0; i--) {
178  //x[i] = d[i] - c[i] * x[i + 1];
179  x[i] = d[i];
180  auto&& x_i = Impl::asVector(x[i]);
181  Impl::asMatrix(c[i]).mmv(Impl::asVector(x[i+1]), x_i);
182  }
183 
184  }
185 
186  private:
187 
188  // ////////////////////////////////////////////////////////////////////////////
189  // The following methods from the base class should now actually be called
190  // ////////////////////////////////////////////////////////////////////////////
191 
192  // createbegin and createend should be in there, too, but I can't get it to compile
193  // BCRSMatrix<B,A>::CreateIterator createbegin () {}
194  // BCRSMatrix<B,A>::CreateIterator createend () {}
195  void setrowsize (size_type i, size_type s) {}
196  void incrementrowsize (size_type i) {}
197  void endrowsizes () {}
198  void addindex (size_type row, size_type col) {}
199  void endindices () {}
200  };
203 } // end namespace Dune
204 
205 #endif
Implementation of the BCRSMatrix class.
Helper functions for determining the vector/matrix block level.
Col col
Definition: matrixmatrix.hh:349
Definition: allocator.hh:9
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:909
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1147
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:528
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1189
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1246
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1970
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:859
A block-tridiagonal matrix.
Definition: btdmatrix.hh:29
static constexpr auto blocklevel
increment block level counter
Definition: btdmatrix.hh:51
BTDMatrix(size_type size)
Definition: btdmatrix.hh:56
void solve(V &x, const V &rhs) const
Use the Thomas algorithm to solve the system Ax=b in O(n) time.
Definition: btdmatrix.hh:123
BTDMatrix & operator=(const BTDMatrix &other)
assignment
Definition: btdmatrix.hh:106
A::size_type size_type
implement row_type with compressed vector
Definition: btdmatrix.hh:47
A allocator_type
export the allocator type
Definition: btdmatrix.hh:41
B block_type
export the type representing the components
Definition: btdmatrix.hh:38
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: btdmatrix.hh:35
BTDMatrix()
Default constructor.
Definition: btdmatrix.hh:54
void setSize(size_type size)
Resize the matrix. Invalidates the content!
Definition: btdmatrix.hh:79