domainDecompositionPreconditioner.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef LINALG_PATCHDOMAINDECOMPOSITIONPPRECONDITIONER_HH
#define LINALG_PATCHDOMAINDECOMPOSITIONPPRECONDITIONER_HH
 
#include <vector>
 
#include <boost/timer/timer.hpp>
 
#include <dune/common/dynvector.hh>
#include "dune/istl/preconditioner.hh"
 
#include "linalg/dynamicMatrix.hh"
#include "linalg/localMatrices.hh"
#include "linalg/symmetricOperators.hh"
#include "linalg/threadedMatrix.hh"
#include "fem/supports.hh"
#include "utilities/scalar.hh"
#include "utilities/threading.hh"
#include "utilities/timing.hh"
 
 
// #define INEXACTSTORAGE
// #include <random> // for mixed precision *testing*
// extern double epsilon;
 
 
 
namespace Kaskade
{
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  class InverseLocalMatrixStorageBase
  {
  public:
    InverseLocalMatrixStorageBase(std::vector<size_t> const& ap)
    : dofs(ap)
    {
    }
 
    std::vector<size_t> const& indices() const
    {
      return dofs;
    }
 
    size_t storage() const
    {
      return storageSize;
    }
 
    size_t flops() const
    {
      return computeFlops;
    }
 
  protected:
    size_t storageSize;
    size_t computeFlops;
 
    template <int m, class LocalMatrix, class Real>
    void regularize(LocalMatrix& aloc, Real r) const
    {
      if (r > 0)
      {
        Real Amax = 0;
        for (int i=0; i<aloc.N(); ++i)
          Amax = std::max(static_cast<Real>(std::sqrt(frobenius_norm2(aloc[i][i]))),Amax);
        for (int i=0; i<aloc.N(); ++i)
          aloc[i][i] += Amax/r * unitMatrix<Real,m>();
      }
    }
 
  private:
    std::vector<size_t> dofs;
  };
 
  // ----------------------------------------------------------------------------------------------
 
  template <class Scalar>
  class DenseInverseStorageTag {};
 
  template <int m, class Scalar=float>
  class DenseInverseLocalMatrixStorage: public InverseLocalMatrixStorageBase
  {
    using StorageEntry = Dune::FieldMatrix<Scalar,m,m>;
 
  public:
 
    using field_type = Scalar;
 
    template <class T, class Index>
    DenseInverseLocalMatrixStorage(std::vector<size_t> const& ap,
                                   NumaBCRSMatrix<Dune::FieldMatrix<T,m,m>,Index> const& A,
                                   field_type regularizeForCondition, DenseInverseStorageTag<Scalar>)
    : InverseLocalMatrixStorageBase(ap)
    {
      alocinv = full(A,indices(),indices());
 
#ifdef INEXACTSTORAGE
      using Real = typename ScalarTraits<Scalar>::Real;
      Real anorm = two_norm(alocinv);
#endif
 
      regularize<m>(alocinv,regularizeForCondition);
 
      invertSpd(alocinv); // invert the local block
 
#ifdef INEXACTSTORAGE
      DynamicMatrix<Dune::FieldMatrix<Scalar,m,m>> delta(alocinv.N(),alocinv.M());
 
      std::random_device rd;
      std::mt19937 gen(rd());
      std::uniform_real_distribution<> dis(-1.0,1.0);
      for (int i=0; i<alocinv.N(); ++i)
        for (int j=0; j<alocinv.M(); ++j)
          for (int k=0; k<m; ++k)
            for (int l=0; l<m; ++l)
              delta[i][j][k][l] = dis(gen);
      Real dnorm = two_norm(delta);
 
      Real epsilon = ::epsilon;
      delta *= epsilon / anorm / dnorm;
      alocinv += delta;
#endif
 
      storageSize = sizeof(StorageEntry)*alocinv.N()*alocinv.M() + sizeof(size_t)*indices().size();
      auto n = alocinv.N() * m;
      computeFlops = n*(2*n-1); // in every row n multiplications and n-1 additions
    }
 
 
    template <class Vector>
    void mv(Vector const& x, Vector& y) const
    {
      alocinv.mv(x,y);
    }
 
 
  private:
    DynamicMatrix<StorageEntry> alocinv;
  };
 
  // ----------------------------------------------------------------------------------------------
 
  std::array<std::vector<size_t>,3> nestedDissection(std::vector<std::pair<size_t,size_t>> const& edges,
                                                     std::vector<size_t> const& edgeStart);
 
  template <class StorageScalar>
  class NestedDissection
  {
  public:
    NestedDissection() = default;
 
    template <class SparseMatrix>
    NestedDissection(SparseMatrix const& A, size_t skipEntries, size_t skipDimension);
 
    template <class FactorizationScalar>
    NestedDissection<StorageScalar>& operator=(NestedDissection<FactorizationScalar>&& nd);
 
    template <class T>
    void solve(T& y) const;
 
    size_t storageSize() const;
    size_t flops() const;
 
    size_t size()
    {
      return sz;
    }
 
  private:
    template <class FactorizationScalar>
    friend class NestedDissection;
 
    std::vector<size_t> v1, v2, s;
    std::unique_ptr<NestedDissection<StorageScalar>> l1FactorND, l2FactorND;
    DynamicMatrix<StorageScalar> r1, r2;
    std::vector<StorageScalar> sFactor;
    size_t sz;
 
    template <int m>
    void trisolve(DynamicMatrix<StorageScalar>& y, char trans) const;
  };
 
  template <class StorageScalar, class FactorizationScalar>
  class NestedDissectionStorageTag
  {
  public:
 
    NestedDissectionStorageTag(size_t skipEntries_=100, size_t skipDimension_=32)
    : skipEntries(skipEntries_), skipDimension(skipDimension_)
    {}
 
    size_t skipEntries;
 
    size_t skipDimension;
  };
 
  template <int m, class StorageScalar, class FactorizationScalar>
  class NestedDissectionCholeskyLocalMatrixStorage: public InverseLocalMatrixStorageBase
  {
  public:
    using field_type = StorageScalar;
 
    template <class T, class Index>
    NestedDissectionCholeskyLocalMatrixStorage(std::vector<size_t> const& ap,
                                               NumaBCRSMatrix<Dune::FieldMatrix<T,m,m>,Index> const& A,
                                               FactorizationScalar regularizeForCondition,
                                               NestedDissectionStorageTag<StorageScalar,FactorizationScalar> storageTag)
    : InverseLocalMatrixStorageBase(ap)
    {
      using LocalMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<T,m,m>>;
 
      LocalMatrix aloc = submatrix<LocalMatrix>(A,indices(),indices());
 
      regularize<m>(aloc,regularizeForCondition);
 
      choleskyFactor = NestedDissection<FactorizationScalar>(aloc,storageTag.skipEntries,
                                                             storageTag.skipDimension);  // factorize the local block and compress
 
      storageSize = choleskyFactor.storageSize();
      computeFlops = choleskyFactor.flops();
    }
 
 
    template <class Vector>
    void mv(Vector const& x, Vector& y) const
    {
      y = x;
      choleskyFactor.solve(y);
    }
 
  private:
    NestedDissection<StorageScalar> choleskyFactor;
  };
 
 
  // ----------------------------------------------------------------------------------------------
 
  namespace DomainDecompositionPreconditionerDetail
  {
    template <class Scalar>
    std::vector<Scalar> choleskyFactor(DynamicMatrix<Scalar> const& A);
 
    template <class Scalar>
    void choleskySolve(int n, int nrhs, std::vector<Scalar> const& L, Scalar* b);
  }
 
  template <class Scalar>
  class DenseCholeskyStorageTag {};
 
  template <int m, class Scalar=float>
  class DenseCholeskyLocalMatrixStorage: public InverseLocalMatrixStorageBase
  {
    using StorageEntry = Dune::FieldMatrix<Scalar,m,m>;
 
  public:
 
    using field_type = Scalar;
 
    template <class T, class Index>
    DenseCholeskyLocalMatrixStorage(std::vector<size_t> const& ap,
                                    NumaBCRSMatrix<Dune::FieldMatrix<T,m,m>,Index> const& A,
                                    field_type regularizeForCondition, DenseCholeskyStorageTag<Scalar>)
    : InverseLocalMatrixStorageBase(ap)
    {
      DynamicMatrix<Scalar> aloc = DynamicMatrixDetail::flatMatrix(full(A,indices(),indices()));
      n = aloc.N();
 
      regularize<1>(aloc,regularizeForCondition);
 
      choleskyFactor = DomainDecompositionPreconditionerDetail::choleskyFactor(aloc); // factorize the local block
      storageSize = sizeof(Scalar)*choleskyFactor.size() + sizeof(size_t)*indices().size();
      computeFlops =  2*n*n;
 
#ifdef INEXACTSTORAGE
      std::random_device rd;
      std::mt19937 gen(rd());
      std::uniform_real_distribution<> dis(-1.0,0.0);
 
      DynamicMatrix<Scalar> L(n,n), delta(n,n);
      L = 0;
      delta = 0;
      auto lpos = begin(choleskyFactor);
      for (int j=0; j<n; ++j)
        for (int i=j; i<n; ++i)
        {
          L[i][j] = *lpos;
          ++lpos;
          delta[i][j] = dis(gen) / (i-j+1);
        }
 
      invert(L);
 
      Scalar epsilon = ::epsilon;
      delta *= epsilon/two_norm(delta)/two_norm(L);
 
      lpos = begin(choleskyFactor);
      for (int j=0; j<n; ++j)
        for (int i=j; i<n; ++i)
        {
          *lpos += delta[i][j];
          ++lpos;
        }
 
#endif
    }
 
 
    template <class Vector>
    void mv(Vector const& x, Vector& y) const
    {
      y = x;
      DomainDecompositionPreconditionerDetail::choleskySolve(n,1,choleskyFactor,&y[0][0]);
    }
 
 
  private:
    std::vector<Scalar> choleskyFactor; // L in LAPACK packed storage
    int n;                              // size of L (with scalar entries, not blocks of size m)
  };
 
 
 
  // ----------------------------------------------------------------------------------------------
  // ----------------------------------------------------------------------------------------------
 
  template <class StorageTag, int components>
  struct PatchDomainDecompositionPreconditionerTraits
  {
    using type = void;
  };
 
  template <class Scalar, int components>
  struct PatchDomainDecompositionPreconditionerTraits<DenseInverseStorageTag<Scalar>,components>
  {
    using type = DenseInverseLocalMatrixStorage<components,Scalar>;
  };
 
  template <class Scalar, int components>
  struct PatchDomainDecompositionPreconditionerTraits<DenseCholeskyStorageTag<Scalar>,components>
  {
    using type = DenseCholeskyLocalMatrixStorage<components,Scalar>;
  };
 
  template <class StorageScalar, class FactorizationScalar, int components>
  struct PatchDomainDecompositionPreconditionerTraits<NestedDissectionStorageTag<StorageScalar,FactorizationScalar>,components>
  {
    using type = NestedDissectionCholeskyLocalMatrixStorage<components,StorageScalar,FactorizationScalar>;
  };
 
  template <class Space, int m, class StorageTag,
            class SparseMatrixIndex = std::size_t>
  class PatchDomainDecompositionPreconditioner: public SymmetricPreconditioner<typename Space::template Element<m>::type::StorageType,
                                                                               typename Space::template Element<m>::type::StorageType>
  {
    using InverseLocalMatrix = typename PatchDomainDecompositionPreconditionerTraits<StorageTag,m>::type;
 
  public:
    using domain_type = typename Space::template Element<m>::type::StorageType;
    using range_type = domain_type;
 
    template <class Matrix>
    PatchDomainDecompositionPreconditioner(Space const& space, Matrix const& A, 
                                           typename Matrix::field_type regularizeForCondition = -1.0,
                                           StorageTag storageTag = StorageTag(),
                                           int nTasks = 0)
    {
      assert(A.N()>0);
      assert(A.N()==A.M());
      Timings& timer = Timings::instance();
      timer.start("make Schwarz preconditioner");
 
      // for each patch around a vertex, store the global dof indices
      std::vector<std::vector<size_t>> patchDofs(space.gridView().size(space.gridView().dimension));
      
      auto patches = computePatches(space.gridView());
      assert(patches.size()==patchDofs.size());
      
      int nThreads = nTasks<=0? 4*NumaThreadPool::instance().cpus()
                              : std::min(nTasks,(int)patches.size());
      inverseLocalMatrices.resize(nThreads);
 
      TaskTiming taskTiming(nThreads);
      
      // Compute the inverse of each block's local stiffness matrix in parallel.
      parallelFor([&patches,&A,&space,regularizeForCondition,storageTag,&taskTiming,this](size_t j, size_t s)
      {
        taskTiming.start(j);
        // For all patches (in our parallel loop range)
        for (size_t i=uniformWeightRangeStart(j,s,patches.size()); i<uniformWeightRangeStart(j+1,s,patches.size()); ++i)
        {
          auto const& ap = algebraicPatch(space,patches[i]);
          // todo: (catch possible exception, store by std::exception_ptr (which needs to be passed from main thread and locked), then rethrow in main thread
          try 
          {
            InverseLocalMatrix aloc = InverseLocalMatrix(ap,A,regularizeForCondition,storageTag);
            this->inverseLocalMatrices[j].push_back(std::move(aloc));
          } 
          catch (DetailedException const& ex)
          {
            std::cout << "exception in inverting block " << i << " in thread " << j << " of " << s << ":\n";
            std::cout << ex.what() << "\n";
          }
        }
        taskTiming.stop(j);
      },nThreads);
 
      timer.stop("make Schwarz preconditioner");
 
std::cerr << "Schwarz preconditioner storage size: " << storage()/(1024*1024) << " MB\n"
          << "                       flop count:   " << flops()/(1024*1024) << " M\n";
 
      std::ofstream f("out.gnu");
      f << taskTiming;
    }
 
    virtual typename Space::Scalar applyDp(domain_type& x, range_type const& y)
    {
      apply(x,y);
      return x*y;
    }
 
    virtual void apply(domain_type& x, range_type const& y)
    {      
      using StorageScalar = typename InverseLocalMatrix::field_type;
      using LocalVector = Dune::BlockVector<Dune::FieldVector<StorageScalar,m>>;
 
      // TODO: consider scattering into local copies of global vector? Maybe dependent on the actual
      // local matrix size? Or range locking?
      // But for 3D order 3&4 locking appears not to be a performance bottleneck on
      // a quadcore (as of 2018-11).
 
      ScopedTimingSection timer("Schwarz preconditioner apply");
      
      std::mutex rhsAccess;
      parallelFor([&rhsAccess,&x,&y,this](int j, int s)
      {
        for (auto const& alocinv: this->inverseLocalMatrices[j])
        {
          auto const& indices = alocinv.indices();
          int const n = indices.size();
 
          LocalVector xloc(n);
          LocalVector yloc(n);
          
          for (int i=0; i<n; ++i)               // gather relevant vector entries from x
            yloc[i] = y[indices[i]];
          
          alocinv.mv(yloc,xloc);                // perform local multiplication
          
          std::lock_guard<std::mutex> lock(rhsAccess);
          for (int i=0; i<n; ++i)               // scatter local result to global vector y
            x[indices[i]] += xloc[i];
        }
      },inverseLocalMatrices.size());
    }
 
    virtual bool requiresInitializedInput() const
    {
      return true;
    }
 
    size_t storage() const
    {
      size_t s = 0;
      for (auto const& as: inverseLocalMatrices)
        for (auto const& a: as)
          s += a.storage();
      return s;
    }
 
    size_t flops() const
    {
      size_t f = 0;
      for (auto const& as: inverseLocalMatrices)
      {
        for (auto const& a: as)
          f += a.flops();
      }
      return f;
    }
 
  private:
    std::vector<std::vector<InverseLocalMatrix>> inverseLocalMatrices;
  };
}
 
#endif