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