Dune Core Modules (2.5.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 <type_traits>
10
11#ifdef __cplusplus
12extern "C"
13{
14#include "ldl.h"
15#include "amd.h"
16}
17#endif
18
20#include <dune/common/unused.hh>
21
22#include <dune/istl/colcompmatrix.hh>
23#include <dune/istl/solvers.hh>
25
26namespace Dune {
38 // forward declarations
39 template<class M, class T, class TM, class TD, class TA>
40 class SeqOverlappingSchwarz;
41
42 template<class T, bool tag>
43 struct SeqOverlappingSchwarzAssemblerHelper;
44
51 template<class Matrix>
52 class LDL
53 {};
54
68 template<typename T, typename A, int n, int m>
69 class LDL<BCRSMatrix<FieldMatrix<T,n,m>,A > >
70 : public InverseOperator<BlockVector<FieldVector<T,m>, typename A::template rebind<FieldVector<T,m> >::other>,
71 BlockVector<FieldVector<T,n>, typename A::template rebind<FieldVector<T,n> >::other> >
72 {
73 public:
82 typedef Dune::BlockVector<FieldVector<T,m>, typename A::template rebind<FieldVector<T,m> >::other> domain_type;
84 typedef Dune::BlockVector<FieldVector<T,n>, typename A::template rebind<FieldVector<T,n> >::other> range_type;
85
95 LDL(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
96 {
97 //check whether T is a supported type
98 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
99 setMatrix(matrix);
100 }
101
111 LDL(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false), verbose_(verbose)
112 {
113 //check whether T is a supported type
114 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
115 setMatrix(matrix);
116 }
117
119 LDL() : matrixIsLoaded_(false), verbose_(0)
120 {}
121
123 virtual ~LDL()
124 {
125 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
126 free();
127 }
128
131 {
132 const int dimMat(ldlMatrix_.N());
133 ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
134 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
135 ldl_dsolve(dimMat, Y_, D_);
136 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
137 ldl_permt(dimMat, reinterpret_cast<double*>(&x[0]), Y_, P_);
138 // this is a direct solver
139 res.iterations = 1;
140 res.converged = true;
141 }
142
144 virtual void apply(domain_type& x, range_type& b, double reduction, InverseOperatorResult& res)
145 {
146 DUNE_UNUSED_PARAMETER(reduction);
147 apply(x,b,res);
148 }
149
155 void apply(T* x, T* b)
156 {
157 const int dimMat(ldlMatrix_.N());
158 ldl_perm(dimMat, Y_, b, P_);
159 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
160 ldl_dsolve(dimMat, Y_, D_);
161 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
162 ldl_permt(dimMat, x, Y_, P_);
163 }
164
165 void setOption(unsigned int option, double value)
166 {
167 DUNE_UNUSED_PARAMETER(option);
169 }
170
172 void setMatrix(const Matrix& matrix)
173 {
174 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
175 free();
176 ldlMatrix_ = matrix;
177 decompose();
178 }
179
180 template<class S>
181 void setSubMatrix(const Matrix& matrix, const S& rowIndexSet)
182 {
183 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
184 free();
185 ldlMatrix_.setMatrix(matrix,rowIndexSet);
186 decompose();
187 }
188
193 inline void setVerbosity(int v)
194 {
195 verbose_=v;
196 }
197
203 {
204 return ldlMatrix_;
205 }
206
211 void free()
212 {
213 delete [] D_;
214 delete [] Y_;
215 delete [] Lp_;
216 delete [] Lx_;
217 delete [] Li_;
218 delete [] P_;
219 delete [] Pinv_;
220 ldlMatrix_.free();
221 matrixIsLoaded_ = false;
222 }
223
225 inline const char* name()
226 {
227 return "LDL";
228 }
229
234 inline double* getD()
235 {
236 return D_;
237 }
238
243 inline int* getLp()
244 {
245 return Lp_;
246 }
247
252 inline int* getLi()
253 {
254 return Li_;
255 }
256
261 inline double* getLx()
262 {
263 return Lx_;
264 }
265
266 private:
267 template<class M,class X, class TM, class TD, class T1>
268 friend class SeqOverlappingSchwarz;
269
270 friend struct SeqOverlappingSchwarzAssemblerHelper<LDL<Matrix>,true>;
271
273 void decompose()
274 {
275 // allocate vectors
276 const int dimMat(ldlMatrix_.N());
277 D_ = new double [dimMat];
278 Y_ = new double [dimMat];
279 Lp_ = new int [dimMat + 1];
280 Parent_ = new int [dimMat];
281 Lnz_ = new int [dimMat];
282 Flag_ = new int [dimMat];
283 Pattern_ = new int [dimMat];
284 P_ = new int [dimMat];
285 Pinv_ = new int [dimMat];
286
287 double Info [AMD_INFO];
288 if(amd_order (dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), P_, (double *) NULL, Info) < AMD_OK)
289 DUNE_THROW(InvalidStateException,"Error: AMD failed!");
290 if(verbose_ > 0)
291 amd_info (Info);
292 // compute the symbolic factorisation
293 ldl_symbolic(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
294 // initialise those entries of additionalVectors_ whose dimension is known only now
295 Lx_ = new double [Lp_[dimMat]];
296 Li_ = new int [Lp_[dimMat]];
297 // compute the numeric factorisation
298 const int rank(ldl_numeric(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), ldlMatrix_.getValues(),
299 Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
300 // free temporary vectors
301 delete [] Flag_;
302 delete [] Pattern_;
303 delete [] Parent_;
304 delete [] Lnz_;
305
306 if(rank!=dimMat)
307 DUNE_THROW(InvalidStateException,"Error: LDL factorisation failed!");
308 }
309
310 LDLMatrix ldlMatrix_;
311 bool matrixIsLoaded_;
312 int verbose_;
313 int* Lp_;
314 int* Parent_;
315 int* Lnz_;
316 int* Flag_;
317 int* Pattern_;
318 int* P_;
319 int* Pinv_;
320 double* D_;
321 double* Y_;
322 double* Lx_;
323 int* Li_;
324 };
325
326 template<typename T, typename A, int n, int m>
327 struct IsDirectSolver<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
328 {
329 enum {value = true};
330 };
331
332 template<typename T, typename A, int n, int m>
333 struct StoresColumnCompressed<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
334 {
335 enum {value = true};
336 };
337
338}
339
340#endif //HAVE_SUITESPARSE_LDL
341#endif //DUNE_ISTL_LDL_HH
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:423
A vector of blocks with memory management.
Definition: bvector.hh:313
A dense n x m matrix.
Definition: fmatrix.hh:68
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:79
Use the LDL package to directly solve linear systems – empty default class.
Definition: ldl.hh:53
A generic dynamic dense matrix.
Definition: matrix.hh:555
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:742
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
int * getLp()
Get factorization Lp.
Definition: ldl.hh:243
LDLMatrix & getInternalMatrix()
Return the column compress matrix.
Definition: ldl.hh:202
Dune::BlockVector< FieldVector< T, n >, typename A::template rebind< FieldVector< T, n > >::other > range_type
The type of the range of the solver.
Definition: ldl.hh:84
double * getLx()
Get factorization Lx.
Definition: ldl.hh:261
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition: ldl.hh:193
void apply(T *x, T *b)
Additional apply method with c-arrays in analogy to superlu.
Definition: ldl.hh:155
ColCompMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A > > MatrixInitializer
Type of an associated initializer class.
Definition: ldl.hh:80
virtual void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: ldl.hh:144
virtual ~LDL()
Default constructor.
Definition: ldl.hh:123
Dune::BlockVector< FieldVector< T, m >, typename A::template rebind< FieldVector< T, m > >::other > domain_type
The type of the domain of the solver.
Definition: ldl.hh:82
void free()
Free allocated space.
Definition: ldl.hh:211
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: ldl.hh:130
int * getLi()
Get factorization Li.
Definition: ldl.hh:252
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: ldl.hh:172
LDL(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: ldl.hh:111
double * getD()
Get factorization diagonal matrix D.
Definition: ldl.hh:234
Dune::ColCompMatrix< Matrix > LDLMatrix
The corresponding SuperLU Matrix type.
Definition: ldl.hh:78
const char * name()
Get method name.
Definition: ldl.hh:225
LDL()
Default constructor.
Definition: ldl.hh:119
LDL(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition: ldl.hh:95
Dune namespace.
Definition: alignment.hh:11
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Inititializer for the ColCompMatrix as needed by OverlappingSchwarz.
Definition: colcompmatrix.hh:154
Statistics about the application of an inverse operator.
Definition: solver.hh:32
int iterations
Number of iterations.
Definition: solver.hh:50
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)