Dune Core Modules (2.6.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
8
14namespace Dune {
24 template <class B, class A=std::allocator<B> >
25 class BTDMatrix : public BCRSMatrix<B,A>
26 {
27 public:
28
29 //===== type definitions and constants
30
32 typedef typename B::field_type field_type;
33
35 typedef B block_type;
36
38 typedef A allocator_type;
39
41 //typedef BCRSMatrix<B,A>::row_type row_type;
42
44 typedef typename A::size_type size_type;
45
47 enum {blocklevel = B::blocklevel+1};
48
50 BTDMatrix() : BCRSMatrix<B,A>() {}
51
52 explicit BTDMatrix(size_type size)
53 : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
54 {
55 // Set number of entries for each row
56 // All rows get three entries, except for the first and the last one
57 for (size_t i=0; i<size; i++)
58 this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
59
61
62 // The actual entries for each row
63 for (size_t i=0; i<size; i++) {
64 if (i>0)
65 this->BCRSMatrix<B,A>::addindex(i, i-1);
66 this->BCRSMatrix<B,A>::addindex(i, i );
67 if (i<size-1)
68 this->BCRSMatrix<B,A>::addindex(i, i+1);
69 }
70
72 }
73
75 void setSize(size_type size)
76 {
77 auto nonZeros = (size==0) ? 0 : size + 2*(size-1);
78 this->BCRSMatrix<B,A>::setSize(size, // rows
79 size, // columns
80 nonZeros);
81
82 // Set number of entries for each row
83 // All rows get three entries, except for the first and the last one
84 for (size_t i=0; i<size; i++)
85 this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
86
88
89 // The actual entries for each row
90 for (size_t i=0; i<size; i++) {
91 if (i>0)
92 this->BCRSMatrix<B,A>::addindex(i, i-1);
93 this->BCRSMatrix<B,A>::addindex(i, i );
94 if (i<size-1)
95 this->BCRSMatrix<B,A>::addindex(i, i+1);
96 }
97
99 }
100
103 this->BCRSMatrix<B,A>::operator=(other);
104 return *this;
105 }
106
110 return *this;
111 }
112
118 template <class V>
119 void solve (V& x, const V& rhs) const {
120
121 // special handling for 1x1 matrices. The generic algorithm doesn't work for them
122 if (this->N()==1) {
123 (*this)[0][0].solve(x[0],rhs[0]);
124 return;
125 }
126
127 // Make copies of the rhs and the right matrix band
128 V d = rhs;
129 std::vector<block_type> c(this->N()-1);
130 for (size_t i=0; i<this->N()-1; i++)
131 c[i] = (*this)[i][i+1];
132
133 /* Modify the coefficients. */
134 block_type a_00_inv = (*this)[0][0];
135 a_00_inv.invert();
136
137 //c[0] /= (*this)[0][0]; /* Division by zero risk. */
138 block_type c_0_tmp = c[0];
139 FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]);
140
141 // d = a^{-1} d /* Division by zero would imply a singular matrix. */
142 typename V::block_type d_0_tmp = d[0];
143 a_00_inv.mv(d_0_tmp,d[0]);
144
145 for (unsigned int i = 1; i < this->N(); i++) {
146
147 // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
148 block_type tmp;
149 FMatrixHelp::multMatrix((*this)[i][i-1],c[i-1], tmp);
150 block_type id = (*this)[i][i];
151 id -= tmp;
152 id.invert(); /* Division by zero risk. */
153
154 if (i<c.size()) {
155 // c[i] *= id
156 tmp = c[i];
157 FMatrixHelp::multMatrix(id,tmp, c[i]); /* Last value calculated is redundant. */
158 }
159
160 // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
161 (*this)[i][i-1].mmv(d[i-1], d[i]);
162 typename V::block_type tmpVec = d[i];
163 id.mv(tmpVec, d[i]);
164 //d[i] *= id;
165
166 }
167
168 /* Now back substitute. */
169 x[this->N() - 1] = d[this->N() - 1];
170 for (int i = this->N() - 2; i >= 0; i--) {
171 //x[i] = d[i] - c[i] * x[i + 1];
172 x[i] = d[i];
173 c[i].mmv(x[i+1], x[i]);
174 }
175
176 }
177
178 private:
179
180 // ////////////////////////////////////////////////////////////////////////////
181 // The following methods from the base class should now actually be called
182 // ////////////////////////////////////////////////////////////////////////////
183
184 // createbegin and createend should be in there, too, but I can't get it to compile
185 // BCRSMatrix<B,A>::CreateIterator createbegin () {}
186 // BCRSMatrix<B,A>::CreateIterator createend () {}
187 void setrowsize (size_type i, size_type s) {}
188 void incrementrowsize (size_type i) {}
189 void endrowsizes () {}
190 void addindex (size_type row, size_type col) {}
191 void endindices () {}
192 };
195} // end namespace Dune
196
197#endif
Implementation of the BCRSMatrix class.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1108
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:489
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1150
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1207
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1894
void setSize(size_type rows, size_type columns, size_type nnz=0)
Set the size of the matrix.
Definition: bcrsmatrix.hh:820
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:870
A block-tridiagonal matrix.
Definition: btdmatrix.hh:26
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:119
A::size_type size_type
implement row_type with compressed vector
Definition: btdmatrix.hh:44
A allocator_type
export the allocator type
Definition: btdmatrix.hh:38
B::field_type field_type
export the type representing the field
Definition: btdmatrix.hh:32
B block_type
export the type representing the components
Definition: btdmatrix.hh:35
BTDMatrix & operator=(const BTDMatrix &other)
assignment
Definition: btdmatrix.hh:102
BTDMatrix()
Default constructor.
Definition: btdmatrix.hh:50
void setSize(size_type size)
Resize the matrix. Invalidates the content!
Definition: btdmatrix.hh:75
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Dune namespace.
Definition: alignedallocator.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)