5#ifndef DUNE_ISTL_DILU_HH
6#define DUNE_ISTL_DILU_HH
17#include "istlexception.hh"
54 for (
auto row = A.begin(); row != endi; ++row)
56 const auto row_i = row.index();
57 const auto col = row->find(row_i);
62 for (
auto row = A.begin(); row != endi; ++row)
64 const auto row_i = row.index();
65 for (
auto a_ij = row->begin(); a_ij.index() < row_i; ++a_ij)
67 const auto col_j = a_ij.index();
68 const auto a_ji = A[col_j].find(row_i);
70 if (a_ji != A[col_j].end())
73 Dinv_[row_i] -= (*a_ij) * Dinv_[col_j] * (*a_ji);
80 Impl::asMatrix(Dinv_[row_i]).invert();
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)
109 using dblock =
typename Y::block_type;
110 using vblock =
typename X::block_type;
114 for (
auto row = A.begin(); row != endi; ++row)
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)
123 const auto col_j = a_ij.index();
124 Impl::asMatrix(*a_ij).mmv(v[col_j], rhs);
128 auto &&vi = Impl::asVector(v[row_i]);
129 Impl::asMatrix(Dinv_[row_i]).mv(rhs, vi);
133 auto rendi = A.beforeBegin();
134 for (
auto row = A.beforeEnd(); row != rendi; --row)
136 const auto row_i = row.index();
139 for (
auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij)
143 const auto col_j = a_ij.index();
144 Impl::asMatrix(*a_ij).umv(v[col_j], rhs);
149 auto &&vi = Impl::asVector(v[row_i]);
150 Impl::asMatrix(Dinv_[row_i]).mmv(rhs, vi);
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.