Skip to content
Snippets Groups Projects
Commit d20a4604 authored by Oliver Sander's avatar Oliver Sander
Browse files

Merge branch 'feature/ilu-naming-conventions' into 'master' 

Feature/ilu naming conventions

See merge request core/dune-istl!355
parents bf0582fa 230357e9
No related branches found
No related tags found
1 merge request!355Feature/ilu naming conventions
Pipeline #32576 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
dune/istl/ilu.hh
+ 215
− 202
View file @ d20a4604
......@@ -9,12 +9,14 @@
#include <vector>
#include <dune/common/fmatrix.hh>
#include <dune/common/deprecated.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include "istlexception.hh"
/** \file
* \brief ???
* \brief The incomplete LU factorization kernels
*/
namespace Dune {
......@@ -23,240 +25,232 @@ namespace Dune {
@{
*/
class MatrixBlockError : public virtual Dune::FMatrixError {
public:
int r, c;
};
namespace ILU {
//! compute ILU decomposition of A. A is overwritten by its decomposition
template<class M>
void bilu0_decomposition (M& A)
{
// iterator types
typedef typename M::RowIterator rowiterator;
typedef typename M::ColIterator coliterator;
typedef typename M::block_type block;
// implement left looking variant with stored inverse
rowiterator endi=A.end();
for (rowiterator i=A.begin(); i!=endi; ++i)
//! compute ILU decomposition of A. A is overwritten by its decomposition
template<class M>
void blockILU0Decomposition (M& A)
{
// coliterator is diagonal after the following loop
coliterator endij=(*i).end(); // end of row i
coliterator ij;
// iterator types
typedef typename M::RowIterator rowiterator;
typedef typename M::ColIterator coliterator;
typedef typename M::block_type block;
// eliminate entries left of diagonal; store L factor
for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
// implement left looking variant with stored inverse
rowiterator endi=A.end();
for (rowiterator i=A.begin(); i!=endi; ++i)
{
// find A_jj which eliminates A_ij
coliterator jj = A[ij.index()].find(ij.index());
// compute L_ij = A_jj^-1 * A_ij
Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
// modify row
coliterator endjk=A[ij.index()].end(); // end of row j
coliterator jk=jj; ++jk;
coliterator ik=ij; ++ik;
while (ik!=endij && jk!=endjk)
if (ik.index()==jk.index())
{
block B(*jk);
Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
*ik -= B;
++ik; ++jk;
}
else
{
if (ik.index()<jk.index())
++ik;
// coliterator is diagonal after the following loop
coliterator endij=(*i).end(); // end of row i
coliterator ij;
// eliminate entries left of diagonal; store L factor
for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
{
// find A_jj which eliminates A_ij
coliterator jj = A[ij.index()].find(ij.index());
// compute L_ij = A_jj^-1 * A_ij
Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
// modify row
coliterator endjk=A[ij.index()].end(); // end of row j
coliterator jk=jj; ++jk;
coliterator ik=ij; ++ik;
while (ik!=endij && jk!=endjk)
if (ik.index()==jk.index())
{
block B(*jk);
Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
*ik -= B;
++ik; ++jk;
}
else
++jk;
}
{
if (ik.index()<jk.index())
++ik;
else
++jk;
}
}
// invert pivot and store it in A
if (ij.index()!=i.index())
DUNE_THROW(ISTLError,"diagonal entry missing");
try {
Impl::asMatrix(*ij).invert(); // compute inverse of diagonal block
}
catch (Dune::FMatrixError & e) {
DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
<< i.index() << "][" << ij.index() << "]" << e.what();
th__ex.r=i.index(); th__ex.c=ij.index(););
}
}
}
// invert pivot and store it in A
if (ij.index()!=i.index())
DUNE_THROW(ISTLError,"diagonal entry missing");
try {
Impl::asMatrix(*ij).invert(); // compute inverse of diagonal block
//! LU backsolve with stored inverse
template<class M, class X, class Y>
void blockILUBacksolve (const M& A, X& v, const Y& d)
{
// iterator types
typedef typename M::ConstRowIterator rowiterator;
typedef typename M::ConstColIterator coliterator;
typedef typename Y::block_type dblock;
typedef typename X::block_type vblock;
// lower triangular solve
rowiterator endi=A.end();
for (rowiterator i=A.begin(); i!=endi; ++i)
{
// We need to be careful here: Directly using
// auto rhs = Impl::asVector(d[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
dblock rhsValue(d[i.index()]);
auto&& rhs = Impl::asVector(rhsValue);
for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
Impl::asVector(v[i.index()]) = rhs; // Lii = I
}
catch (Dune::FMatrixError & e) {
DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
<< i.index() << "][" << ij.index() << "]" << e.what();
th__ex.r=i.index(); th__ex.c=ij.index(););
// upper triangular solve
rowiterator rendi=A.beforeBegin();
for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
{
// We need to be careful here: Directly using
// auto rhs = Impl::asVector(v[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
vblock rhsValue(v[i.index()]);
auto&& rhs = Impl::asVector(rhsValue);
coliterator j;
for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
auto&& vi = Impl::asVector(v[i.index()]);
Impl::asMatrix(*j).mv(rhs,vi); // diagonal stores inverse!
}
}
}
//! LU backsolve with stored inverse
template<class M, class X, class Y>
void bilu_backsolve (const M& A, X& v, const Y& d)
{
// iterator types
typedef typename M::ConstRowIterator rowiterator;
typedef typename M::ConstColIterator coliterator;
typedef typename Y::block_type dblock;
typedef typename X::block_type vblock;
// lower triangular solve
rowiterator endi=A.end();
for (rowiterator i=A.begin(); i!=endi; ++i)
// recursive function template to access first entry of a matrix
template<class M>
typename M::field_type& firstMatrixElement (M& A, typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
{
// We need to be careful here: Directly using
// auto rhs = Impl::asVector(d[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
dblock rhsValue(d[i.index()]);
auto&& rhs = Impl::asVector(rhsValue);
for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
Impl::asVector(v[i.index()]) = rhs; // Lii = I
return firstMatrixElement(*(A.begin()->begin()));
}
// upper triangular solve
rowiterator rendi=A.beforeBegin();
for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
template<class K>
K& firstMatrixElement (K& A, typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
{
// We need to be careful here: Directly using
// auto rhs = Impl::asVector(v[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
vblock rhsValue(v[i.index()]);
auto&& rhs = Impl::asVector(rhsValue);
coliterator j;
for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
auto&& vi = Impl::asVector(v[i.index()]);
Impl::asMatrix(*j).mv(rhs,vi); // diagonal stores inverse!
return A;
}
}
template<class K, int n, int m>
K& firstMatrixElement (FieldMatrix<K,n,m>& A)
{
return A[0][0];
}
// recursive function template to access first entry of a matrix
template<class M>
typename M::field_type& firstmatrixelement (M& A, typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
{
return firstmatrixelement(*(A.begin()->begin()));
}
template<class K>
K& firstmatrixelement (K& A, typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
{
return A;
}
template<class K, int n, int m>
K& firstmatrixelement (FieldMatrix<K,n,m>& A)
{
return A[0][0];
}
/*! ILU decomposition of order n
Computes ILU decomposition of order n. The matrix ILU should
be an empty matrix in row_wise creation mode. This allows the user
to either specify the number of nonzero elements or to
determine it automatically at run-time.
*/
template<class M>
void bilu_decomposition (const M& A, int n, M& ILU)
{
// iterator types
typedef typename M::ColIterator coliterator;
typedef typename M::ConstRowIterator crowiterator;
typedef typename M::ConstColIterator ccoliterator;
typedef typename M::CreateIterator createiterator;
typedef typename M::field_type K;
typedef std::map<size_t, int> map;
typedef typename map::iterator mapiterator;
// symbolic factorization phase, store generation number in first matrix element
crowiterator endi=A.end();
createiterator ci=ILU.createbegin();
for (crowiterator i=A.begin(); i!=endi; ++i)
/*! ILU decomposition of order n
Computes ILU decomposition of order n. The matrix ILU should
be an empty matrix in row_wise creation mode. This allows the user
to either specify the number of nonzero elements or to
determine it automatically at run-time.
*/
template<class M>
void blockILUDecomposition (const M& A, int n, M& ILU)
{
map rowpattern; // maps column index to generation
// iterator types
typedef typename M::ColIterator coliterator;
typedef typename M::ConstRowIterator crowiterator;
typedef typename M::ConstColIterator ccoliterator;
typedef typename M::CreateIterator createiterator;
typedef typename M::field_type K;
typedef std::map<size_t, int> map;
typedef typename map::iterator mapiterator;
// symbolic factorization phase, store generation number in first matrix element
crowiterator endi=A.end();
createiterator ci=ILU.createbegin();
for (crowiterator i=A.begin(); i!=endi; ++i)
{
map rowpattern; // maps column index to generation
// initialize pattern with row of A
for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
rowpattern[j.index()] = 0;
// initialize pattern with row of A
for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
rowpattern[j.index()] = 0;
// eliminate entries in row which are to the left of the diagonal
for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
{
if ((*ik).second<n)
// eliminate entries in row which are to the left of the diagonal
for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
{
coliterator endk = ILU[(*ik).first].end(); // end of row k
coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
for (++kj; kj!=endk; ++kj) // row k eliminates in row i
if ((*ik).second<n)
{
// we misuse the storage to store an int. If the field_type is std::complex, we have to access the real/abs part
// starting from C++11, we can use std::abs to always return a real value, even if it is double/float
using std::abs;
int generation = (int) Simd::lane(0,abs( firstmatrixelement(*kj) ));
if (generation<n)
coliterator endk = ILU[(*ik).first].end(); // end of row k
coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
for (++kj; kj!=endk; ++kj) // row k eliminates in row i
{
mapiterator ij = rowpattern.find(kj.index());
if (ij==rowpattern.end())
// we misuse the storage to store an int. If the field_type is std::complex, we have to access the real/abs part
// starting from C++11, we can use std::abs to always return a real value, even if it is double/float
using std::abs;
int generation = (int) Simd::lane(0, abs( firstMatrixElement(*kj) ));
if (generation<n)
{
rowpattern[kj.index()] = generation+1;
mapiterator ij = rowpattern.find(kj.index());
if (ij==rowpattern.end())
{
rowpattern[kj.index()] = generation+1;
}
}
}
}
}
}
// create row
for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
ci.insert((*ik).first);
++ci; // now row i exist
// create row
for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
ci.insert((*ik).first);
++ci; // now row i exist
// write generation index into entries
coliterator endILUij = ILU[i.index()].end();;
for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
Simd::lane(0,firstmatrixelement(*ILUij)) = (Simd::Scalar<K>) rowpattern[ILUij.index()];
}
// write generation index into entries
coliterator endILUij = ILU[i.index()].end();;
for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
Simd::lane(0, firstMatrixElement(*ILUij)) = (Simd::Scalar<K>) rowpattern[ILUij.index()];
}
// copy entries of A
for (crowiterator i=A.begin(); i!=endi; ++i)
{
coliterator ILUij;
coliterator endILUij = ILU[i.index()].end();;
for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
(*ILUij) = 0; // clear row
ccoliterator Aij = (*i).begin();
ccoliterator endAij = (*i).end();
ILUij = ILU[i.index()].begin();
while (Aij!=endAij && ILUij!=endILUij)
// copy entries of A
for (crowiterator i=A.begin(); i!=endi; ++i)
{
if (Aij.index()==ILUij.index())
{
*ILUij = *Aij;
++Aij; ++ILUij;
}
else
coliterator ILUij;
coliterator endILUij = ILU[i.index()].end();;
for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
(*ILUij) = 0; // clear row
ccoliterator Aij = (*i).begin();
ccoliterator endAij = (*i).end();
ILUij = ILU[i.index()].begin();
while (Aij!=endAij && ILUij!=endILUij)
{
if (Aij.index()<ILUij.index())
++Aij;
if (Aij.index()==ILUij.index())
{
*ILUij = *Aij;
++Aij; ++ILUij;
}
else
++ILUij;
{
if (Aij.index()<ILUij.index())
++Aij;
else
++ILUij;
}
}
}
}
// call decomposition on pattern
bilu0_decomposition(ILU);
}
namespace ILU {
// call decomposition on pattern
blockILU0Decomposition(ILU);
}
//! a simple compressed row storage matrix class
template <class B, class Alloc = std::allocator<B>>
struct CRS
{
......@@ -371,13 +365,17 @@ namespace Dune {
}
} // end convertToCRS
template<class CRS, class InvVector, class X, class Y>
DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use blockILUBacksolve.")
void bilu_backsolve (const CRS& lower, const CRS& upper, const InvVector& inv, X& v, const Y& d)
{ blockILUBacksolve(lower, upper, inv, v, d); }
//! LU backsolve with stored inverse in CRS format for lower and upper triangular
template<class CRS, class InvVector, class X, class Y>
void bilu_backsolve (const CRS& lower,
const CRS& upper,
const InvVector& inv,
X& v, const Y& d)
void blockILUBacksolve (const CRS& lower,
const CRS& upper,
const InvVector& inv,
X& v, const Y& d)
{
// iterator types
typedef typename Y :: block_type dblock;
......@@ -388,7 +386,7 @@ namespace Dune {
const size_type lastRow = iEnd - 1;
if( iEnd != upper.rows() )
{
DUNE_THROW(ISTLError,"ILU::bilu_backsolve: lower and upper rows must be the same");
DUNE_THROW(ISTLError,"ILU::blockILUBacksolve: lower and upper rows must be the same");
}
// lower triangular solve
......@@ -424,9 +422,24 @@ namespace Dune {
} // end namespace ILU
/** @} end documentation */
template<class M>
DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILU0Decomposition.")
void bilu0_decomposition (M& A) { ILU::blockILU0Decomposition(A); }
template<class M, class X, class Y>
DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILUBacksolve.")
void bilu_backsolve (const M& A, X& v, const Y& d) { ILU::blockILUBacksolve(A, v, d); }
template<class M>
DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::firstMatrixElement.")
decltype(auto) firstmatrixelement (M& A) { return ILU::firstMatrixElement(A); }
template<class M>
DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILUDecomposition.")
void bilu_decomposition (const M& A, int n, M& ilu) { return ILU::blockILUDecomposition(A, n, ilu); }
} // end namespace
#endif
dune/istl/ilusubdomainsolver.hh
+ 4
− 4
View file @ d20a4604
......@@ -89,7 +89,7 @@ namespace Dune {
*/
void apply (X& v, const Y& d)
{
bilu_backsolve(this->ILU,v,d);
ILU::blockILUBacksolve(this->ILU,v,d);
}
/**
* @brief Set the data of the local problem.
......@@ -121,7 +121,7 @@ namespace Dune {
*/
void apply (X& v, const Y& d)
{
bilu_backsolve(RILU,v,d);
ILU::blockILUBacksolve(RILU,v,d);
}
/**
......@@ -218,7 +218,7 @@ namespace Dune {
void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
{
this->copyToLocalMatrix(A,rowSet);
bilu0_decomposition(this->ILU);
ILU::blockILU0Decomposition(this->ILU);
}
template<class M, class X, class Y>
......@@ -228,7 +228,7 @@ namespace Dune {
std::size_t offset=copyToLocalMatrix(A,rowSet);
RILU.setSize(rowSet.size(),rowSet.size(), (1+2*offset)*rowSet.size());
RILU.setBuildMode(matrix_type::row_wise);
bilu_decomposition(this->ILU, (offset+1)/2, RILU);
ILU::blockILUDecomposition(this->ILU, (offset+1)/2, RILU);
}
/** @} */
......
dune/istl/istlexception.hh
+ 10
− 0
View file @ d20a4604
......@@ -4,6 +4,7 @@
#define DUNE_ISTL_ISTLEXCEPTION_HH
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
namespace Dune {
......@@ -51,6 +52,15 @@ namespace Dune {
*/
class SolverAbort : public ISTLError {};
//! Error when performing an operation on a matrix block
/**
* For example an error in a block LU decomposition
*/
class MatrixBlockError : public virtual Dune::FMatrixError {
public:
int r, c; // row and column index of the entry from which the error resulted
};
/** @} end documentation */
} // end namespace
......
dune/istl/preconditioners.hh
+ 4
− 4
View file @ d20a4604
......@@ -627,14 +627,14 @@ namespace Dune {
// copy A
ILU_.reset( new matrix_type( A ) );
// create ILU(0) decomposition
bilu0_decomposition( *ILU_ );
ILU::blockILU0Decomposition( *ILU_ );
}
else
{
// create matrix in build mode
ILU_.reset( new matrix_type( A.N(), A.M(), matrix_type::row_wise) );
// create ILU(n) decomposition
bilu_decomposition( A, n, *ILU_ );
ILU::blockILUDecomposition( A, n, *ILU_ );
}
if( resort )
......@@ -665,11 +665,11 @@ namespace Dune {
{
if( ILU_ )
{
bilu_backsolve( *ILU_, v, d);
ILU::blockILUBacksolve( *ILU_, v, d);
}
else
{
ILU::bilu_backsolve(lower_, upper_, inv_, v, d);
ILU::blockILUBacksolve(lower_, upper_, inv_, v, d);
}
if( wNotIdentity_ )
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment