[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

@@ -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 */
 

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
      */

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.
 

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);
 }