Dune Core Modules (2.8.0)

ildl.hh
Go to the documentation of this file.
1#ifndef DUNE_ISTL_ILDL_HH
2#define DUNE_ISTL_ILDL_HH
3
6#include "ilu.hh"
7
15namespace Dune
16{
17
18 // bildl_subtractBCT
19 // -----------------
20
21 template< class K, int m, int n >
22 inline static void bildl_subtractBCT ( const FieldMatrix< K, m, n > &B, const FieldMatrix< K, m, n > &CT, FieldMatrix< K, m, n > &A )
23 {
24 for( int i = 0; i < m; ++i )
25 {
26 for( int j = 0; j < n; ++j )
27 {
28 for( int k = 0; k < n; ++k )
29 A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
30 }
31 }
32 }
33
34 template< class K >
35 inline static void bildl_subtractBCT ( const K &B, const K &CT, K &A,
36 typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr )
37 {
38 A -= B * CT;
39 }
40
41 template< class Matrix >
42 inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matrix &A,
43 typename std::enable_if_t<!Dune::IsNumber<Matrix>::value>* sfinae = nullptr )
44 {
45 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
46 {
47 auto &&A_i = *i;
48 auto &&B_i = B[ i.index() ];
49 const auto ikend = B_i.end();
50 for( auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
51 {
52 auto &&A_ij = *j;
53 auto &&CT_j = CT[ j.index() ];
54 const auto jkend = CT_j.end();
55 for( auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
56 {
57 if( ik.index() == jk.index() )
58 {
59 bildl_subtractBCT( *ik, *jk, A_ij );
60 ++ik; ++jk;
61 }
62 else if( ik.index() < jk.index() )
63 ++ik;
64 else
65 ++jk;
66 }
67 }
68 }
69 }
70
71
72
73 // bildl_decompose
74 // ---------------
75
85 template< class Matrix >
86 inline void bildl_decompose ( Matrix &A )
87 {
88 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
89 {
90 auto &&A_i = *i;
91
92 auto ij = A_i.begin();
93 for( ; ij.index() < i.index(); ++ij )
94 {
95 auto &&A_ij = *ij;
96 auto &&A_j = A[ ij.index() ];
97
98 // store L_ij Dj in A_ij (note: for k < i: A_kj = L_kj)
99 // L_ij Dj = A_ij - \sum_{k < j} (L_ik D_k) L_jk^T
100 auto ik = A_i.begin();
101 auto jk = A_j.begin();
102 while( (ik != ij) && (jk.index() < ij.index()) )
103 {
104 if( ik.index() == jk.index() )
105 {
106 bildl_subtractBCT(*ik, *jk, A_ij);
107 ++ik; ++jk;
108 }
109 else if( ik.index() < jk.index() )
110 ++ik;
111 else
112 ++jk;
113 }
114 }
115
116 if( ij.index() != i.index() )
117 DUNE_THROW( ISTLError, "diagonal entry missing" );
118
119 // update diagonal and multiply A_ij by D_j^{-1}
120 auto &&A_ii = *ij;
121 for( auto ik = A_i.begin(); ik != ij; ++ik )
122 {
123 auto &&A_ik = *ik;
124 const auto &A_k = A[ ik.index() ];
125
126 auto B = A_ik;
127 Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
128 bildl_subtractBCT( B, A_ik, A_ii );
129 }
130 try
131 {
132 Impl::asMatrix(A_ii).invert();
133 }
134 catch( const Dune::FMatrixError &e )
135 {
136 DUNE_THROW( MatrixBlockError, "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what(); th__ex.r = i.index(); th__ex.c = ij.index() );
137 }
138 }
139 }
140
141
142
143 // bildl_backsolve
144 // ---------------
145
146 template< class Matrix, class X, class Y >
147 inline void bildl_backsolve ( const Matrix &A, X &v, const Y &d, bool isLowerTriangular = false )
148 {
149 // solve L v = d, note: Lii = I
150 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
151 {
152 const auto &A_i = *i;
153 v[ i.index() ] = d[ i.index() ];
154 for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
155 {
156 auto&& vi = Impl::asVector( v[ i.index() ] );
157 Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
158 }
159 }
160
161 // solve D w = v, note: diagonal stores Dii^{-1}
162 if( isLowerTriangular )
163 {
164 // The matrix is lower triangular, so the diagonal entry is the
165 // last one in each row.
166 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
167 {
168 const auto &A_i = *i;
169 const auto ii = A_i.beforeEnd();
170 assert( ii.index() == i.index() );
171 // We need to be careful here: Directly using
172 // auto rhs = Impl::asVector(v[ i.index() ]);
173 // is not OK in case this is a proxy. Hence
174 // we first have to copy the value. Notice that
175 // this is still not OK, if the vector type itself returns
176 // proxy references.
177 auto rhsValue = v[ i.index() ];
178 auto&& rhs = Impl::asVector(rhsValue);
179 auto&& vi = Impl::asVector( v[ i.index() ] );
180 Impl::asMatrix(*ii).mv(rhs, vi);
181 }
182 }
183 else
184 {
185 // Without assumptions on the sparsity pattern we have to search
186 // for the diagonal entry in each row.
187 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
188 {
189 const auto &A_i = *i;
190 const auto ii = A_i.find( i.index() );
191 assert( ii.index() == i.index() );
192 // We need to be careful here: Directly using
193 // auto rhs = Impl::asVector(v[ i.index() ]);
194 // is not OK in case this is a proxy. Hence
195 // we first have to copy the value. Notice that
196 // this is still not OK, if the vector type itself returns
197 // proxy references.
198 auto rhsValue = v[ i.index() ];
199 auto&& rhs = Impl::asVector(rhsValue);
200 auto&& vi = Impl::asVector( v[ i.index() ] );
201 Impl::asMatrix(*ii).mv(rhs, vi);
202 }
203 }
204
205 // solve L^T v = w, note: only L is stored
206 // note: we perform the operation column-wise from right to left
207 for( auto i = A.beforeEnd(), iend = A.beforeBegin(); i != iend; --i )
208 {
209 const auto &A_i = *i;
210 for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
211 {
212 auto&& vij = Impl::asVector( v[ ij.index() ] );
213 Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
214 }
215 }
216 }
217
218} // namespace Dune
219
220#endif // #ifndef DUNE_ISTL_ILDL_HH
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:151
derive error class from the base class in common
Definition: istlexception.hh:17
Error when performing an operation on a matrix block.
Definition: istlexception.hh:59
A generic dynamic dense matrix.
Definition: matrix.hh:559
RowIterator beforeBegin()
Definition: matrix.hh:632
RowIterator beforeEnd()
Definition: matrix.hh:625
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:618
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:612
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:35
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
The incomplete LU factorization kernels.
Dune namespace.
Definition: alignedallocator.hh:11
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition: ildl.hh:86
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)