All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
ParallelOverlappingILU0.hpp
1 /*
2  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
3  Copyright 2015 Statoil AS
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 #ifndef OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
21 #define OPM_PARALLELOVERLAPPINGILU0_HEADER_INCLUDED
22 
23 #include <opm/common/Exceptions.hpp>
24 
25 #include <dune/istl/preconditioner.hh>
26 #include <dune/istl/paamg/smoother.hh>
27 #include <dune/istl/paamg/pinfo.hh>
28 
29 namespace Opm
30 {
31 
32 //template<class M, class X, class Y, class C>
33 //class ParallelOverlappingILU0;
34 template<class Matrix, class Domain, class Range, class ParallelInfo = Dune::Amg::SequentialInformation>
36 
37 } // end namespace Opm
38 
39 namespace Dune
40 {
41 
42 namespace Amg
43 {
44 
51 template<class Matrix, class Domain, class Range, class ParallelInfo>
52 struct ConstructionTraits<Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo> >
53 {
54  typedef Dune::SeqILU0<Matrix,Domain,Range> T;
55  typedef DefaultParallelConstructionArgs<T,ParallelInfo> Arguments;
56  typedef ConstructionTraits<T> SeqConstructionTraits;
57  static inline Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* construct(Arguments& args)
58  {
60  args.getComm(),
61  args.getArgs().relaxationFactor);
62  }
63 
64  static inline void deconstruct(Opm::ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfo>* bp)
65  {
66  delete bp;
67  }
68 
69 };
70 
71 } // end namespace Amg
72 
73 
74 
75 } // end namespace Dune
76 
77 namespace Opm
78 {
79  namespace detail
80  {
82  template<class M, class CRS, class InvVector>
83  void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
84  {
85  typedef typename M :: size_type size_type;
86 
87  lower.resize( A.N() );
88  upper.resize( A.N() );
89  inv.resize( A.N() );
90 
91  lower.reserveAdditional( 2*A.N() );
92 
93  // implement left looking variant with stored inverse
94  const auto endi = A.end();
95  size_type row = 0;
96  size_type colcount = 0;
97  lower.rows_[ 0 ] = colcount;
98  for (auto i=A.begin(); i!=endi; ++i, ++row)
99  {
100  const size_type iIndex = i.index();
101  lower.reserveAdditional( (*i).size() );
102 
103  // eliminate entries left of diagonal; store L factor
104  for (auto j=(*i).begin(); j.index() < iIndex; ++j )
105  {
106  lower.push_back( (*j), j.index() );
107  ++colcount;
108  }
109  lower.rows_[ iIndex+1 ] = colcount;
110  }
111 
112  const auto rendi = A.beforeBegin();
113  row = 0;
114  colcount = 0;
115  upper.rows_[ 0 ] = colcount ;
116 
117  upper.reserveAdditional( lower.nonZeros() + A.N() );
118 
119  // NOTE: upper and inv store entries in reverse order, reverse here
120  // relative to ILU
121  for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
122  {
123  const size_type iIndex = i.index();
124  upper.reserveAdditional( (*i).size() );
125 
126  // store in reverse row order
127  // eliminate entries left of diagonal; store L factor
128  for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
129  {
130  const size_type jIndex = j.index();
131  if( j.index() == iIndex )
132  {
133  inv[ row ] = (*j);
134  break;
135  }
136  else if ( j.index() >= i.index() )
137  {
138  upper.push_back( (*j), jIndex );
139  ++colcount ;
140  }
141  }
142  upper.rows_[ row+1 ] = colcount;
143  }
144  }
145  } // end namespace detail
146 
162 template<class Matrix, class Domain, class Range, class ParallelInfoT>
163 class ParallelOverlappingILU0
164  : public Dune::Preconditioner<Domain,Range>
165 {
166  typedef ParallelInfoT ParallelInfo;
167 
168 
169 public:
171  typedef typename std::remove_const<Matrix>::type matrix_type;
173  typedef Domain domain_type;
175  typedef Range range_type;
177  typedef typename Domain::field_type field_type;
178 
179  typedef typename matrix_type::block_type block_type;
180  typedef typename matrix_type::size_type size_type;
181 
182 protected:
183  struct CRS
184  {
185  CRS() : nRows_( 0 ) {}
186 
187  size_type rows() const { return nRows_; }
188 
189  size_type nonZeros() const
190  {
191  assert( rows_[ rows() ] != size_type(-1) );
192  return rows_[ rows() ];
193  }
194 
195  void resize( const size_type nRows )
196  {
197  if( nRows_ != nRows )
198  {
199  nRows_ = nRows ;
200  rows_.resize( nRows_+1, size_type(-1) );
201  }
202  }
203 
204  void reserveAdditional( const size_type nonZeros )
205  {
206  const size_type needed = values_.size() + nonZeros ;
207  if( values_.capacity() < needed )
208  {
209  const size_type estimate = needed * 1.1;
210  values_.reserve( estimate );
211  cols_.reserve( estimate );
212  }
213  }
214 
215  void push_back( const block_type& value, const size_type index )
216  {
217  values_.push_back( value );
218  cols_.push_back( index );
219  }
220 
221  std::vector< size_type > rows_;
222  std::vector< block_type > values_;
223  std::vector< size_type > cols_;
224  size_type nRows_;
225  };
226 
227 public:
228  // define the category
229  enum {
231  category = std::is_same<ParallelInfoT, Dune::Amg::SequentialInformation>::value ?
232  Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping
233  };
234 
242  template<class BlockType, class Alloc>
243  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
244  const int n, const field_type w )
245  : lower_(),
246  upper_(),
247  inv_(),
248  comm_(nullptr), w_(w),
249  relaxation_( std::abs( w - 1.0 ) > 1e-15 )
250  {
251  // BlockMatrix is a Subclass of FieldMatrix that just adds
252  // methods. Therefore this cast should be safe.
253  init( reinterpret_cast<const Matrix&>(A), n );
254  }
255 
262  template<class BlockType, class Alloc>
263  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
264  const field_type w)
265  : ParallelOverlappingILU0( A, 0, w )
266  {
267  }
268 
276  template<class BlockType, class Alloc>
277  ParallelOverlappingILU0 (const Dune::BCRSMatrix<BlockType,Alloc>& A,
278  const ParallelInfo& comm, const field_type w)
279  : lower_(),
280  upper_(),
281  inv_(),
282  comm_(&comm), w_(w),
283  relaxation_( std::abs( w - 1.0 ) > 1e-15 )
284  {
285  // BlockMatrix is a Subclass of FieldMatrix that just adds
286  // methods. Therefore this cast should be safe.
287  init( reinterpret_cast<const Matrix&>(A), 0 );
288  }
289 
295  virtual void pre (Domain& x, Range& b)
296  {
297  DUNE_UNUSED_PARAMETER(x);
298  DUNE_UNUSED_PARAMETER(b);
299  }
300 
306  virtual void apply (Domain& v, const Range& d)
307  {
308  Range& md = const_cast<Range&>(d);
309  copyOwnerToAll( md );
310 
311  // iterator types
312  typedef typename Range ::block_type dblock;
313  typedef typename Domain::block_type vblock;
314 
315  const size_type iEnd = lower_.rows();
316  const size_type lastRow = iEnd - 1;
317  if( iEnd != upper_.rows() )
318  {
319  std::abort();
320  // OPM_THROW(std::logic_error,"ILU: lower and upper rows must be the same");
321  }
322 
323  // lower triangular solve
324  for( size_type i=0; i<iEnd; ++ i )
325  {
326  dblock rhs( d[ i ] );
327  const size_type rowI = lower_.rows_[ i ];
328  const size_type rowINext = lower_.rows_[ i+1 ];
329 
330  for( size_type col = rowI; col < rowINext; ++ col )
331  {
332  lower_.values_[ col ].mmv( v[ lower_.cols_[ col ] ], rhs );
333  }
334 
335  v[ i ] = rhs; // Lii = I
336  }
337 
338  copyOwnerToAll( v );
339 
340  for( size_type i=0; i<iEnd; ++ i )
341  {
342  vblock& vBlock = v[ lastRow - i ];
343  vblock rhs ( vBlock );
344  const size_type rowI = upper_.rows_[ i ];
345  const size_type rowINext = upper_.rows_[ i+1 ];
346 
347  for( size_type col = rowI; col < rowINext; ++ col )
348  {
349  upper_.values_[ col ].mmv( v[ upper_.cols_[ col ] ], rhs );
350  }
351 
352  // apply inverse and store result
353  inv_[ i ].mv( rhs, vBlock);
354  }
355 
356  copyOwnerToAll( v );
357 
358  if( relaxation_ ) {
359  v *= w_;
360  }
361  }
362 
363  template <class V>
364  void copyOwnerToAll( V& v ) const
365  {
366  if( comm_ ) {
367  comm_->copyOwnerToAll(v, v);
368  }
369  }
370 
376  virtual void post (Range& x)
377  {
378  DUNE_UNUSED_PARAMETER(x);
379  }
380 
381 protected:
382  void init( const Matrix& A, const int iluIteration )
383  {
384  int ilu_setup_successful = 1;
385  std::string message;
386  const int rank = ( comm_ ) ? comm_->communicator().rank() : 0;
387 
388  std::unique_ptr< Matrix > ILU;
389 
390  try
391  {
392  if( iluIteration == 0 ) {
393  // create ILU-0 decomposition
394  ILU.reset( new Matrix( A ) );
395  bilu0_decomposition( *ILU );
396  }
397  else {
398  // create ILU-n decomposition
399  ILU.reset( new Matrix( A.N(), A.M(), Matrix::row_wise) );
400  bilu_decomposition( A, iluIteration, *ILU );
401  }
402  }
403  catch ( Dune::MatrixBlockError error )
404  {
405  message = error.what();
406  std::cerr<<"Exception occured on process " << rank << " during " <<
407  "setup of ILU0 preconditioner with message: " <<
408  message<<std::endl;
409  ilu_setup_successful = 0;
410  }
411 
412  // Check whether there was a problem on some process
413  if ( comm_ && comm_->communicator().min(ilu_setup_successful) == 0 )
414  {
415  throw Dune::MatrixBlockError();
416  }
417 
418  // store ILU in simple CRS format
419  detail::convertToCRS( *ILU, lower_, upper_, inv_ );
420  }
421 
422 protected:
425  CRS upper_;
426  std::vector< block_type > inv_;
427 
428  const ParallelInfo* comm_;
430  const field_type w_;
431  const bool relaxation_;
432 
433 };
434 
435 } // end namespace Opm
436 #endif
virtual void pre(Domain &x, Range &b)
Prepare the preconditioner.
Definition: ParallelOverlappingILU0.hpp:295
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:35
Range range_type
The range type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:175
const field_type w_
The relaxation factor to use.
Definition: ParallelOverlappingILU0.hpp:430
CRS lower_
The ILU0 decomposition of the matrix.
Definition: ParallelOverlappingILU0.hpp:424
Domain domain_type
The domain type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:173
std::remove_const< Matrix >::type matrix_type
The matrix type the preconditioner is for.
Definition: ParallelOverlappingILU0.hpp:171
The category the preconditioner is part of.
Definition: ParallelOverlappingILU0.hpp:231
virtual void apply(Domain &v, const Range &d)
Apply the preconditoner.
Definition: ParallelOverlappingILU0.hpp:306
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const int n, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:243
Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:177
Definition: ParallelOverlappingILU0.hpp:183
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:263
virtual void post(Range &x)
Clean up.
Definition: ParallelOverlappingILU0.hpp:376
ParallelOverlappingILU0(const Dune::BCRSMatrix< BlockType, Alloc > &A, const ParallelInfo &comm, const field_type w)
Constructor.
Definition: ParallelOverlappingILU0.hpp:277