Dune Core Modules (2.8.0)

ldl.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_LDL_HH
4#define DUNE_ISTL_LDL_HH
5
6#if HAVE_SUITESPARSE_LDL || defined DOXYGEN
7
8#include <iostream>
9#include <memory>
10#include <type_traits>
11
12#ifdef __cplusplus
13extern "C"
14{
15#include "ldl.h"
16#include "amd.h"
17}
18#endif
19
21
22#include <dune/istl/bccsmatrixinitializer.hh>
23#include <dune/istl/solvers.hh>
25#include <dune/istl/solverfactory.hh>
26
27namespace Dune {
39 // forward declarations
40 template<class M, class T, class TM, class TD, class TA>
41 class SeqOverlappingSchwarz;
42
43 template<class T, bool tag>
44 struct SeqOverlappingSchwarzAssemblerHelper;
45
52 template<class Matrix>
53 class LDL
54 {};
55
69 template<typename T, typename A, int n, int m>
70 class LDL<BCRSMatrix<FieldMatrix<T,n,m>,A > >
71 : public InverseOperator<BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >,
72 BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > >
73 {
74 public:
79 typedef Dune::ISTL::Impl::BCCSMatrix<T,int> LDLMatrix;
81 typedef ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A>, int> MatrixInitializer;
83 typedef Dune::BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > > domain_type;
85 typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;
86
89 {
90 return SolverCategory::Category::sequential;
91 }
92
102 LDL(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
103 {
104 //check whether T is a supported type
105 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
106 setMatrix(matrix);
107 }
108
118 LDL(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false), verbose_(verbose)
119 {
120 //check whether T is a supported type
121 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
122 setMatrix(matrix);
123 }
124
134 LDL(const Matrix& matrix, const ParameterTree& config)
135 : LDL(matrix, config.get<int>("verbose", 0))
136 {}
137
139 LDL() : matrixIsLoaded_(false), verbose_(0)
140 {}
141
143 virtual ~LDL()
144 {
145 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
146 free();
147 }
148
151 {
152 const int dimMat(ldlMatrix_.N());
153 ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
154 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
155 ldl_dsolve(dimMat, Y_, D_);
156 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
157 ldl_permt(dimMat, reinterpret_cast<double*>(&x[0]), Y_, P_);
158 // this is a direct solver
159 res.iterations = 1;
160 res.converged = true;
161 }
162
164 virtual void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
165 {
166 apply(x,b,res);
167 }
168
174 void apply(T* x, T* b)
175 {
176 const int dimMat(ldlMatrix_.N());
177 ldl_perm(dimMat, Y_, b, P_);
178 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
179 ldl_dsolve(dimMat, Y_, D_);
180 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
181 ldl_permt(dimMat, x, Y_, P_);
182 }
183
184 void setOption([[maybe_unused]] unsigned int option, [[maybe_unused]] double value)
185 {}
186
188 void setMatrix(const Matrix& matrix)
189 {
190 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
191 free();
192
193 if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
194 ldlMatrix_.free();
195 ldlMatrix_.setSize(MatrixDimension<Matrix>::rowdim(matrix),
196 MatrixDimension<Matrix>::coldim(matrix));
197 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
198
199 copyToBCCSMatrix(initializer, matrix);
200
201 decompose();
202 }
203
204 template<class S>
205 void setSubMatrix(const Matrix& matrix, const S& rowIndexSet)
206 {
207 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
208 free();
209
210 if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
211 ldlMatrix_.free();
212
213 ldlMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(matrix) / matrix.N(),
214 rowIndexSet.size()*MatrixDimension<Matrix>::coldim(matrix) / matrix.M());
215 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
216
217 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(matrix,rowIndexSet));
218
219 decompose();
220 }
221
226 inline void setVerbosity(int v)
227 {
228 verbose_=v;
229 }
230
236 {
237 return ldlMatrix_;
238 }
239
244 void free()
245 {
246 delete [] D_;
247 delete [] Y_;
248 delete [] Lp_;
249 delete [] Lx_;
250 delete [] Li_;
251 delete [] P_;
252 delete [] Pinv_;
253 ldlMatrix_.free();
254 matrixIsLoaded_ = false;
255 }
256
258 inline const char* name()
259 {
260 return "LDL";
261 }
262
267 inline double* getD()
268 {
269 return D_;
270 }
271
276 inline int* getLp()
277 {
278 return Lp_;
279 }
280
285 inline int* getLi()
286 {
287 return Li_;
288 }
289
294 inline double* getLx()
295 {
296 return Lx_;
297 }
298
299 private:
300 template<class M,class X, class TM, class TD, class T1>
301 friend class SeqOverlappingSchwarz;
302
303 friend struct SeqOverlappingSchwarzAssemblerHelper<LDL<Matrix>,true>;
304
306 void decompose()
307 {
308 // allocate vectors
309 const int dimMat(ldlMatrix_.N());
310 D_ = new double [dimMat];
311 Y_ = new double [dimMat];
312 Lp_ = new int [dimMat + 1];
313 Parent_ = new int [dimMat];
314 Lnz_ = new int [dimMat];
315 Flag_ = new int [dimMat];
316 Pattern_ = new int [dimMat];
317 P_ = new int [dimMat];
318 Pinv_ = new int [dimMat];
319
320 double Info [AMD_INFO];
321 if(amd_order (dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), P_, (double *) NULL, Info) < AMD_OK)
322 DUNE_THROW(InvalidStateException,"Error: AMD failed!");
323 if(verbose_ > 0)
324 amd_info (Info);
325 // compute the symbolic factorisation
326 ldl_symbolic(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
327 // initialise those entries of additionalVectors_ whose dimension is known only now
328 Lx_ = new double [Lp_[dimMat]];
329 Li_ = new int [Lp_[dimMat]];
330 // compute the numeric factorisation
331 const int rank(ldl_numeric(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), ldlMatrix_.getValues(),
332 Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
333 // free temporary vectors
334 delete [] Flag_;
335 delete [] Pattern_;
336 delete [] Parent_;
337 delete [] Lnz_;
338
339 if(rank!=dimMat)
340 DUNE_THROW(InvalidStateException,"Error: LDL factorisation failed!");
341 }
342
343 LDLMatrix ldlMatrix_;
344 bool matrixIsLoaded_;
345 int verbose_;
346 int* Lp_;
347 int* Parent_;
348 int* Lnz_;
349 int* Flag_;
350 int* Pattern_;
351 int* P_;
352 int* Pinv_;
353 double* D_;
354 double* Y_;
355 double* Lx_;
356 int* Li_;
357 };
358
359 template<typename T, typename A, int n, int m>
360 struct IsDirectSolver<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
361 {
362 enum {value = true};
363 };
364
365 template<typename T, typename A, int n, int m>
366 struct StoresColumnCompressed<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
367 {
368 enum {value = true};
369 };
370
371 struct LDLCreator {
372 template<class F> struct isValidBlock : std::false_type{};
373 template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
374
375 template<typename TL, typename M>
376 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
377 typename Dune::TypeListElement<2, TL>::type>>
378 operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
379 std::enable_if_t<
380 isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
381 {
382 int verbose = config.get("verbose", 0);
383 return std::make_shared<Dune::LDL<M>>(mat,verbose);
384 }
385
386 // second version with SFINAE to validate the template parameters of LDL
387 template<typename TL, typename M>
388 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
389 typename Dune::TypeListElement<2, TL>::type>>
390 operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
391 std::enable_if_t<
392 !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
393 {
394 DUNE_THROW(UnsupportedType,
395 "Unsupported Type in LDL (only double and std::complex<double> supported)");
396 }
397 };
398 DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator());
399
400} // end namespace Dune
401
402
403#endif //HAVE_SUITESPARSE_LDL
404#endif //DUNE_ISTL_LDL_HH
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
A vector of blocks with memory management.
Definition: bvector.hh:393
A dense n x m matrix.
Definition: fmatrix.hh:69
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:279
Abstract base class for all solvers.
Definition: solver.hh:97
Use the LDL package to directly solve linear systems – empty default class.
Definition: ldl.hh:54
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
Hierarchical structure of string parameters.
Definition: parametertree.hh:35
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:183
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
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:83
int * getLp()
Get factorization Lp.
Definition: ldl.hh:276
LDLMatrix & getInternalMatrix()
Return the column compress matrix.
Definition: ldl.hh:235
Dune::ISTL::Impl::BCCSMatrix< T, int > LDLMatrix
The corresponding SuperLU Matrix type.
Definition: ldl.hh:79
double * getLx()
Get factorization Lx.
Definition: ldl.hh:294
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition: ldl.hh:226
void apply(T *x, T *b)
Additional apply method with c-arrays in analogy to superlu.
Definition: ldl.hh:174
virtual void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: ldl.hh:164
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:85
virtual ~LDL()
Default constructor.
Definition: ldl.hh:143
void free()
Free allocated space.
Definition: ldl.hh:244
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: ldl.hh:150
int * getLi()
Get factorization Li.
Definition: ldl.hh:285
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: ldl.hh:188
LDL(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: ldl.hh:118
double * getD()
Get factorization diagonal matrix D.
Definition: ldl.hh:267
LDL(const Matrix &matrix, const ParameterTree &config)
Constructs the LDL solver.
Definition: ldl.hh:134
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: ldl.hh:88
ISTL::Impl::BCCSMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A >, int > MatrixInitializer
Type of an associated initializer class.
Definition: ldl.hh:81
const char * name()
Get method name.
Definition: ldl.hh:258
LDL()
Default constructor.
Definition: ldl.hh:139
LDL(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition: ldl.hh:102
Dune namespace.
Definition: alignedallocator.hh:11
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Statistics about the application of an inverse operator.
Definition: solver.hh:46
int iterations
Number of iterations.
Definition: solver.hh:65
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Category
Definition: solvercategory.hh:21
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)