My Project
ParallelOverlappingILU0_impl.hpp
1 /*
2  Copyright 2015, 2022 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 
21 #include <opm/simulators/linalg/ParallelOverlappingILU0.hpp>
22 
23 #include <dune/istl/ilu.hh>
24 #include <dune/istl/owneroverlapcopy.hh>
25 
26 #include <opm/common/ErrorMacros.hpp>
27 
28 #include <opm/simulators/linalg/GraphColoring.hpp>
29 #include <opm/simulators/linalg/matrixblock.hh>
30 
31 namespace Opm
32 {
33 namespace detail
34 {
35 
37 template<class M>
38 void ghost_last_bilu0_decomposition (M& A, size_t interiorSize)
39 {
40  // iterator types
41  using rowiterator = typename M::RowIterator;
42  using coliterator = typename M::ColIterator;
43  using block = typename M::block_type;
44 
45  // implement left looking variant with stored inverse
46  for (rowiterator i = A.begin(); i.index() < interiorSize; ++i)
47  {
48  // coliterator is diagonal after the following loop
49  coliterator endij=(*i).end(); // end of row i
50  coliterator ij;
51 
52  // eliminate entries left of diagonal; store L factor
53  for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
54  {
55  // find A_jj which eliminates A_ij
56  coliterator jj = A[ij.index()].find(ij.index());
57 
58  // compute L_ij = A_jj^-1 * A_ij
59  (*ij).rightmultiply(*jj);
60 
61  // modify row
62  coliterator endjk=A[ij.index()].end(); // end of row j
63  coliterator jk=jj; ++jk;
64  coliterator ik=ij; ++ik;
65  while (ik!=endij && jk!=endjk)
66  if (ik.index()==jk.index())
67  {
68  block B(*jk);
69  B.leftmultiply(*ij);
70  *ik -= B;
71  ++ik; ++jk;
72  }
73  else
74  {
75  if (ik.index()<jk.index())
76  ++ik;
77  else
78  ++jk;
79  }
80  }
81 
82  // invert pivot and store it in A
83  if (ij.index()!=i.index())
84  DUNE_THROW(Dune::ISTLError,"diagonal entry missing");
85  try {
86  (*ij).invert(); // compute inverse of diagonal block
87  }
88  catch (Dune::FMatrixError & e) {
89  DUNE_THROW(Dune::ISTLError,"ILU failed to invert matrix block");
90  }
91  }
92 }
93 
95 template<class M, class CRS, class InvVector>
96 void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv)
97 {
98  // No need to do anything for 0 rows. Return to prevent indexing a
99  // a zero sized array.
100  if ( A.N() == 0 )
101  {
102  return;
103  }
104 
105  using size_type = typename M :: size_type;
106 
107  lower.clear();
108  upper.clear();
109  inv.clear();
110  lower.resize( A.N() );
111  upper.resize( A.N() );
112  inv.resize( A.N() );
113 
114  // Count the lower and upper matrix entries.
115  size_type numLower = 0;
116  size_type numUpper = 0;
117  const auto endi = A.end();
118  for (auto i = A.begin(); i != endi; ++i) {
119  const size_type iIndex = i.index();
120  size_type numLowerRow = 0;
121  for (auto j = (*i).begin(); j.index() < iIndex; ++j) {
122  ++numLowerRow;
123  }
124  numLower += numLowerRow;
125  numUpper += (*i).size() - numLowerRow - 1;
126  }
127  assert(numLower + numUpper + A.N() == A.nonzeroes());
128 
129  lower.reserveAdditional( numLower );
130 
131  // implement left looking variant with stored inverse
132  size_type row = 0;
133  size_type colcount = 0;
134  lower.rows_[ 0 ] = colcount;
135  for (auto i=A.begin(); i!=endi; ++i, ++row)
136  {
137  const size_type iIndex = i.index();
138 
139  // eliminate entries left of diagonal; store L factor
140  for (auto j=(*i).begin(); j.index() < iIndex; ++j )
141  {
142  lower.push_back( (*j), j.index() );
143  ++colcount;
144  }
145  lower.rows_[ iIndex+1 ] = colcount;
146  }
147 
148  assert(colcount == numLower);
149 
150  const auto rendi = A.beforeBegin();
151  row = 0;
152  colcount = 0;
153  upper.rows_[ 0 ] = colcount ;
154 
155  upper.reserveAdditional( numUpper );
156 
157  // NOTE: upper and inv store entries in reverse order, reverse here
158  // relative to ILU
159  for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
160  {
161  const size_type iIndex = i.index();
162 
163  // store in reverse row order
164  // eliminate entries left of diagonal; store L factor
165  for (auto j=(*i).beforeEnd(); j.index()>=iIndex; --j )
166  {
167  const size_type jIndex = j.index();
168  if( j.index() == iIndex )
169  {
170  inv[ row ] = (*j);
171  break;
172  }
173  else if ( j.index() >= i.index() )
174  {
175  upper.push_back( (*j), jIndex );
176  ++colcount ;
177  }
178  }
179  upper.rows_[ row+1 ] = colcount;
180  }
181  assert(colcount == numUpper);
182 }
183 
184 } // end namespace detail
185 
186 
187 template<class Matrix, class Domain, class Range, class ParallelInfoT>
188 Dune::SolverCategory::Category
189 ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::category() const
190 {
191  return std::is_same_v<ParallelInfoT, Dune::Amg::SequentialInformation> ?
192  Dune::SolverCategory::sequential : Dune::SolverCategory::overlapping;
193 }
194 
195 template<class Matrix, class Domain, class Range, class ParallelInfoT>
197 ParallelOverlappingILU0(const Matrix& A,
198  const int n, const field_type w,
199  MILU_VARIANT milu, bool redblack,
200  bool reorder_sphere)
201  : lower_(),
202  upper_(),
203  inv_(),
204  comm_(nullptr), w_(w),
205  relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
206  A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
207  milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
208 {
209  interiorSize_ = A.N();
210  // BlockMatrix is a Subclass of FieldMatrix that just adds
211  // methods. Therefore this cast should be safe.
212  update();
213 }
214 
215 template<class Matrix, class Domain, class Range, class ParallelInfoT>
217 ParallelOverlappingILU0(const Matrix& A,
218  const ParallelInfo& comm, const int n, const field_type w,
219  MILU_VARIANT milu, bool redblack,
220  bool reorder_sphere)
221  : lower_(),
222  upper_(),
223  inv_(),
224  comm_(&comm), w_(w),
225  relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
226  A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(n),
227  milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
228 {
229  interiorSize_ = A.N();
230  // BlockMatrix is a Subclass of FieldMatrix that just adds
231  // methods. Therefore this cast should be safe.
232  update();
233 }
234 
235 template<class Matrix, class Domain, class Range, class ParallelInfoT>
237 ParallelOverlappingILU0(const Matrix& A,
238  const field_type w, MILU_VARIANT milu, bool redblack,
239  bool reorder_sphere)
240  : ParallelOverlappingILU0( A, 0, w, milu, redblack, reorder_sphere )
241 {}
242 
243 template<class Matrix, class Domain, class Range, class ParallelInfoT>
245 ParallelOverlappingILU0(const Matrix& A,
246  const ParallelInfo& comm, const field_type w,
247  MILU_VARIANT milu, bool redblack,
248  bool reorder_sphere)
249  : lower_(),
250  upper_(),
251  inv_(),
252  comm_(&comm), w_(w),
253  relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
254  A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
255  milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
256 {
257  interiorSize_ = A.N();
258  // BlockMatrix is a Subclass of FieldMatrix that just adds
259  // methods. Therefore this cast should be safe.
260  update();
261 }
262 
263 template<class Matrix, class Domain, class Range, class ParallelInfoT>
265 ParallelOverlappingILU0(const Matrix& A,
266  const ParallelInfo& comm,
267  const field_type w, MILU_VARIANT milu,
268  size_type interiorSize, bool redblack,
269  bool reorder_sphere)
270  : lower_(),
271  upper_(),
272  inv_(),
273  comm_(&comm), w_(w),
274  relaxation_( std::abs( w - 1.0 ) > 1e-15 ),
275  interiorSize_(interiorSize),
276  A_(&reinterpret_cast<const Matrix&>(A)), iluIteration_(0),
277  milu_(milu), redBlack_(redblack), reorderSphere_(reorder_sphere)
278 {
279  // BlockMatrix is a Subclass of FieldMatrix that just adds
280  // methods. Therefore this cast should be safe.
281  update( );
282 }
283 
284 template<class Matrix, class Domain, class Range, class ParallelInfoT>
286 apply (Domain& v, const Range& d)
287 {
288  Range& md = reorderD(d);
289  Domain& mv = reorderV(v);
290 
291  // iterator types
292  using dblock = typename Range ::block_type;
293  using vblock = typename Domain::block_type;
294 
295  const size_type iEnd = lower_.rows();
296  const size_type lastRow = iEnd - 1;
297  size_type upperLoopStart = iEnd - interiorSize_;
298  size_type lowerLoopEnd = interiorSize_;
299  if (iEnd != upper_.rows())
300  {
301  OPM_THROW(std::logic_error,"ILU: number of lower and upper rows must be the same");
302  }
303 
304  // lower triangular solve
305  for (size_type i = 0; i < lowerLoopEnd; ++i)
306  {
307  dblock rhs( md[ i ] );
308  const size_type rowI = lower_.rows_[ i ];
309  const size_type rowINext = lower_.rows_[ i+1 ];
310 
311  for (size_type col = rowI; col < rowINext; ++col)
312  {
313  lower_.values_[ col ].mmv( mv[ lower_.cols_[ col ] ], rhs );
314  }
315 
316  mv[ i ] = rhs; // Lii = I
317  }
318 
319  for (size_type i = upperLoopStart; i < iEnd; ++i)
320  {
321  vblock& vBlock = mv[ lastRow - i ];
322  vblock rhs ( vBlock );
323  const size_type rowI = upper_.rows_[ i ];
324  const size_type rowINext = upper_.rows_[ i+1 ];
325 
326  for (size_type col = rowI; col < rowINext; ++col)
327  {
328  upper_.values_[ col ].mmv( mv[ upper_.cols_[ col ] ], rhs );
329  }
330 
331  // apply inverse and store result
332  inv_[ i ].mv( rhs, vBlock);
333  }
334 
335  copyOwnerToAll( mv );
336 
337  if( relaxation_ ) {
338  mv *= w_;
339  }
340  reorderBack(mv, v);
341 }
342 
343 template<class Matrix, class Domain, class Range, class ParallelInfoT>
344 template<class V>
346 copyOwnerToAll(V& v) const
347 {
348  if( comm_ ) {
349  comm_->copyOwnerToAll(v, v);
350  }
351 }
352 
353 template<class Matrix, class Domain, class Range, class ParallelInfoT>
354 void ParallelOverlappingILU0<Matrix,Domain,Range,ParallelInfoT>::
355 update()
356 {
357  // (For older DUNE versions the communicator might be
358  // invalid if redistribution in AMG happened on the coarset level.
359  // Therefore we check for nonzero size
360  if (comm_ && comm_->communicator().size() <= 0)
361  {
362  if (A_->N() > 0)
363  {
364  OPM_THROW(std::logic_error, "Expected a matrix with zero rows for an invalid communicator.");
365  }
366  else
367  {
368  // simply set the communicator to null
369  comm_ = nullptr;
370  }
371  }
372 
373  int ilu_setup_successful = 1;
374  std::string message;
375  const int rank = comm_ ? comm_->communicator().rank() : 0;
376 
377  std::unique_ptr<Matrix> ILU;
378 
379  if (redBlack_)
380  {
381  using Graph = Dune::Amg::MatrixGraph<const Matrix>;
382  Graph graph(*A_);
383  auto colorsTuple = colorVerticesWelshPowell(graph);
384  const auto& colors = std::get<0>(colorsTuple);
385  const auto& verticesPerColor = std::get<2>(colorsTuple);
386  auto noColors = std::get<1>(colorsTuple);
387  if ( reorderSphere_ )
388  {
389  ordering_ = reorderVerticesSpheres(colors, noColors, verticesPerColor,
390  graph, 0);
391  }
392  else
393  {
394  ordering_ = reorderVerticesPreserving(colors, noColors, verticesPerColor,
395  graph);
396  }
397  }
398 
399  std::vector<std::size_t> inverseOrdering(ordering_.size());
400  std::size_t index = 0;
401  for (const auto newIndex : ordering_)
402  {
403  inverseOrdering[newIndex] = index++;
404  }
405 
406  try
407  {
408  if (iluIteration_ == 0) {
409  // create ILU-0 decomposition
410  if (ordering_.empty())
411  {
412  ILU = std::make_unique<Matrix>(*A_);
413  }
414  else
415  {
416  ILU = std::make_unique<Matrix>(A_->N(), A_->M(),
417  A_->nonzeroes(), Matrix::row_wise);
418  auto& newA = *ILU;
419  // Create sparsity pattern
420  for (auto iter = newA.createbegin(), iend = newA.createend(); iter != iend; ++iter)
421  {
422  const auto& row = (*A_)[inverseOrdering[iter.index()]];
423  for (auto col = row.begin(), cend = row.end(); col != cend; ++col)
424  {
425  iter.insert(ordering_[col.index()]);
426  }
427  }
428  // Copy values
429  for (auto iter = A_->begin(), iend = A_->end(); iter != iend; ++iter)
430  {
431  auto& newRow = newA[ordering_[iter.index()]];
432  for (auto col = iter->begin(), cend = iter->end(); col != cend; ++col)
433  {
434  newRow[ordering_[col.index()]] = *col;
435  }
436  }
437  }
438 
439  switch (milu_)
440  {
442  detail::milu0_decomposition ( *ILU);
443  break;
445  detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
446  detail::signFunctor<typename Matrix::field_type> );
447  break;
449  detail::milu0_decomposition ( *ILU, detail::absFunctor<typename Matrix::field_type>,
450  detail::signFunctor<typename Matrix::field_type> );
451  break;
453  detail::milu0_decomposition ( *ILU, detail::identityFunctor<typename Matrix::field_type>,
454  detail::isPositiveFunctor<typename Matrix::field_type> );
455  break;
456  default:
457  if (interiorSize_ == A_->N())
458 #if DUNE_VERSION_LT(DUNE_GRID, 2, 8)
459  bilu0_decomposition( *ILU );
460 #else
461  Dune::ILU::blockILU0Decomposition( *ILU );
462 #endif
463  else
464  detail::ghost_last_bilu0_decomposition(*ILU, interiorSize_);
465  break;
466  }
467  }
468  else {
469  // create ILU-n decomposition
470  ILU = std::make_unique<Matrix>(A_->N(), A_->M(), Matrix::row_wise);
471  std::unique_ptr<detail::Reorderer> reorderer, inverseReorderer;
472  if (ordering_.empty())
473  {
474  reorderer.reset(new detail::NoReorderer());
475  inverseReorderer.reset(new detail::NoReorderer());
476  }
477  else
478  {
479  reorderer.reset(new detail::RealReorderer(ordering_));
480  inverseReorderer.reset(new detail::RealReorderer(inverseOrdering));
481  }
482 
483  milun_decomposition( *A_, iluIteration_, milu_, *ILU, *reorderer, *inverseReorderer );
484  }
485  }
486  catch (const Dune::MatrixBlockError& error)
487  {
488  message = error.what();
489  std::cerr << "Exception occurred on process " << rank << " during " <<
490  "setup of ILU0 preconditioner with message: "
491  << message<<std::endl;
492  ilu_setup_successful = 0;
493  }
494 
495  // Check whether there was a problem on some process
496  const bool parallel_failure = comm_ && comm_->communicator().min(ilu_setup_successful) == 0;
497  const bool local_failure = ilu_setup_successful == 0;
498  if (local_failure || parallel_failure)
499  {
500  throw Dune::MatrixBlockError();
501  }
502 
503  // store ILU in simple CRS format
504  detail::convertToCRS(*ILU, lower_, upper_, inv_);
505 }
506 
507 template<class Matrix, class Domain, class Range, class ParallelInfoT>
509 reorderD(const Range& d)
510 {
511  if (ordering_.empty())
512  {
513  // As d is non-const in the apply method of the
514  // solver casting away constness in this particular
515  // setting is not undefined. It is ugly though but due
516  // to the preconditioner interface of dune-istl.
517  return const_cast<Range&>(d);
518  }
519  else
520  {
521  reorderedD_.resize(d.size());
522  std::size_t i = 0;
523  for (const auto index : ordering_)
524  {
525  reorderedD_[index] = d[i++];
526  }
527  return reorderedD_;
528  }
529 }
530 
531 template<class Matrix, class Domain, class Range, class ParallelInfoT>
533 reorderV(Domain& v)
534 {
535  if (ordering_.empty())
536  {
537  return v;
538  }
539  else
540  {
541  reorderedV_.resize(v.size());
542  std::size_t i = 0;
543  for (const auto index : ordering_)
544  {
545  reorderedV_[index] = v[i++];
546  }
547  return reorderedV_;
548  }
549 }
550 
551 template<class Matrix, class Domain, class Range, class ParallelInfoT>
553 reorderBack(const Range& reorderedV, Range& v)
554 {
555  if (!ordering_.empty())
556  {
557  std::size_t i = 0;
558  for (const auto index : ordering_)
559  {
560  v[i++] = reorderedV[index];
561  }
562  }
563 }
564 
565 } // end namespace Opm
A two-step version of an overlapping Schwarz preconditioner using one step ILU0 as.
Definition: ParallelOverlappingILU0.hpp:145
ParallelOverlappingILU0(const Matrix &A, const int n, const field_type w, MILU_VARIANT milu, bool redblack=false, bool reorder_sphere=true)
Constructor.
Definition: ParallelOverlappingILU0_impl.hpp:197
Domain & reorderV(Domain &v)
Reorder V if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:533
Range & reorderD(const Range &d)
Reorder D if needed and return a reference to it.
Definition: ParallelOverlappingILU0_impl.hpp:509
void apply(Domain &v, const Range &d) override
Apply the preconditoner.
Definition: ParallelOverlappingILU0_impl.hpp:286
typename Domain::field_type field_type
The field type of the preconditioner.
Definition: ParallelOverlappingILU0.hpp:156
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
MILU_VARIANT
Definition: MILU.hpp:34
@ MILU_1
sum(dropped entries)
@ MILU_2
sum(dropped entries)
@ MILU_3
sum(|dropped entries|)
@ MILU_4
sum(dropped entries)
@ ILU
Do not perform modified ILU.
std::vector< std::size_t > reorderVerticesPreserving(const std::vector< int > &colors, int noColors, const std::vector< std::size_t > &verticesPerColor, const Graph &graph)
! Reorder colored graph preserving order of vertices with the same color.
Definition: GraphColoring.hpp:186
std::vector< std::size_t > reorderVerticesSpheres(const std::vector< int > &colors, int noColors, const std::vector< std::size_t > &verticesPerColor, const Graph &graph, typename Graph::VertexDescriptor root)
! Reorder Vetrices in spheres
Definition: GraphColoring.hpp:206
std::tuple< std::vector< int >, int, std::vector< std::size_t > > colorVerticesWelshPowell(const Graph &graph)
Color the vertices of graph.
Definition: GraphColoring.hpp:118