dune-fem 2.8.0
Loading...
Searching...
No Matches
ldlsolver.hh
Go to the documentation of this file.
1#ifndef DUNE_FEM_LDLSOLVER_HH
2#define DUNE_FEM_LDLSOLVER_HH
3
4#include <algorithm>
5#include <iostream>
6#include <tuple>
7#include <vector>
8#include <utility>
9
10#include <dune/common/exceptions.hh>
11#include <dune/common/hybridutilities.hh>
19
20#if HAVE_SUITESPARSE_LDL
21
22#ifdef __cplusplus
23extern "C"
24{
25#include "ldl.h"
26#include "amd.h"
27}
28#endif
29
30namespace Dune
31{
32namespace Fem
33{
34
45template< class DiscreteFunction,
46 class Matrix = SparseRowMatrix< typename DiscreteFunction::DiscreteFunctionSpaceType::RangeFieldType > >
47class LDLInverseOperator;
48
49template< class DiscreteFunction, class Matrix >
50struct LDLInverseOperatorTraits
51{
52 typedef DiscreteFunction DiscreteFunctionType;
53 typedef AdaptiveDiscreteFunction< typename DiscreteFunction::DiscreteFunctionSpaceType > SolverDiscreteFunctionType;
54
56 typedef OperatorType PreconditionerType;
57
58 typedef Fem::SparseRowLinearOperator< DiscreteFunction, DiscreteFunction, Matrix > AssembledOperatorType;
59
60 typedef LDLInverseOperator< DiscreteFunction, Matrix > InverseOperatorType;
61 typedef SolverParameter SolverParameterType;
62};
63
64
65
73template< class DF, class Matrix >
74class LDLInverseOperator : public InverseOperatorInterface< LDLInverseOperatorTraits< DF, Matrix > >
75{
76 typedef LDLInverseOperatorTraits< DF, Matrix > Traits;
77 typedef InverseOperatorInterface< Traits > BaseType;
78
79 friend class InverseOperatorInterface< Traits >;
80public:
81 typedef LDLInverseOperator< DF, Matrix > ThisType;
82
83 typedef typename BaseType :: SolverDiscreteFunctionType
84 SolverDiscreteFunctionType;
85
86 typedef typename BaseType :: OperatorType OperatorType;
87 typedef typename BaseType :: AssembledOperatorType AssembledOperatorType;
88
89 // \brief The column-compressed matrix type.
90 typedef ColCompMatrix<typename AssembledOperatorType::MatrixType::MatrixBaseType,int> CCSMatrixType;
91
92 typedef typename SolverDiscreteFunctionType::DofType DofType;
93 typedef typename SolverDiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
94
101 LDLInverseOperator(const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
102 const ParameterReader &parameter = Parameter::container() ) :
103 verbose_(verbose && Dune::Fem::Parameter::verbose()), ccsmat_()
104 {}
105
111 LDLInverseOperator(const double& redEps, const double& absLimit, const int& maxIter,
112 const ParameterReader &parameter = Parameter::container() ) :
113 verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_()
114 {}
115
116 LDLInverseOperator(const SolverParameter &parameter = SolverParameter(Parameter::container()) )
117 : BaseType(parameter), verbose_(BaseType::verbose())
118 {}
119
127 LDLInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter, const bool& verbose,
128 const ParameterReader &parameter = Parameter::container() ) :
129 verbose_(verbose), ccsmat_(), isloaded_(false)
130 {
131 bind(op);
132 }
133
140 LDLInverseOperator(const OperatorType& op, const double& redEps, const double& absLimit, const int& maxIter,
141 const ParameterReader &parameter = Parameter::container() ) :
142 verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
143 {
144 bind(op);
145 }
146
147 LDLInverseOperator(const OperatorType& op, const ParameterReader &parameter = Parameter::container() ) :
148 verbose_(parameter.getValue<bool>("fem.solver.verbose",false)), ccsmat_(), isloaded_(false)
149 {
150 bind(op);
151 }
152
153 // \brief Destructor.
154 ~LDLInverseOperator()
155 {
156 unbind();
157 }
158
159 void bind( const OperatorType& op )
160 {
161 // clear old storage
162 finalize();
163 BaseType::bind( op );
164 }
165
166 void unbind ()
167 {
168 finalize();
169 BaseType::unbind();
170 }
171
172 // \brief Decompose matrix.
173 template<typename... A>
174 void prepare(A... ) const
175 {
176 if( ! assembledOperator_ )
177 DUNE_THROW(NotImplemented,"LDLInverseOperator only works for assembled systems!");
178
179 if(!isloaded_)
180 {
181 ccsmat_ = assembledOperator_->exportMatrix();
182 decompose();
183 isloaded_ = true;
184 }
185 }
186
187 // \brief Free allocated memory.
188 virtual void finalize()
189 {
190 if(isloaded_)
191 {
192 ccsmat_.free();
193 delete [] D_;
194 delete [] Y_;
195 delete [] Lp_;
196 delete [] Lx_;
197 delete [] Li_;
198 delete [] P_;
199 delete [] Pinv_;
200 isloaded_ = false;
201 }
202 }
203
210 template<typename... DFs>
211 void apply(const TupleDiscreteFunction<DFs...>& arg,TupleDiscreteFunction<DFs...>& dest) const
212 {
213 // copy DOF's arg into a consecutive vector
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());
219 // apply operator
220 apply(vecArg.data(),vecDest.data());
221 // copy back solution into dest
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++));});
224 }
225
226protected:
227 void printTexInfo(std::ostream& out) const
228 {
229 out<<"Solver: LDL direct solver"<<std::endl;
230 }
231
232 // \brief Print some statistics about the LDL decomposition.
233 void printDecompositionInfo() const
234 {
235 amd_info(info_);
236 }
237
241 DofType* getD()
242 {
243 return D_;
244 }
245
249 int* getLp()
250 {
251 return Lp_;
252 }
253
257 int* getLi()
258 {
259 return Li_;
260 }
261
265 DofType* getLx()
266 {
267 return Lx_;
268 }
269
273 CCSMatrixType& getCCSMatrix()
274 {
275 return ccsmat_;
276 }
277
278protected:
285 void apply(const DofType* arg, DofType* dest) const
286 {
287 prepare();
288
289 // apply part of the call
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_);
296
297 const_cast<ThisType*>(this)->finalize();
298 }
299
306 int apply(const SolverDiscreteFunctionType& arg,
307 SolverDiscreteFunctionType& dest) const
308 {
309 apply(arg.leakPointer(),dest.leakPointer());
310 return 0;
311 }
312
313
314protected:
315 using BaseType :: assembledOperator_;
316 const bool verbose_;
317
318 mutable CCSMatrixType ccsmat_;
319 mutable bool isloaded_ = false;
320 mutable int* Lp_;
321 mutable int* Parent_;
322 mutable int* Lnz_;
323 mutable int* Flag_;
324 mutable int* Pattern_;
325 mutable int* P_;
326 mutable int* Pinv_;
327 mutable DofType* D_;
328 mutable DofType* Y_;
329 mutable DofType* Lx_;
330 mutable int* Li_;
331 mutable double info_[AMD_INFO];
332
333 // \brief Computes the LDL decomposition.
334 void decompose() const
335 {
336 // allocate vectors
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];
347
348 if(amd_order (dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), P_, (DofType *) NULL, info_) < AMD_OK)
349 DUNE_THROW(InvalidStateException,"LDL Error: AMD failed!");
350 if(verbose_)
351 printDecompositionInfo();
352 // compute the symbolic factorisation
353 ldl_symbolic(dimMat, ccsmat_.getColStart(), ccsmat_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
354 // initialise those entries of additionalVectors_ whose dimension is known only now
355 Lx_ = new DofType [Lp_[dimMat]];
356 Li_ = new int [Lp_[dimMat]];
357 // compute the numeric factorisation
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_));
360 // free temporary vectors
361 delete [] Flag_;
362 delete [] Pattern_;
363 delete [] Parent_;
364 delete [] Lnz_;
365
366 if(k!=dimMat)
367 {
368 std::cerr<<"LDL Error: D("<<k<<","<<k<<") is zero!"<<std::endl;
369 DUNE_THROW(InvalidStateException,"LDL Error: factorisation failed!");
370 }
371 }
372};
373
374// deprecated old type
375template<class DF, class Op, bool symmetric=false>
376using LDLOp = LDLInverseOperator< DF, typename Op::MatrixType >;
377
378} // end namespace Fem
379} // end namespace Dune
380
381#endif // #if HAVE_SUITESPARSE_LDL
382
383#endif // #ifndef DUNE_FEM_LDLSOLVER_HH
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76