DUNE PDELab (git)

btdmatrix.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_BTDMATRIX_HH
6#define DUNE_ISTL_BTDMATRIX_HH
7
13
19namespace Dune {
29 template <class B, class A=std::allocator<B> >
30 class BTDMatrix : public BCRSMatrix<B,A>
31 {
32 public:
33
34 //===== type definitions and constants
35
37 using field_type = typename Imp::BlockTraits<B>::field_type;
38
40 typedef B block_type;
41
43 typedef A allocator_type;
44
46 //typedef BCRSMatrix<B,A>::row_type row_type;
47
49 typedef typename A::size_type size_type;
50
52 BTDMatrix() : BCRSMatrix<B,A>() {}
53
54 explicit BTDMatrix(size_type size)
55 : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
56 {
57 // Set number of entries for each row
58 // All rows get three entries, except for the first and the last one
59 for (size_t i=0; i<size; i++)
60 this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
61
63
64 // The actual entries for each row
65 for (size_t i=0; i<size; i++) {
66 if (i>0)
67 this->BCRSMatrix<B,A>::addindex(i, i-1);
68 this->BCRSMatrix<B,A>::addindex(i, i );
69 if (i<size-1)
70 this->BCRSMatrix<B,A>::addindex(i, i+1);
71 }
72
74 }
75
77 void setSize(size_type size)
78 {
79 auto nonZeros = (size==0) ? 0 : size + 2*(size-1);
80 this->BCRSMatrix<B,A>::setSize(size, // rows
81 size, // columns
82 nonZeros);
83
84 // Set number of entries for each row
85 // All rows get three entries, except for the first and the last one
86 for (size_t i=0; i<size; i++)
87 this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
88
90
91 // The actual entries for each row
92 for (size_t i=0; i<size; i++) {
93 if (i>0)
94 this->BCRSMatrix<B,A>::addindex(i, i-1);
95 this->BCRSMatrix<B,A>::addindex(i, i );
96 if (i<size-1)
97 this->BCRSMatrix<B,A>::addindex(i, i+1);
98 }
99
101 }
102
105 this->BCRSMatrix<B,A>::operator=(other);
106 return *this;
107 }
108
112 return *this;
113 }
114
120 template <class V>
121 void solve (V& x, const V& rhs) const {
122
123 // special handling for 1x1 matrices. The generic algorithm doesn't work for them
124 if (this->N()==1) {
125 auto&& x0 = Impl::asVector(x[0]);
126 auto&& rhs0 = Impl::asVector(rhs[0]);
127 Impl::asMatrix((*this)[0][0]).solve(x0, rhs0);
128 return;
129 }
130
131 // Make copies of the rhs and the right matrix band
132 V d = rhs;
133 std::vector<block_type> c(this->N()-1);
134 for (size_t i=0; i<this->N()-1; i++)
135 c[i] = (*this)[i][i+1];
136
137 /* Modify the coefficients. */
138 block_type a_00_inv = (*this)[0][0];
139 Impl::asMatrix(a_00_inv).invert();
140
141 //c[0] /= (*this)[0][0]; /* Division by zero risk. */
142 block_type tmp = a_00_inv;
143 Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[0]));
144 c[0] = tmp;
145
146 // d = a^{-1} d /* Division by zero would imply a singular matrix. */
147 auto d_0_tmp = d[0];
148 auto&& d_0 = Impl::asVector(d[0]);
149 Impl::asMatrix(a_00_inv).mv(Impl::asVector(d_0_tmp),d_0);
150
151 for (unsigned int i = 1; i < this->N(); i++) {
152
153 // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
154 block_type tmp;
155 tmp = (*this)[i][i-1];
156 Impl::asMatrix(tmp).rightmultiply(Impl::asMatrix(c[i-1]));
157
158 block_type id = (*this)[i][i];
159 id -= tmp;
160 Impl::asMatrix(id).invert(); /* Division by zero risk. */
161
162 if (i<c.size()) {
163 Impl::asMatrix(c[i]).leftmultiply(Impl::asMatrix(id)); /* Last value calculated is redundant. */
164 }
165
166 // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
167 auto&& d_i = Impl::asVector(d[i]);
168 Impl::asMatrix((*this)[i][i-1]).mmv(Impl::asVector(d[i-1]), d_i);
169 auto tmpVec = d[i];
170 Impl::asMatrix(id).mv(Impl::asVector(tmpVec), d_i);
171 }
172
173 /* Now back substitute. */
174 x[this->N() - 1] = d[this->N() - 1];
175 for (int i = this->N() - 2; i >= 0; i--) {
176 //x[i] = d[i] - c[i] * x[i + 1];
177 x[i] = d[i];
178 auto&& x_i = Impl::asVector(x[i]);
179 Impl::asMatrix(c[i]).mmv(Impl::asVector(x[i+1]), x_i);
180 }
181
182 }
183
184 private:
185
186 // ////////////////////////////////////////////////////////////////////////////
187 // The following methods from the base class should now actually be called
188 // ////////////////////////////////////////////////////////////////////////////
189
190 // createbegin and createend should be in there, too, but I can't get it to compile
191 // BCRSMatrix<B,A>::CreateIterator createbegin () {}
192 // BCRSMatrix<B,A>::CreateIterator createend () {}
193 void setrowsize (size_type i, size_type s) {}
194 void incrementrowsize (size_type i) {}
195 void endrowsizes () {}
196 void addindex (size_type row, size_type col) {}
197 void endindices () {}
198 };
199
200 template<typename B, typename A>
201 struct FieldTraits< BTDMatrix<B, A> >
202 {
203 using field_type = typename BTDMatrix<B, A>::field_type;
204 using real_type = typename FieldTraits<field_type>::real_type;
205 };
206
209} // end namespace Dune
210
211#endif
Helper functions for determining the vector/matrix block level.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1146
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:526
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1188
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1269
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:2001
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:858
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:908
A block-tridiagonal matrix.
Definition: btdmatrix.hh:31
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:121
A::size_type size_type
implement row_type with compressed vector
Definition: btdmatrix.hh:49
A allocator_type
export the allocator type
Definition: btdmatrix.hh:43
B block_type
export the type representing the components
Definition: btdmatrix.hh:40
typename Imp::BlockTraits< B >::field_type field_type
export the type representing the field
Definition: btdmatrix.hh:37
BTDMatrix & operator=(const BTDMatrix &other)
assignment
Definition: btdmatrix.hh:104
BTDMatrix()
Default constructor.
Definition: btdmatrix.hh:52
void setSize(size_type size)
Resize the matrix. Invalidates the content!
Definition: btdmatrix.hh:77
Implementation of the BCRSMatrix class.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
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 (Nov 24, 23:30, 2024)