DUNE PDELab (git)

dilu.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// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_DILU_HH
6#define DUNE_ISTL_DILU_HH
7
8#include <cmath>
9#include <complex>
10#include <map>
11#include <vector>
12
16
17#include "istlexception.hh"
18
23namespace Dune
24{
25
30 namespace DILU
31 {
32
50 template <class M>
51 void blockDILUDecomposition(M &A, std::vector<typename M::block_type> &Dinv_)
52 {
53 auto endi = A.end();
54 for (auto row = A.begin(); row != endi; ++row)
55 {
56 const auto row_i = row.index();
57 const auto col = row->find(row_i);
58 // initialise Dinv[i] = A[i, i]
59 Dinv_[row_i] = *col;
60 }
61
62 for (auto row = A.begin(); row != endi; ++row)
63 {
64 const auto row_i = row.index();
65 for (auto a_ij = row->begin(); a_ij.index() < row_i; ++a_ij)
66 {
67 const auto col_j = a_ij.index();
68 const auto a_ji = A[col_j].find(row_i);
69 // if A[i, j] != 0 and A[j, i] != 0
70 if (a_ji != A[col_j].end())
71 {
72 // Dinv[i] -= A[i, j] * d[j] * A[j, i]
73 Dinv_[row_i] -= (*a_ij) * Dinv_[col_j] * (*a_ji);
74 }
75 }
76
77 // store the inverse
78 try
79 {
80 Impl::asMatrix(Dinv_[row_i]).invert(); // compute inverse of diagonal block
81 }
82 catch (Dune::FMatrixError &e)
83 {
84 DUNE_THROW(MatrixBlockError, "DILU failed to invert matrix block D[" << row_i << "]"
85 << e.what();
86 th__ex.r = row_i;);
87 }
88 }
89 }
90
106 template <class M, class X, class Y>
107 void blockDILUBacksolve(const M &A, const std::vector<typename M::block_type> Dinv_, X &v, const Y &d)
108 {
109 using dblock = typename Y::block_type;
110 using vblock = typename X::block_type;
111
112 // lower triangular solve: (D + L_A) y = d
113 auto endi = A.end();
114 for (auto row = A.begin(); row != endi; ++row)
115 {
116 const auto row_i = row.index();
117 dblock rhsValue(d[row_i]);
118 auto &&rhs = Impl::asVector(rhsValue);
119 for (auto a_ij = (*row).begin(); a_ij.index() < row_i; ++a_ij)
120 {
121 // if A[i][j] != 0
122 // rhs -= A[i][j]* y[j], where v_j stores y_j
123 const auto col_j = a_ij.index();
124 Impl::asMatrix(*a_ij).mmv(v[col_j], rhs);
125 }
126 // y_i = Dinv_i * rhs
127 // storing y_i in v_i
128 auto &&vi = Impl::asVector(v[row_i]);
129 Impl::asMatrix(Dinv_[row_i]).mv(rhs, vi); // (D + L_A)_ii = D_i
130 }
131
132 // upper triangular solve: (D + U_A) v = Dy
133 auto rendi = A.beforeBegin();
134 for (auto row = A.beforeEnd(); row != rendi; --row)
135 {
136 const auto row_i = row.index();
137 // rhs = 0
138 vblock rhs(0.0);
139 for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij)
140 {
141 // if A[i][j] != 0
142 // rhs += A[i][j]*v[j]
143 const auto col_j = a_ij.index();
144 Impl::asMatrix(*a_ij).umv(v[col_j], rhs);
145 }
146 // calculate update v = M^-1*d
147 // v_i = y_i - Dinv_i*rhs
148 // before update v_i is y_i
149 auto &&vi = Impl::asVector(v[row_i]);
150 Impl::asMatrix(Dinv_[row_i]).mmv(rhs, vi);
151 }
152 }
153 } // end namespace DILU
154
157} // end namespace
158
159#endif
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:126
Error when performing an operation on a matrix block.
Definition: istlexception.hh:52
void blockDILUDecomposition(M &A, std::vector< typename M::block_type > &Dinv_)
Definition: dilu.hh:51
void blockDILUBacksolve(const M &A, const std::vector< typename M::block_type > Dinv_, X &v, const Y &d)
Definition: dilu.hh:107
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
Dune namespace.
Definition: alignedallocator.hh:13
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)