ldl.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_LDL_HH
#define DUNE_ISTL_LDL_HH
 
#if HAVE_SUITESPARSE_LDL || defined DOXYGEN
 
#include <iostream>
#include <memory>
#include <type_traits>
 
#ifdef __cplusplus
extern "C"
{
#include <ldl.h>
#include <amd.h>
}
#endif
 
#include <dune/common/exceptions.hh>
 
#include <dune/istl/bccsmatrixinitializer.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/solvertype.hh>
#include <dune/istl/solverfactory.hh>
 
namespace Dune {
  // forward declarations
  template<class M, class T, class TM, class TD, class TA>
  class SeqOverlappingSchwarz;
 
  template<class T, bool tag>
  struct SeqOverlappingSchwarzAssemblerHelper;
 
  template<class Matrix>
  class LDL
  {};
 
  template<typename T, typename A, int n, int m>
  class LDL<BCRSMatrix<FieldMatrix<T,n,m>,A > >
    : public InverseOperator<BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >,
                             BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > >
  {
    public:
    typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
    typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
    typedef Dune::ISTL::Impl::BCCSMatrix<T,int> LDLMatrix;
    typedef ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A>, int> MatrixInitializer;
    typedef Dune::BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > > domain_type;
    typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::Category::sequential;
    }
 
    LDL(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
    {
      //check whether T is a supported type
      static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
      setMatrix(matrix);
    }
 
    LDL(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false), verbose_(verbose)
    {
      //check whether T is a supported type
      static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
      setMatrix(matrix);
    }
 
    LDL(const Matrix& matrix, const ParameterTree& config)
      : LDL(matrix, config.get<int>("verbose", 0))
    {}
 
    LDL() : matrixIsLoaded_(false), verbose_(0)
    {}
 
    virtual ~LDL()
    {
      if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
        free();
    }
 
    virtual void apply(domain_type& x, range_type& b, InverseOperatorResult& res)
    {
      const int dimMat(ldlMatrix_.N());
      ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
      ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
      ldl_dsolve(dimMat, Y_, D_);
      ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
      ldl_permt(dimMat, reinterpret_cast<double*>(&x[0]), Y_, P_);
      // this is a direct solver
      res.iterations = 1;
      res.converged = true;
    }
 
    virtual void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
    {
      apply(x,b,res);
    }
 
    void apply(T* x, T* b)
    {
      const int dimMat(ldlMatrix_.N());
      ldl_perm(dimMat, Y_, b, P_);
      ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
      ldl_dsolve(dimMat, Y_, D_);
      ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
      ldl_permt(dimMat, x, Y_, P_);
    }
 
    void setOption([[maybe_unused]] unsigned int option, [[maybe_unused]] double value)
    {}
 
    void setMatrix(const Matrix& matrix)
    {
      if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
        free();
 
      if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
        ldlMatrix_.free();
      ldlMatrix_.setSize(MatrixDimension<Matrix>::rowdim(matrix),
                         MatrixDimension<Matrix>::coldim(matrix));
      ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
 
      copyToBCCSMatrix(initializer, matrix);
 
      decompose();
    }
 
    template<class S>
    void setSubMatrix(const Matrix& matrix, const S& rowIndexSet)
    {
      if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
        free();
 
      if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
        ldlMatrix_.free();
 
      ldlMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(matrix) / matrix.N(),
                         rowIndexSet.size()*MatrixDimension<Matrix>::coldim(matrix) / matrix.M());
      ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
 
      copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(matrix,rowIndexSet));
 
      decompose();
    }
 
    inline void setVerbosity(int v)
    {
      verbose_=v;
    }
 
    inline LDLMatrix& getInternalMatrix()
    {
      return ldlMatrix_;
    }
 
    void free()
    {
      delete [] D_;
      delete [] Y_;
      delete [] Lp_;
      delete [] Lx_;
      delete [] Li_;
      delete [] P_;
      delete [] Pinv_;
      ldlMatrix_.free();
      matrixIsLoaded_ = false;
    }
 
    inline const char* name()
    {
      return "LDL";
    }
 
    inline double* getD()
    {
      return D_;
    }
 
    inline int* getLp()
    {
      return Lp_;
    }
 
    inline int* getLi()
    {
      return Li_;
    }
 
    inline double* getLx()
    {
      return Lx_;
    }
 
    private:
    template<class M,class X, class TM, class TD, class T1>
    friend class SeqOverlappingSchwarz;
 
    friend struct SeqOverlappingSchwarzAssemblerHelper<LDL<Matrix>,true>;
 
    void decompose()
    {
      // allocate vectors
      const int dimMat(ldlMatrix_.N());
      D_ = new double [dimMat];
      Y_ = new double [dimMat];
      Lp_ = new int [dimMat + 1];
      Parent_ = new int [dimMat];
      Lnz_ = new int [dimMat];
      Flag_ = new int [dimMat];
      Pattern_ = new int [dimMat];
      P_ = new int [dimMat];
      Pinv_ = new int [dimMat];
 
      double Info [AMD_INFO];
      if(amd_order (dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), P_, (double *) NULL, Info) < AMD_OK)
        DUNE_THROW(InvalidStateException,"Error: AMD failed!");
      if(verbose_ > 0)
        amd_info (Info);
      // compute the symbolic factorisation
      ldl_symbolic(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
      // initialise those entries of additionalVectors_ whose dimension is known only now
      Lx_ = new double [Lp_[dimMat]];
      Li_ = new int [Lp_[dimMat]];
      // compute the numeric factorisation
      const int rank(ldl_numeric(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), ldlMatrix_.getValues(),
                                 Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
      // free temporary vectors
      delete [] Flag_;
      delete [] Pattern_;
      delete [] Parent_;
      delete [] Lnz_;
 
      if(rank!=dimMat)
        DUNE_THROW(InvalidStateException,"Error: LDL factorisation failed!");
    }
 
    LDLMatrix ldlMatrix_;
    bool matrixIsLoaded_;
    int verbose_;
    int* Lp_;
    int* Parent_;
    int* Lnz_;
    int* Flag_;
    int* Pattern_;
    int* P_;
    int* Pinv_;
    double* D_;
    double* Y_;
    double* Lx_;
    int* Li_;
  };
 
  template<typename T, typename A, int n, int m>
  struct IsDirectSolver<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
  {
    enum {value = true};
  };
 
  template<typename T, typename A, int n, int m>
  struct StoresColumnCompressed<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
  {
    enum {value = true};
  };
 
  struct LDLCreator {
    template<class F> struct isValidBlock : std::false_type{};
    template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
 
    template<typename TL, typename M>
    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
                                          typename Dune::TypeListElement<2, TL>::type>>
    operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
      std::enable_if_t<
                isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
    {
      int verbose = config.get("verbose", 0);
      return std::make_shared<Dune::LDL<M>>(mat,verbose);
    }
 
    // second version with SFINAE to validate the template parameters of LDL
    template<typename TL, typename M>
    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
                                          typename Dune::TypeListElement<2, TL>::type>>
    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
      std::enable_if_t<
                !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
    {
      DUNE_THROW(UnsupportedType,
        "Unsupported Type in LDL (only double and std::complex<double> supported)");
    }
  };
  DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator());
 
} // end namespace Dune
 
 
#endif //HAVE_SUITESPARSE_LDL
#endif //DUNE_ISTL_LDL_HH