Loading...
Searching...
No Matches
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
solvers.hh
Implementations of the inverse operator interface.
solvertype.hh
Templates characterizing the type of a solver.
bccsmatrixinitializer.hh
solverregistry.hh
DUNE_REGISTER_SOLVER
#define DUNE_REGISTER_SOLVER(name,...)
Definition solverregistry.hh:16
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::SeqOverlappingSchwarz
friend class SeqOverlappingSchwarz
Definition ldl.hh:303
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::domain_type
Dune::BlockVector< FieldVector< T, m >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, m > > > domain_type
The type of the domain of the solver.
Definition ldl.hh:85
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::Matrix
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > Matrix
The matrix type.
Definition ldl.hh:78
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::getLp
int * getLp()
Get factorization Lp.
Definition ldl.hh:278
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::getInternalMatrix
LDLMatrix & getInternalMatrix()
Return the column compress matrix.
Definition ldl.hh:237
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::apply
void apply(domain_type &x, range_type &b, InverseOperatorResult &res) override
Apply inverse operator,.
Definition ldl.hh:152
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::LDLMatrix
Dune::ISTL::Impl::BCCSMatrix< T, int > LDLMatrix
The corresponding SuperLU Matrix type.
Definition ldl.hh:81
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::getLx
double * getLx()
Get factorization Lx.
Definition ldl.hh:296
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::setVerbosity
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition ldl.hh:228
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::apply
void apply(T *x, T *b)
Additional apply method with c-arrays in analogy to superlu.
Definition ldl.hh:176
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::range_type
Dune::BlockVector< FieldVector< T, n >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, n > > > range_type
The type of the range of the solver.
Definition ldl.hh:87
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::apply
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res) override
apply inverse operator, with given convergence criteria.
Definition ldl.hh:166
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::~LDL
virtual ~LDL()
Default constructor.
Definition ldl.hh:145
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::free
void free()
Free allocated space.
Definition ldl.hh:246
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::getLi
int * getLi()
Get factorization Li.
Definition ldl.hh:287
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::matrix_type
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > matrix_type
Definition ldl.hh:79
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::setMatrix
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition ldl.hh:190
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::setOption
void setOption(unsigned int option, double value)
Definition ldl.hh:186
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::LDL
LDL(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition ldl.hh:120
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::getD
double * getD()
Get factorization diagonal matrix D.
Definition ldl.hh:269
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::category
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category).
Definition ldl.hh:90
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::LDL
LDL(const Matrix &matrix, const ParameterTree &config)
Constructs the LDL solver.
Definition ldl.hh:136
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::MatrixInitializer
ISTL::Impl::BCCSMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A >, int > MatrixInitializer
Type of an associated initializer class.
Definition ldl.hh:83
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::name
const char * name()
Get method name.
Definition ldl.hh:260
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::setSubMatrix
void setSubMatrix(const Matrix &matrix, const S &rowIndexSet)
Definition ldl.hh:207
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::LDL
LDL()
Default constructor.
Definition ldl.hh:141
Dune::LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > >::LDL
LDL(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition ldl.hh:104
Dune::IsDirectSolver< LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > > >::value
@ value
Definition ldl.hh:364
Dune::StoresColumnCompressed< LDL< BCRSMatrix< FieldMatrix< T, n, m >, A > > >::value
@ value
Definition ldl.hh:370
Dune
Definition allocator.hh:11
Dune::get
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
Dune::MatrixDimension::coldim
static auto coldim(const M &)
Definition matrixutils.hh:219
Dune::MatrixDimension::rowdim
static auto rowdim(const M &)
Definition matrixutils.hh:214
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:467
Dune::BCRSMatrix::M
size_type M() const
number of columns (counted in blocks)
Definition bcrsmatrix.hh:2010
Dune::BCRSMatrix::N
size_type N() const
number of rows (counted in blocks)
Definition bcrsmatrix.hh:2004
Dune::BlockVector
A vector of blocks with memory management.
Definition bvector.hh:392
Dune::SeqOverlappingSchwarz
Sequential overlapping Schwarz preconditioner.
Definition overlappingschwarz.hh:755
Dune::SeqOverlappingSchwarzAssemblerHelper
Definition overlappingschwarz.hh:694
Dune::LDL
Use the LDL package to directly solve linear systems – empty default class.
Definition ldl.hh:56
Dune::FieldMatrix
Definition matrixutils.hh:27
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition solver.hh:50
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition solver.hh:69
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition solver.hh:75
Dune::InverseOperator
Abstract base class for all solvers.
Definition solver.hh:101
Dune::SolverCategory::Category
Category
Definition solvercategory.hh:23
Dune::SolverCategory::sequential
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25
Dune::UnsupportedType
Definition solverregistry.hh:97
Dune::IsDirectSolver
Definition solvertype.hh:16
Dune::StoresColumnCompressed
Definition solvertype.hh:30