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/solverregistry.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;

    SolverCategory::Category category() const override
    {
      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();
    }

    void apply(domain_type& x, range_type& b, InverseOperatorResult& res) override
    {
      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;
    }

    void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
    {
      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};
  };

  DUNE_REGISTER_SOLVER("ldl",
                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
                       {
                         using OpTraits = decltype(opTraits);
                         using M = typename OpTraits::matrix_type;
                         // works only for sequential operators
                         if constexpr (OpTraits::isParallel){
                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
                             DUNE_THROW(Dune::InvalidStateException, "LDL works only for sequential operators.");
                         }
                         // check if LDL<M>* is convertible to
                         // InverseOperator*. This allows only the explicit
                         // specialized variants of LDL
                         if constexpr (std::is_convertible_v<LDL<M>*,
                                       Dune::InverseOperator<typename OpTraits::domain_type,
                                       typename OpTraits::range_type>*> &&
                                       std::is_same_v<typename FieldTraits<M>::field_type, double>
                                       ){
                           const auto& A = opTraits.getAssembledOpOrThrow(op);
                           const M& mat = A->getmat();
                           int verbose = config.get("verbose", 0);
                           return std::make_shared<LDL<M>>(mat,verbose);
                         }
                         DUNE_THROW(UnsupportedType,
                                    "Unsupported Type in LDL (only FieldMatrix<double,...> supported)");
                       });

} // end namespace Dune


#endif //HAVE_SUITESPARSE_LDL
#endif //DUNE_ISTL_LDL_HH