Dune Core Modules (2.6.0)

ildl.hh
Go to the documentation of this file.
1#ifndef DUNE_ISTL_ILDL_HH
2#define DUNE_ISTL_ILDL_HH
3
4#include "ilu.hh"
5
13namespace Dune
14{
15
16 // bildl_subtractBCT
17 // -----------------
18
19 template< class K, int m, int n >
20 inline static void bildl_subtractBCT ( const FieldMatrix< K, m, n > &B, const FieldMatrix< K, m, n > &CT, FieldMatrix< K, m, n > &A )
21 {
22 for( int i = 0; i < m; ++i )
23 {
24 for( int j = 0; j < n; ++j )
25 {
26 for( int k = 0; k < n; ++k )
27 A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
28 }
29 }
30 }
31
32 template< class Matrix >
33 inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matrix &A )
34 {
35 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
36 {
37 auto &&A_i = *i;
38 auto &&B_i = B[ i.index() ];
39 const auto ikend = B_i.end();
40 for( auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
41 {
42 auto &&A_ij = *j;
43 auto &&CT_j = CT[ j.index() ];
44 const auto jkend = CT_j.end();
45 for( auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
46 {
47 if( ik.index() == jk.index() )
48 {
49 bildl_subtractBCT( *ik, *jk, A_ij );
50 ++ik; ++jk;
51 }
52 else if( ik.index() < jk.index() )
53 ++ik;
54 else
55 ++jk;
56 }
57 }
58 }
59 }
60
61
62
63 // bildl_decompose
64 // ---------------
65
75 template< class Matrix >
76 inline void bildl_decompose ( Matrix &A )
77 {
78 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
79 {
80 auto &&A_i = *i;
81
82 auto ij = A_i.begin();
83 for( ; ij.index() < i.index(); ++ij )
84 {
85 auto &&A_ij = *ij;
86 auto &&A_j = A[ ij.index() ];
87
88 // store L_ij Dj in A_ij (note: for k < i: A_kj = L_kj)
89 // L_ij Dj = A_ij - \sum_{k < j} (L_ik D_k) L_jk^T
90 auto ik = A_i.begin();
91 auto jk = A_j.begin();
92 while( (ik != ij) && (jk.index() < ij.index()) )
93 {
94 if( ik.index() == jk.index() )
95 {
96 bildl_subtractBCT(*ik, *jk, A_ij);
97 ++ik; ++jk;
98 }
99 else if( ik.index() < jk.index() )
100 ++ik;
101 else
102 ++jk;
103 }
104 }
105
106 if( ij.index() != i.index() )
107 DUNE_THROW( ISTLError, "diagonal entry missing" );
108
109 // update diagonal and multiply A_ij by D_j^{-1}
110 auto &&A_ii = *ij;
111 for( auto ik = A_i.begin(); ik != ij; ++ik )
112 {
113 auto &&A_ik = *ik;
114 const auto &A_k = A[ ik.index() ];
115
116 auto B = A_ik;
117 A_ik.rightmultiply( *A_k.find( ik.index() ) );
118 bildl_subtractBCT( B, A_ik, A_ii );
119 }
120 try
121 {
122 A_ii.invert();
123 }
124 catch( const Dune::FMatrixError &e )
125 {
126 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() );
127 }
128 }
129 }
130
131
132
133 // bildl_backsolve
134 // ---------------
135
136 template< class Matrix, class X, class Y >
137 inline void bildl_backsolve ( const Matrix &A, X &v, const Y &d, bool isLowerTriangular = false )
138 {
139 // solve L v = d, note: Lii = I
140 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
141 {
142 const auto &A_i = *i;
143 v[ i.index() ] = d[ i.index() ];
144 for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
145 (*ij).mmv( v[ ij.index() ], v[ i.index() ] );
146 }
147
148 // solve D w = v, note: diagonal stores Dii^{-1}
149 if( isLowerTriangular )
150 {
151 // The matrix is lower triangular, so the diagonal entry is the
152 // last one in each row.
153 for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
154 {
155 const auto &A_i = *i;
156 const auto ii = A_i.beforeEnd();
157 assert( ii.index() == i.index() );
158 auto rhs = v[ i.index() ];
159 ii->mv( rhs, v[ i.index() ] );
160 }
161 }
162 else
163 {
164 // Without assumptions on the sparsity pattern we have to search
165 // for the diagonal entry 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.find( i.index() );
170 assert( ii.index() == i.index() );
171 auto rhs = v[ i.index() ];
172 ii->mv( rhs, v[ i.index() ] );
173 }
174 }
175
176 // solve L^T v = w, note: only L is stored
177 // note: we perform the operation column-wise from right to left
178 for( auto i = A.beforeEnd(), iend = A.beforeBegin(); i != iend; --i )
179 {
180 const auto &A_i = *i;
181 for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
182 (*ij).mmtv( v[ i.index() ], v[ ij.index() ] );
183 }
184 }
185
186} // namespace Dune
187
188#endif // #ifndef DUNE_ISTL_ILDL_HH
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:146
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:555
RowIterator beforeBegin()
Definition: matrix.hh:629
RowIterator beforeEnd()
Definition: matrix.hh:622
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:35
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune namespace.
Definition: alignedallocator.hh:10
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition: ildl.hh:76
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)