Merge pull request #1218 from blattms/faster-ilu-old-order

This is the faster ilu pull request but with the traversal of the columns of upper in the original ordering

opm/autodiff/ISTLSolver.hpp

@@ -289,13 +289,14 @@ namespace Opm
             }
         }
 
-        typedef Dune::SeqILU0<Matrix, Vector, Vector> SeqPreconditioner;
+        typedef ParallelOverlappingILU0<Matrix,Vector,Vector> SeqPreconditioner;
 
         template <class Operator>
         std::unique_ptr<SeqPreconditioner> constructPrecond(Operator& opA, const Dune::Amg::SequentialInformation&) const
         {
-            const double relax = 0.9;
+            const double relax   = parameters_.ilu_relaxation_;
-            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), relax));
+            const int ilu_fillin = parameters_.ilu_fillin_level_;
+            std::unique_ptr<SeqPreconditioner> precond(new SeqPreconditioner(opA.getmat(), ilu_fillin, relax));
             return precond;
         }
 
@@ -307,7 +308,7 @@ namespace Opm
         constructPrecond(Operator& opA, const Comm& comm) const
         {
             typedef std::unique_ptr<ParPreconditioner> Pointer;
-            const double relax = 0.9;
+            const double relax  = parameters_.ilu_relaxation_;
             return Pointer(new ParPreconditioner(opA.getmat(), comm, relax));
         }
 #endif

opm/autodiff/NewtonIterationBlackoilInterleaved.hpp

@@ -37,9 +37,11 @@ namespace Opm
     struct NewtonIterationBlackoilInterleavedParameters
     {
         double linear_solver_reduction_;
+        double ilu_relaxation_;
         int    linear_solver_maxiter_;
         int    linear_solver_restart_;
         int    linear_solver_verbosity_;
+        int    ilu_fillin_level_;
         bool   newton_use_gmres_;
         bool   require_full_sparsity_pattern_;
         bool   ignoreConvergenceFailure_;
@@ -61,6 +63,8 @@ namespace Opm
             require_full_sparsity_pattern_ = param.getDefault("require_full_sparsity_pattern", require_full_sparsity_pattern_);
             ignoreConvergenceFailure_ = param.getDefault("linear_solver_ignoreconvergencefailure", ignoreConvergenceFailure_);
             linear_solver_use_amg_    = param.getDefault("linear_solver_use_amg", linear_solver_use_amg_ );
+            ilu_relaxation_           = param.getDefault("ilu_relaxation", ilu_relaxation_ );
+            ilu_fillin_level_         = param.getDefault("ilu_fillin_level",  ilu_fillin_level_ );
         }
 
         // set default values
@@ -74,6 +78,8 @@ namespace Opm
             require_full_sparsity_pattern_ = false;
             ignoreConvergenceFailure_ = false;
             linear_solver_use_amg_    = false;
+            ilu_fillin_level_         = 0;
+            ilu_relaxation_           = 0.9;
         }
     };
 

opm/autodiff/ParallelOverlappingILU0.hpp

@@ -20,13 +20,18 @@
 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
 
+#include <opm/common/Exceptions.hpp>
+
 #include <dune/istl/preconditioner.hh>
 #include <dune/istl/paamg/smoother.hh>
+#include <dune/istl/paamg/pinfo.hh>
 
 namespace Opm
 {
 
-template<class M, class X, class Y, class C>
+//template<class M, class X, class Y, class C>
+//class ParallelOverlappingILU0;
+template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
 class ParallelOverlappingILU0;
 
 } // end namespace Opm
@@ -65,10 +70,79 @@ struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,Paral
 
 } // end namespace Amg
 
+
+
 } // end namespace Dune
 
 namespace Opm
 {
+    namespace detail
+    {
+      //! compute ILU decomposition of A. A is overwritten by its decomposition
+      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() );
+
+        // implement left looking variant with stored inverse
+        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() );
+
+          // eliminate entries left of diagonal; store L factor
+          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 size_type iIndex = i.index();
+          upper.reserveAdditional( (*i).size() );
+
+          // store in reverse row order
+          // eliminate entries left of diagonal; store L factor
+          for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
+          {
+            const size_type jIndex = j.index();
+            if( j.index() == iIndex )
+            {
+              inv[ row ] = (*j);
+	      break;
+            }
+            else if ( j.index() >= i.index() )
+            {
+              upper.push_back( (*j), jIndex );
+              ++colcount ;
+            }
+          }
+          upper.rows_[ row+1 ] = colcount;
+        }
+      }
+    } // end namespace detail
 
 /// \brief A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as
 ///
@@ -85,9 +159,12 @@ namespace Opm
 /// \tparam Range The type of the Vector representing the range.
 /// \tparam ParallelInfo The type of the parallel information object
 ///         used, e.g. Dune::OwnerOverlapCommunication
-template<class Matrix, class Domain, class Range, class ParallelInfo>
+template<class Matrix, class Domain, class Range, class ParallelInfoT>
 class ParallelOverlappingILU0
-    : public Dune::Preconditioner<Domain,Range> {
+    : public Dune::Preconditioner<Domain,Range>
+{
+    typedef ParallelInfoT ParallelInfo;
+
 public:
     //! \brief The matrix type the preconditioner is for.
     typedef typename Dune::remove_const<Matrix>::type matrix_type;
@@ -98,42 +175,105 @@ public:
     //! \brief The field type of the preconditioner.
     typedef typename Domain::field_type field_type;
 
+    typedef typename matrix_type::block_type  block_type;
+    typedef typename matrix_type::size_type   size_type;
+
+protected:
+    struct CRS
+    {
+      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_;
+    };
+
+public:
     // define the category
     enum {
         //! \brief The category the preconditioner is part of.
-        category=Dune::SolverCategory::overlapping
+        category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
+            Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
     };
 
+    /*! \brief Constructor.
+
+      Constructor gets all parameters to operate the prec.
+      \param A The matrix to operate on.
+      \param n ILU fill in level (for testing). This does not work in parallel.
+      \param w The relaxation factor.
+    */
+    ParallelOverlappingILU0 (const Matrix& A, const int n, const field_type w )
+        : lower_(),
+          upper_(),
+          inv_(),
+          comm_(nullptr), w_(w),
+          relaxation_( std::abs( w - 1.0 ) > 1e-15 )
+    {
+        init( A, n );
+    }
+
     /*! \brief Constructor.
 
       Constructor gets all parameters to operate the prec.
       \param A The matrix to operate on.
       \param w The relaxation factor.
     */
-    ParallelOverlappingILU0 (const Matrix& A, const ParallelInfo& comm,
+    ParallelOverlappingILU0 (const Matrix& A, const field_type w)
-                             field_type w)
+        : ParallelOverlappingILU0( A, 0, w )
-        : ilu_(A), comm_(comm), w_(w)
     {
-        int ilu_setup_successful = 1;
+    }
-        std::string message;
+
-        try
+    /*! \brief Constructor.
-        {
+
-            bilu0_decomposition(ilu_);
+      Constructor gets all parameters to operate the prec.
-        }
+      \param A      The matrix to operate on.
-        catch ( Dune::MatrixBlockError error )
+      \param comm   communication object, e.g. Dune::OwnerOverlapCopyCommunication
-        {
+      \param w      The relaxation factor.
-            message = error.what();
+    */
-            std::cerr<<"Exception occured on process " <<
+    ParallelOverlappingILU0 (const Matrix& A, const ParallelInfo& comm, const field_type w)
-                comm.communicator().rank() << " during " <<
+        : lower_(),
-                "setup of ILU0 preconditioner with message: " <<
+          upper_(),
-                message<<std::endl;
+          inv_(),
-            ilu_setup_successful = 0;
+          comm_(&comm), w_(w),
-        }
+          relaxation_( std::abs( w - 1.0 ) > 1e-15 )
-        // Check whether there was a problem on some process
+    {
-        if ( comm.communicator().min(ilu_setup_successful) == 0 )
+        init( A, 0 );
-        {
-            throw Dune::MatrixBlockError();
-        }
     }
 
     /*!
@@ -155,32 +295,66 @@ public:
     virtual void apply (Domain& v, const Range& d)
     {
         Range& md = const_cast<Range&>(d);
-        comm_.copyOwnerToAll(md,md);
+        copyOwnerToAll( md );
-        auto endrow=ilu_.end();
+
-        for ( auto row = ilu_.begin(); row != endrow; ++row )
+        // iterator types
+        typedef typename Range ::block_type  dblock;
+        typedef typename Domain::block_type  vblock;
+
+        const size_type iEnd = lower_.rows();
+        const size_type lastRow = iEnd - 1;
+        if( iEnd != upper_.rows() )
         {
-            auto rhs(d[row.index()]);
+            std::abort();
-            for ( auto col = row->begin(); col.index() < row.index(); ++col )
+           // OPM_THROW(std::logic_error,"ILU: lower and upper rows must be the same");
-            {
-                col->mmv(v[col.index()],rhs);
-            }
-            v[row.index()] = rhs;
         }
-        comm_.copyOwnerToAll(v, v);
+
-        auto rendrow = ilu_.beforeBegin();
+        // lower triangular solve
-        for( auto row = ilu_.beforeEnd(); row != rendrow; --row)
+        for( size_type i=0; i<iEnd; ++ i )
         {
-            auto rhs(v[row.index()]);
+          dblock rhs( d[ i ] );
-            auto col = row->beforeEnd();
+          const size_type rowI     = lower_.rows_[ i ];
-            for( ; col.index() > row.index(); --col)
+          const size_type rowINext = lower_.rows_[ i+1 ];
-            {
+
-                col->mmv(v[col.index()], rhs);
+          for( size_type col = rowI; col < rowINext; ++ col )
-            }
+          {
-            v[row.index()] = 0;
+            lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
-            col->umv(rhs, v[row.index()]);
+          }
+
+          v[ i ] = rhs;  // Lii = I
+        }
+
+        copyOwnerToAll( v );
+
+        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);
+        }
+
+        copyOwnerToAll( v );
+
+        if( relaxation_ ) {
+            v *= w_;
+        }
+    }
+
+    template <class V>
+    void copyOwnerToAll( V& v ) const
+    {
+        if( comm_ ) {
+            comm_->copyOwnerToAll(v, v);
         }
-        comm_.copyOwnerToAll(v, v);
-        v *= w_;
     }
 
     /*!
@@ -193,12 +367,57 @@ public:
         DUNE_UNUSED_PARAMETER(x);
     }
 
-private:
+protected:
+    void init( const Matrix& A, const int iluIteration )
+    {
+        int ilu_setup_successful = 1;
+        std::string message;
+        const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
+
+        std::unique_ptr< Matrix > ILU;
+
+        try
+        {
+            if( iluIteration == 0 ) {
+                // create ILU-0 decomposition
+                ILU.reset( new Matrix( A ) );
+                bilu0_decomposition( *ILU );
+            }
+            else {
+                // create ILU-n decomposition
+                ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
+                bilu_decomposition( A, iluIteration, *ILU );
+            }
+        }
+        catch ( Dune::MatrixBlockError error )
+        {
+            message = error.what();
+            std::cerr<<"Exception occured on process " << rank << " during " <<
+                "setup of ILU0 preconditioner with message: " <<
+                message<<std::endl;
+            ilu_setup_successful = 0;
+        }
+
+        // Check whether there was a problem on some process
+        if ( comm_ && comm_->communicator().min(ilu_setup_successful) == 0 )
+        {
+            throw Dune::MatrixBlockError();
+        }
+
+        // store ILU in simple CRS format
+        detail::convertToCRS( *ILU, lower_, upper_, inv_ );
+    }
+
+protected:
     //! \brief The ILU0 decomposition of the matrix.
-    Matrix ilu_;
+    CRS lower_;
-    const ParallelInfo& comm_;
+    CRS upper_;
+    std::vector< block_type > inv_;
+
+    const ParallelInfo* comm_;
     //! \brief The relaxation factor to use.
-    field_type w_;
+    const field_type w_;
+    const bool relaxation_;
 
 };