1#ifndef DUNE_FEM_LDLSOLVER_HH
2#define DUNE_FEM_LDLSOLVER_HH
10#include <dune/common/exceptions.hh>
11#include <dune/common/hybridutilities.hh>
20#if HAVE_SUITESPARSE_LDL
45template<
class DiscreteFunction,
46 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
47class LDLInverseOperator;
49template<
class DiscreteFunction,
class Matrix >
50struct LDLInverseOperatorTraits
52 typedef DiscreteFunction DiscreteFunctionType;
53 typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
56 typedef OperatorType PreconditionerType;
58 typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
60 typedef LDLInverseOperator< DiscreteFunction, Matrix > InverseOperatorType;
61 typedef SolverParameter SolverParameterType;
73template<
class DF,
class Matrix >
74class LDLInverseOperator :
public InverseOperatorInterface< LDLInverseOperatorTraits< DF, Matrix > >
76 typedef LDLInverseOperatorTraits< DF, Matrix > Traits;
77 typedef InverseOperatorInterface< Traits > BaseType;
79 friend class InverseOperatorInterface< Traits >;
81 typedef LDLInverseOperator< DF, Matrix > ThisType;
83 typedef typename BaseType :: SolverDiscreteFunctionType
84 SolverDiscreteFunctionType;
86 typedef typename BaseType :: OperatorType OperatorType;
87 typedef typename BaseType :: AssembledOperatorType AssembledOperatorType;
90 typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,int> CCSMatrixType;
92 typedef typename SolverDiscreteFunctionType::DofType DofType;
93 typedef typename SolverDiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
101 LDLInverseOperator(
const double& redEps,
const double& absLimit,
const int& maxIter,
const bool& verbose,
102 const ParameterReader ¶meter = Parameter::container() ) :
103 verbose_(verbose &&
Dune::Fem::Parameter::verbose()), ccsmat_()
111 LDLInverseOperator(
const double& redEps,
const double& absLimit,
const int& maxIter,
112 const ParameterReader ¶meter = Parameter::container() ) :
113 verbose_(parameter.getValue<bool>(
"fem.solver.verbose",false)), ccsmat_()
116 LDLInverseOperator(
const SolverParameter ¶meter = SolverParameter(Parameter::container()) )
117 : BaseType(parameter), verbose_(BaseType::verbose())
127 LDLInverseOperator(
const OperatorType& op,
const double& redEps,
const double& absLimit,
const int& maxIter,
const bool& verbose,
128 const ParameterReader ¶meter = Parameter::container() ) :
129 verbose_(verbose), ccsmat_(), isloaded_(false)
140 LDLInverseOperator(
const OperatorType& op,
const double& redEps,
const double& absLimit,
const int& maxIter,
141 const ParameterReader ¶meter = Parameter::container() ) :
142 verbose_(parameter.getValue<bool>(
"fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
147 LDLInverseOperator(
const OperatorType& op,
const ParameterReader ¶meter = Parameter::container() ) :
148 verbose_(parameter.getValue<bool>(
"fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
154 ~LDLInverseOperator()
159 void bind(
const OperatorType& op )
163 BaseType::bind( op );
173 template<
typename... A>
174 void prepare(A... )
const
176 if( ! assembledOperator_ )
177 DUNE_THROW(NotImplemented,
"LDLInverseOperator only works for assembled systems!");
181 ccsmat_ = assembledOperator_->exportMatrix();
188 virtual void finalize()
210 template<
typename... DFs>
211 void apply(
const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest)
const
214 std::vector<DofType> vecArg(arg.size());
215 auto vecArgIt(vecArg.begin());
216 Hybrid::forEach(std::make_index_sequence<
sizeof...(DFs)>{},
217 [&](
auto i){vecArgIt=std::copy(std::get<i>(arg).dbegin(),std::get<i>(arg).dend(),vecArgIt);});
218 std::vector<DofType> vecDest(dest.size());
220 apply(vecArg.data(),vecDest.data());
222 auto vecDestIt(vecDest.begin());
223 Hybrid::forEach(std::make_index_sequence<
sizeof...(DFs)>{},[&](
auto i){
for(
auto& dof:
dofs(std::get<i>(dest))) dof=(*(vecDestIt++));});
227 void printTexInfo(std::ostream& out)
const
229 out<<
"Solver: LDL direct solver"<<std::endl;
233 void printDecompositionInfo()
const
273 CCSMatrixType& getCCSMatrix()
285 void apply(
const DofType* arg, DofType* dest)
const
290 const std::size_t dimMat(ccsmat_.N());
291 ldl_perm(dimMat, Y_,
const_cast<DofType*
>(arg), P_);
292 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
293 ldl_dsolve(dimMat, Y_, D_);
294 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
295 ldl_permt(dimMat, dest, Y_, P_);
297 const_cast<ThisType*
>(
this)->finalize();
306 int apply(
const SolverDiscreteFunctionType& arg,
307 SolverDiscreteFunctionType& dest)
const
309 apply(arg.leakPointer(),dest.leakPointer());
315 using BaseType :: assembledOperator_;
318 mutable CCSMatrixType ccsmat_;
319 mutable bool isloaded_ =
false;
321 mutable int* Parent_;
324 mutable int* Pattern_;
329 mutable DofType* Lx_;
331 mutable double info_[AMD_INFO];
334 void decompose()
const
337 const std::size_t dimMat(ccsmat_.N());
338 D_ =
new DofType [dimMat];
339 Y_ =
new DofType [dimMat];
340 Lp_ =
new int [dimMat + 1];
341 Parent_ =
new int [dimMat];
342 Lnz_ =
new int [dimMat];
343 Flag_ =
new int [dimMat];
344 Pattern_ =
new int [dimMat];
345 P_ =
new int [dimMat];
346 Pinv_ =
new int [dimMat];
348 if(amd_order (dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), P_, (DofType *) NULL, info_) < AMD_OK)
349 DUNE_THROW(InvalidStateException,
"LDL Error: AMD failed!");
351 printDecompositionInfo();
353 ldl_symbolic(dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
355 Lx_ =
new DofType [Lp_[dimMat]];
356 Li_ =
new int [Lp_[dimMat]];
358 const std::size_t k(ldl_numeric(dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), ccsmat_.getValues(),
359 Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
368 std::cerr<<
"LDL Error: D("<<k<<
","<<k<<
") is zero!"<<std::endl;
369 DUNE_THROW(InvalidStateException,
"LDL Error: factorisation failed!");
375template<
class DF,
class Op,
bool symmetric=false>
376using LDLOp = LDLInverseOperator< DF, typename Op::MatrixType >;
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76