DUNE-FEM (unstable)

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