From be38b8de9b9af39c01fde9b3c89f647a498c2037 Mon Sep 17 00:00:00 2001
From: Robert K <robertk@posteo.org>
Date: Tue, 30 May 2017 16:13:20 +0200
Subject: [PATCH] [feature][SeqILU] faster implementation of ILU
preconditioners using a simple Compressed Row Storage for the lower and upper
triangular matrices. Also, SeqILU unifies both, SeqILU0 and SeqILUn.
---
dune/istl/ilu.hh | 176 ++++++++++++++++++++++++++++++++-
dune/istl/paamg/smoother.hh | 48 ++++++++-
dune/istl/preconditioners.hh | 129 ++++++++++++++++++++++++
dune/istl/test/multirhstest.cc | 5 +
4 files changed, 354 insertions(+), 4 deletions(-)
diff --git a/dune/istl/ilu.hh b/dune/istl/ilu.hh
index 4b3f19aff..7f9b1003b 100644
--- a/dune/istl/ilu.hh
+++ b/dune/istl/ilu.hh
@@ -10,6 +10,7 @@
#include <string>
#include <set>
#include <map>
+#include <vector>
#include <dune/common/fmatrix.hh>
#include "istlexception.hh"
@@ -169,8 +170,8 @@ namespace Dune {
createiterator ci=ILU.createbegin();
for (crowiterator i=A.begin(); i!=endi; ++i)
{
- // std::cout << "in row " << i.index() << std::endl;
- map rowpattern; // maps column index to generation
+ // std::cout << "in row " << i.index() << std::endl;
+ map rowpattern; // maps column index to generation
// initialize pattern with row of A
for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
@@ -218,7 +219,7 @@ namespace Dune {
firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
}
- // printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
+ // printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
// copy entries of A
for (crowiterator i=A.begin(); i!=endi; ++i)
@@ -251,6 +252,175 @@ namespace Dune {
bilu0_decomposition(ILU);
}
+ namespace ILU {
+
+ template <class B>
+ struct CRS
+ {
+ typedef B block_type;
+ typedef size_t size_type;
+
+ CRS() : nRows_( 0 ) {}
+
+ size_type rows() const { return nRows_; }
+
+ size_type nonZeros() const
+ {
+ assert( rows_[ rows() ] != size_type(-1) );
+ return rows_[ rows() ];
+ }
+
+ void resize( const size_type nRows )
+ {
+ if( nRows_ != nRows )
+ {
+ nRows_ = nRows ;
+ rows_.resize( nRows_+1, size_type(-1) );
+ }
+ }
+
+ void reserveAdditional( const size_type nonZeros )
+ {
+ const size_type needed = values_.size() + nonZeros ;
+ if( values_.capacity() < needed )
+ {
+ const size_type estimate = needed * 1.1;
+ values_.reserve( estimate );
+ cols_.reserve( estimate );
+ }
+ }
+
+ void push_back( const block_type& value, const size_type index )
+ {
+ values_.push_back( value );
+ cols_.push_back( index );
+ }
+
+ std::vector< size_type > rows_;
+ std::vector< block_type > values_;
+ std::vector< size_type > cols_;
+ size_type nRows_;
+ };
+
+ //! convert ILU decomposition into CRS format for lower and upper triangular and inverse.
+ template<class M, class CRS, class InvVector>
+ void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
+ {
+ typedef typename M :: size_type size_type;
+
+ lower.resize( A.N() );
+ upper.resize( A.N() );
+ inv.resize( A.N() );
+
+ lower.reserveAdditional( 2*A.N() );
+
+ const auto endi = A.end();
+ size_type row = 0;
+ size_type colcount = 0;
+ lower.rows_[ 0 ] = colcount;
+ for (auto i=A.begin(); i!=endi; ++i, ++row)
+ {
+ const size_type iIndex = i.index();
+ lower.reserveAdditional( (*i).size() );
+
+ // store entries left of diagonal
+ for (auto j=(*i).begin(); j.index() < iIndex; ++j )
+ {
+ lower.push_back( (*j), j.index() );
+ ++colcount;
+ }
+ lower.rows_[ iIndex+1 ] = colcount;
+ }
+
+ const auto rendi = A.beforeBegin();
+ row = 0;
+ colcount = 0;
+ upper.rows_[ 0 ] = colcount ;
+
+ upper.reserveAdditional( lower.nonZeros() + A.N() );
+
+ // NOTE: upper and inv store entries in reverse order, reverse here
+ // relative to ILU
+ for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
+ {
+ const auto endij=(*i).end(); // end of row i
+
+ const size_type iIndex = i.index();
+ upper.reserveAdditional( (*i).size() );
+
+ // store in reverse row order for faster access during backsolve
+ for (auto j=(*i).begin(); j != endij; ++j )
+ {
+ const size_type jIndex = j.index();
+ if( j.index() == iIndex )
+ {
+ inv[ row ] = (*j);
+ }
+ else if ( j.index() >= i.index() )
+ {
+ upper.push_back( (*j), jIndex );
+ ++colcount ;
+ }
+ }
+ upper.rows_[ row+1 ] = colcount;
+ }
+ } // end convertToCRS
+
+
+ //! LU backsolve with stored inverse in CRS format for lower and upper triangular
+ template<class CRS, class mblock, class X, class Y>
+ void bilu_backsolve (const CRS& lower,
+ const CRS& upper,
+ const std::vector< mblock >& inv,
+ X& v, const Y& d)
+ {
+ // iterator types
+ typedef typename Y :: block_type dblock;
+ typedef typename X :: block_type vblock;
+ typedef typename X :: size_type size_type ;
+
+ const size_type iEnd = lower.rows();
+ const size_type lastRow = iEnd - 1;
+ if( iEnd != upper.rows() )
+ {
+ DUNE_THROW(ISTLError,"ILU::bilu_backsolve: lower and upper rows must be the same");
+ }
+
+ // lower triangular solve
+ for( size_type i=0; i<iEnd; ++ i )
+ {
+ dblock rhs( d[ i ] );
+ const size_type rowI = lower.rows_[ i ];
+ const size_type rowINext = lower.rows_[ i+1 ];
+
+ for( size_type col = rowI; col < rowINext; ++ col )
+ {
+ lower.values_[ col ].mmv( v[ lower.cols_[ col ] ], rhs );
+ }
+
+ v[ i ] = rhs; // Lii = I
+ }
+
+ // upper triangular solve
+ for( size_type i=0; i<iEnd; ++ i )
+ {
+ vblock& vBlock = v[ lastRow - i ];
+ vblock rhs ( vBlock );
+ const size_type rowI = upper.rows_[ i ];
+ const size_type rowINext = upper.rows_[ i+1 ];
+
+ for( size_type col = rowI; col < rowINext; ++ col )
+ {
+ upper.values_[ col ].mmv( v[ upper.cols_[ col ] ], rhs );
+ }
+
+ // apply inverse and store result
+ inv[ i ].mv( rhs, vBlock);
+ }
+ }
+
+ } // end namespace ILU
+
/** @} end documentation */
diff --git a/dune/istl/paamg/smoother.hh b/dune/istl/paamg/smoother.hh
index 1f689ed64..9debb1dea 100644
--- a/dune/istl/paamg/smoother.hh
+++ b/dune/istl/paamg/smoother.hh
@@ -279,7 +279,7 @@ namespace Dune
/**
- * @brief Policy for the construction of the SeqJac smoother
+ * @brief Policy for the construction of the SeqILUn smoother
*/
template<class M, class X, class Y>
struct ConstructionTraits<SeqILUn<M,X,Y> >
@@ -299,6 +299,52 @@ namespace Dune
};
+
+ template<class M, class X, class Y>
+ class ConstructionArgs<SeqILU<M,X,Y> >
+ : public DefaultConstructionArgs<SeqILU<M,X,Y> >
+ {
+ public:
+ ConstructionArgs(int n=0)
+ : n_(n)
+ {}
+
+ void setN(int n)
+ {
+ n_ = n;
+ }
+
+ int getN()
+ {
+ return n_;
+ }
+
+ private:
+ int n_;
+ };
+
+
+ /**
+ * @brief Policy for the construction of the SeqILU smoother
+ */
+ template<class M, class X, class Y>
+ struct ConstructionTraits<SeqILU<M,X,Y> >
+ {
+ typedef ConstructionArgs<SeqILU<M,X,Y> > Arguments;
+
+ static inline SeqILU<M,X,Y>* construct(Arguments& args)
+ {
+ return new SeqILU<M,X,Y>(args.getMatrix(), args.getN(),
+ args.getArgs().relaxationFactor);
+ }
+
+ static void deconstruct(SeqILU<M,X,Y>* ilu)
+ {
+ delete ilu;
+ }
+
+ };
+
/**
* @brief Policy for the construction of the ParSSOR smoother
*/
diff --git a/dune/istl/preconditioners.hh b/dune/istl/preconditioners.hh
index 84acf024c..6e883c9ba 100644
--- a/dune/istl/preconditioners.hh
+++ b/dune/istl/preconditioners.hh
@@ -7,6 +7,7 @@
#include <complex>
#include <iostream>
#include <iomanip>
+#include <memory>
#include <string>
#include <dune/common/unused.hh>
@@ -487,6 +488,134 @@ namespace Dune {
+ /*!
+ \brief Sequential ILU preconditioner.
+
+ Wraps the naked ISTL generic ILU preconditioner into the solver framework.
+
+ \tparam M The matrix type to operate on
+ \tparam X Type of the update
+ \tparam Y Type of the defect
+ \tparam l Ignored. Just there to have the same number of template arguments
+ as other preconditioners.
+ */
+ template<class M, class X, class Y, int l=1>
+ class SeqILU : public Preconditioner<X,Y> {
+ public:
+ //! \brief The matrix type the preconditioner is for.
+ typedef typename std::remove_const<M>::type matrix_type;
+ //! block type of matrix
+ typedef typename matrix_type :: block_type block_type;
+ //! \brief The domain type of the preconditioner.
+ typedef X domain_type;
+ //! \brief The range type of the preconditioner.
+ typedef Y range_type;
+ //! \brief The field type of the preconditioner.
+ typedef typename X::field_type field_type;
+
+ //! \brief type of ILU storage
+ typedef typename ILU::CRS< block_type > CRS;
+
+ /*! \brief Constructor.
+
+ Constructor invoking ILU(0) gets all parameters to operate the prec.
+ \param A The matrix to operate on.
+ \param w The relaxation factor.
+ */
+ SeqILU (const M& A, field_type w)
+ : SeqILU( A, 0, w ) // construct ILU(0)
+ {
+ }
+
+ /*! \brief Constructor.
+
+ Constructor invoking ILU(n).
+ \param A The matrix to operate on.
+ \param n The number of iterations to perform.
+ \param w The relaxation factor.
+ */
+ SeqILU (const M& A, int n, field_type w)
+ : lower_(),
+ upper_(),
+ inv_(),
+ w_(w),
+ wNotIdentity_( std::abs( w_ - field_type(1) ) > 1e-15 )
+ {
+ std::unique_ptr< matrix_type > ILUA;
+
+ if( n == 0 )
+ {
+ // copy A
+ ILUA.reset( new matrix_type( A ) );
+ // create ILU(0) decomposition
+ bilu0_decomposition( *ILUA );
+ }
+ else
+ {
+ // create matrix in build mode
+ ILUA.reset( new matrix_type( A.N(), A.M(), matrix_type::row_wise) );
+ // create ILU(n) decomposition
+ bilu_decomposition( A, n, *ILUA );
+ }
+
+ // store ILU in simple CRS format
+ ILU::convertToCRS( *ILUA, lower_, upper_, inv_ );
+ }
+
+ /*!
+ \brief Prepare the preconditioner.
+
+ \copydoc Preconditioner::pre(X&,Y&)
+ */
+ virtual void pre (X& x, Y& b)
+ {
+ DUNE_UNUSED_PARAMETER(x);
+ DUNE_UNUSED_PARAMETER(b);
+ }
+
+ /*!
+ \brief Apply the preconditoner.
+
+ \copydoc Preconditioner::apply(X&,const Y&)
+ */
+ virtual void apply (X& v, const Y& d)
+ {
+ ILU::bilu_backsolve(lower_, upper_, inv_, v, d);
+ if( wNotIdentity_ )
+ {
+ v *= w_;
+ }
+ }
+
+ /*!
+ \brief Clean up.
+
+ \copydoc Preconditioner::post(X&)
+ */
+ virtual void post (X& x)
+ {
+ DUNE_UNUSED_PARAMETER(x);
+ }
+
+ //! Category of the preconditioner (see SolverCategory::Category)
+ virtual SolverCategory::Category category() const
+ {
+ return SolverCategory::sequential;
+ }
+
+ protected:
+ //! \brief The ILU(n) decomposition of the matrix.
+ CRS lower_;
+ CRS upper_;
+ std::vector< block_type > inv_;
+
+ //! \brief The relaxation factor to use.
+ const field_type w_;
+ //! \brief true if w != 1.0
+ const bool wNotIdentity_;
+ };
+
+
/*!
\brief Sequential ILU0 preconditioner.
diff --git a/dune/istl/test/multirhstest.cc b/dune/istl/test/multirhstest.cc
index da3a66670..3aca3f274 100644
--- a/dune/istl/test/multirhstest.cc
+++ b/dune/istl/test/multirhstest.cc
@@ -142,6 +142,9 @@ void test_all(unsigned int Runs = 1)
Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,0.1); // preconditioner object
Dune::SeqILUn<Matrix,Vector,Vector> ilu1(A,1,0.1); // preconditioner object
+ Dune::SeqILU<Matrix,Vector,Vector> ilu_0(A,0.1); // preconditioner object
+ Dune::SeqILU<Matrix,Vector,Vector> ilu_1(A,1,0.1); // preconditioner object
+
// AMG
typedef Dune::Amg::RowSum Norm;
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Matrix,Norm> >
@@ -173,6 +176,8 @@ void test_all(unsigned int Runs = 1)
test_all_solvers("SSOR", op,ssor,N,Runs);
test_all_solvers("ILU0", op,ilu0,N,Runs);
test_all_solvers("ILU1", op,ilu1,N,Runs);
+ test_all_solvers("ILU(0)", op,ilu_0,N,Runs);
+ test_all_solvers("ILU(1)", op,ilu_1,N,Runs);
test_all_solvers("AMG", op,amg,N,Runs);
}
--
GitLab