Dune Core Modules (2.4.1)

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(int size)
53 : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
54 {
55 // special handling for 1x1 matrices
56 if (size==1) {
57
60
61 this->BCRSMatrix<B,A>::addindex(0, 0);
63
64 return;
65 }
66
67 // Set number of entries for each row
69
70 for (int i=1; i<size-1; i++)
72
73 this->BCRSMatrix<B,A>::setrowsize(size-1, 2);
74
76
77 // The actual entries for each row
78 this->BCRSMatrix<B,A>::addindex(0, 0);
79 this->BCRSMatrix<B,A>::addindex(0, 1);
80
81 for (int i=1; i<size-1; i++) {
82 this->BCRSMatrix<B,A>::addindex(i, i-1);
83 this->BCRSMatrix<B,A>::addindex(i, i );
84 this->BCRSMatrix<B,A>::addindex(i, i+1);
85 }
86
87 this->BCRSMatrix<B,A>::addindex(size-1, size-2);
88 this->BCRSMatrix<B,A>::addindex(size-1, size-1);
89
91
92 }
93
96 this->BCRSMatrix<B,A>::operator=(other);
97 return *this;
98 }
99
103 return *this;
104 }
105
111 template <class V>
112 void solve (V& x, const V& rhs) const {
113
114 // special handling for 1x1 matrices. The generic algorithm doesn't work for them
115 if (this->N()==1) {
116 (*this)[0][0].solve(x[0],rhs[0]);
117 return;
118 }
119
120 // Make copies of the rhs and the right matrix band
121 V d = rhs;
122 std::vector<block_type> c(this->N()-1);
123 for (size_t i=0; i<this->N()-1; i++)
124 c[i] = (*this)[i][i+1];
125
126 /* Modify the coefficients. */
127 block_type a_00_inv = (*this)[0][0];
128 a_00_inv.invert();
129
130 //c[0] /= (*this)[0][0]; /* Division by zero risk. */
131 block_type c_0_tmp = c[0];
132 FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]);
133
134 // d = a^{-1} d /* Division by zero would imply a singular matrix. */
135 typename V::block_type d_0_tmp = d[0];
136 a_00_inv.mv(d_0_tmp,d[0]);
137
138 for (unsigned int i = 1; i < this->N(); i++) {
139
140 // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
141 block_type tmp;
142 FMatrixHelp::multMatrix((*this)[i][i-1],c[i-1], tmp);
143 block_type id = (*this)[i][i];
144 id -= tmp;
145 id.invert(); /* Division by zero risk. */
146
147 if (i<c.size()) {
148 // c[i] *= id
149 tmp = c[i];
150 FMatrixHelp::multMatrix(id,tmp, c[i]); /* Last value calculated is redundant. */
151 }
152
153 // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
154 (*this)[i][i-1].mmv(d[i-1], d[i]);
155 typename V::block_type tmpVec = d[i];
156 id.mv(tmpVec, d[i]);
157 //d[i] *= id;
158
159 }
160
161 /* Now back substitute. */
162 x[this->N() - 1] = d[this->N() - 1];
163 for (int i = this->N() - 2; i >= 0; i--) {
164 //x[i] = d[i] - c[i] * x[i + 1];
165 x[i] = d[i];
166 c[i].mmv(x[i+1], x[i]);
167 }
168
169 }
170
171 private:
172
173 // ////////////////////////////////////////////////////////////////////////////
174 // The following methods from the base class should now actually be called
175 // ////////////////////////////////////////////////////////////////////////////
176
177 // createbegin and createend should be in there, too, but I can't get it to compile
178 // BCRSMatrix<B,A>::CreateIterator createbegin () {}
179 // BCRSMatrix<B,A>::CreateIterator createend () {}
180 void setrowsize (size_type i, size_type s) {}
181 void incrementrowsize (size_type i) {}
182 void endrowsizes () {}
183 void addindex (size_type row, size_type col) {}
184 void endindices () {}
185 };
188} // end namespace Dune
189
190#endif
Implementation of the BCRSMatrix class.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
void endrowsizes()
indicate that size of all rows is defined
Definition: bcrsmatrix.hh:1109
void setrowsize(size_type i, size_type s)
set number of indices in row i to s
Definition: bcrsmatrix.hh:1077
@ random
Build entries randomly.
Definition: bcrsmatrix.hh:479
void addindex(size_type row, size_type col)
add index (row,col) to the matrix
Definition: bcrsmatrix.hh:1151
void endindices()
indicate that all indices are defined, check consistency
Definition: bcrsmatrix.hh:1208
size_type N() const
number of rows (counted in blocks)
Definition: bcrsmatrix.hh:1849
BCRSMatrix & operator=(const BCRSMatrix &Mat)
assignment
Definition: bcrsmatrix.hh:868
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:112
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:95
BTDMatrix()
Default constructor.
Definition: btdmatrix.hh:50
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Dune namespace.
Definition: alignment.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)