Dune Core Modules (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
11
17namespace 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
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.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
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
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:909
A block-tridiagonal matrix.
Definition: btdmatrix.hh:29
static constexpr auto blocklevel
increment block level counter
Definition: btdmatrix.hh:51
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
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 & operator=(const BTDMatrix &other)
assignment
Definition: btdmatrix.hh:106
BTDMatrix()
Default constructor.
Definition: btdmatrix.hh:54
void setSize(size_type size)
Resize the matrix. Invalidates the content!
Definition: btdmatrix.hh:79
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Dune namespace.
Definition: alignedallocator.hh:11
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)