ilu.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_ILU_HH
#define DUNE_ISTL_ILU_HH
 
#include <cmath>
#include <complex>
#include <map>
#include <vector>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
 
#include "istlexception.hh"
 
namespace Dune {
 
  namespace ILU {
 
    template<class M>
    void blockILU0Decomposition (M& A)
    {
      // iterator types
      typedef typename M::RowIterator rowiterator;
      typedef typename M::ColIterator coliterator;
      typedef typename M::block_type block;
 
      // implement left looking variant with stored inverse
      rowiterator endi=A.end();
      for (rowiterator i=A.begin(); i!=endi; ++i)
      {
        // coliterator is diagonal after the following loop
        coliterator endij=(*i).end();           // end of row i
        coliterator ij;
 
        // eliminate entries left of diagonal; store L factor
        for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
        {
          // find A_jj which eliminates A_ij
          coliterator jj = A[ij.index()].find(ij.index());
 
          // compute L_ij = A_jj^-1 * A_ij
          Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
 
          // modify row
          coliterator endjk=A[ij.index()].end();                 // end of row j
          coliterator jk=jj; ++jk;
          coliterator ik=ij; ++ik;
          while (ik!=endij && jk!=endjk)
            if (ik.index()==jk.index())
            {
              block B(*jk);
              Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
              *ik -= B;
              ++ik; ++jk;
            }
            else
            {
              if (ik.index()<jk.index())
                ++ik;
              else
                ++jk;
            }
        }
 
        // invert pivot and store it in A
        if (ij.index()!=i.index())
          DUNE_THROW(ISTLError,"diagonal entry missing");
        try {
          Impl::asMatrix(*ij).invert();   // compute inverse of diagonal block
        }
        catch (Dune::FMatrixError & e) {
          DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
                     << i.index() << "][" << ij.index() << "]" << e.what();
                     th__ex.r=i.index(); th__ex.c=ij.index(););
        }
      }
    }
 
    template<class M, class X, class Y>
    void blockILUBacksolve (const M& A, X& v, const Y& d)
    {
      // iterator types
      typedef typename M::ConstRowIterator rowiterator;
      typedef typename M::ConstColIterator coliterator;
      typedef typename Y::block_type dblock;
      typedef typename X::block_type vblock;
 
      // lower triangular solve
      rowiterator endi=A.end();
      for (rowiterator i=A.begin(); i!=endi; ++i)
      {
        // We need to be careful here: Directly using
        // auto rhs = Impl::asVector(d[ i.index() ]);
        // is not OK in case this is a proxy. Hence
        // we first have to copy the value. Notice that
        // this is still not OK, if the vector type itself returns
        // proxy references.
        dblock rhsValue(d[i.index()]);
        auto&& rhs = Impl::asVector(rhsValue);
        for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
          Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
        Impl::asVector(v[i.index()]) = rhs;           // Lii = I
      }
 
      // upper triangular solve
      rowiterator rendi=A.beforeBegin();
      for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
      {
        // We need to be careful here: Directly using
        // auto rhs = Impl::asVector(v[ i.index() ]);
        // is not OK in case this is a proxy. Hence
        // we first have to copy the value. Notice that
        // this is still not OK, if the vector type itself returns
        // proxy references.
        vblock rhsValue(v[i.index()]);
        auto&& rhs = Impl::asVector(rhsValue);
        coliterator j;
        for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
          Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
        auto&& vi = Impl::asVector(v[i.index()]);
        Impl::asMatrix(*j).mv(rhs,vi);           // diagonal stores inverse!
      }
    }
 
    // recursive function template to access first entry of a matrix
    template<class M>
    typename M::field_type& firstMatrixElement (M& A,
                                                [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
    {
      return firstMatrixElement(*(A.begin()->begin()));
    }
 
    template<class K>
    K& firstMatrixElement (K& A,
                           [[maybe_unused]] typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
    {
      return<