DUNE-FEM (unstable)

ldl.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_LDL_HH
6#define DUNE_ISTL_LDL_HH
7
8#if HAVE_SUITESPARSE_LDL || defined DOXYGEN
9
10#include <iostream>
11#include <memory>
12#include <type_traits>
13
14#ifdef __cplusplus
15extern "C"
16{
17#include <ldl.h>
18#include <amd.h>
19}
20#endif
21
23
24#include <dune/istl/bccsmatrixinitializer.hh>
25#include <dune/istl/solvers.hh>
27#include <dune/istl/solverregistry.hh>
28
29namespace Dune {
41 // forward declarations
42 template<class M, class T, class TM, class TD, class TA>
43 class SeqOverlappingSchwarz;
44
45 template<class T, bool tag>
46 struct SeqOverlappingSchwarzAssemblerHelper;
47
54 template<class Matrix>
55 class LDL
56 {};
57
71 template<typename T, typename A, int n, int m>
72 class LDL<BCRSMatrix<FieldMatrix<T,n,m>,A > >
73 : public InverseOperator<BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >,
74 BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > >
75 {
76 public:
81 typedef Dune::ISTL::Impl::BCCSMatrix<T,int> LDLMatrix;
83 typedef ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A>, int> MatrixInitializer;
85 typedef Dune::BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > > domain_type;
87 typedef Dune::BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > range_type;
88
91 {
92 return SolverCategory::Category::sequential;
93 }
94
104 LDL(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
105 {
106 //check whether T is a supported type
107 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
108 setMatrix(matrix);
109 }
110
120 LDL(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false), verbose_(verbose)
121 {
122 //check whether T is a supported type
123 static_assert(std::is_same<T,double>::value,"Unsupported Type in LDL (only double supported)");
124 setMatrix(matrix);
125 }
126
136 LDL(const Matrix& matrix, const ParameterTree& config)
137 : LDL(matrix, config.get<int>("verbose", 0))
138 {}
139
141 LDL() : matrixIsLoaded_(false), verbose_(0)
142 {}
143
145 virtual ~LDL()
146 {
147 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
148 free();
149 }
150
153 {
154 const int dimMat(ldlMatrix_.N());
155 ldl_perm(dimMat, Y_, reinterpret_cast<double*>(&b[0]), P_);
156 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
157 ldl_dsolve(dimMat, Y_, D_);
158 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
159 ldl_permt(dimMat, reinterpret_cast<double*>(&x[0]), Y_, P_);
160 // this is a direct solver
161 res.iterations = 1;
162 res.converged = true;
163 }
164
166 void apply(domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
167 {
168 apply(x,b,res);
169 }
170
176 void apply(T* x, T* b)
177 {
178 const int dimMat(ldlMatrix_.N());
179 ldl_perm(dimMat, Y_, b, P_);
180 ldl_lsolve(dimMat, Y_, Lp_, Li_, Lx_);
181 ldl_dsolve(dimMat, Y_, D_);
182 ldl_ltsolve(dimMat, Y_, Lp_, Li_, Lx_);
183 ldl_permt(dimMat, x, Y_, P_);
184 }
185
186 void setOption([[maybe_unused]] unsigned int option, [[maybe_unused]] double value)
187 {}
188
190 void setMatrix(const Matrix& matrix)
191 {
192 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
193 free();
194
195 if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
196 ldlMatrix_.free();
197 ldlMatrix_.setSize(MatrixDimension<Matrix>::rowdim(matrix),
198 MatrixDimension<Matrix>::coldim(matrix));
199 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
200
201 copyToBCCSMatrix(initializer, matrix);
202
203 decompose();
204 }
205
206 template<class S>
207 void setSubMatrix(const Matrix& matrix, const S& rowIndexSet)
208 {
209 if ((ldlMatrix_.N() + ldlMatrix_.M() > 0) || matrixIsLoaded_)
210 free();
211
212 if (ldlMatrix_.N() + ldlMatrix_.M() + ldlMatrix_.nonzeroes() != 0)
213 ldlMatrix_.free();
214
215 ldlMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(matrix) / matrix.N(),
216 rowIndexSet.size()*MatrixDimension<Matrix>::coldim(matrix) / matrix.M());
217 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(ldlMatrix_);
218
219 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(matrix,rowIndexSet));
220
221 decompose();
222 }
223
228 inline void setVerbosity(int v)
229 {
230 verbose_=v;
231 }
232
238 {
239 return ldlMatrix_;
240 }
241
246 void free()
247 {
248 delete [] D_;
249 delete [] Y_;
250 delete [] Lp_;
251 delete [] Lx_;
252 delete [] Li_;
253 delete [] P_;
254 delete [] Pinv_;
255 ldlMatrix_.free();
256 matrixIsLoaded_ = false;
257 }
258
260 inline const char* name()
261 {
262 return "LDL";
263 }
264
269 inline double* getD()
270 {
271 return D_;
272 }
273
278 inline int* getLp()
279 {
280 return Lp_;
281 }
282
287 inline int* getLi()
288 {
289 return Li_;
290 }
291
296 inline double* getLx()
297 {
298 return Lx_;
299 }
300
301 private:
302 template<class M,class X, class TM, class TD, class T1>
303 friend class SeqOverlappingSchwarz;
304
305 friend struct SeqOverlappingSchwarzAssemblerHelper<LDL<Matrix>,true>;
306
308 void decompose()
309 {
310 // allocate vectors
311 const int dimMat(ldlMatrix_.N());
312 D_ = new double [dimMat];
313 Y_ = new double [dimMat];
314 Lp_ = new int [dimMat + 1];
315 Parent_ = new int [dimMat];
316 Lnz_ = new int [dimMat];
317 Flag_ = new int [dimMat];
318 Pattern_ = new int [dimMat];
319 P_ = new int [dimMat];
320 Pinv_ = new int [dimMat];
321
322 double Info [AMD_INFO];
323 if(amd_order (dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), P_, (double *) NULL, Info) < AMD_OK)
324 DUNE_THROW(InvalidStateException,"Error: AMD failed!");
325 if(verbose_ > 0)
326 amd_info (Info);
327 // compute the symbolic factorisation
328 ldl_symbolic(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), Lp_, Parent_, Lnz_, Flag_, P_, Pinv_);
329 // initialise those entries of additionalVectors_ whose dimension is known only now
330 Lx_ = new double [Lp_[dimMat]];
331 Li_ = new int [Lp_[dimMat]];
332 // compute the numeric factorisation
333 const int rank(ldl_numeric(dimMat, ldlMatrix_.getColStart(), ldlMatrix_.getRowIndex(), ldlMatrix_.getValues(),
334 Lp_, Parent_, Lnz_, Li_, Lx_, D_, Y_, Pattern_, Flag_, P_, Pinv_));
335 // free temporary vectors
336 delete [] Flag_;
337 delete [] Pattern_;
338 delete [] Parent_;
339 delete [] Lnz_;
340
341 if(rank!=dimMat)
342 DUNE_THROW(InvalidStateException,"Error: LDL factorisation failed!");
343 }
344
345 LDLMatrix ldlMatrix_;
346 bool matrixIsLoaded_;
347 int verbose_;
348 int* Lp_;
349 int* Parent_;
350 int* Lnz_;
351 int* Flag_;
352 int* Pattern_;
353 int* P_;
354 int* Pinv_;
355 double* D_;
356 double* Y_;
357 double* Lx_;
358 int* Li_;
359 };
360
361 template<typename T, typename A, int n, int m>
362 struct IsDirectSolver<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
363 {
364 enum {value = true};
365 };
366
367 template<typename T, typename A, int n, int m>
368 struct StoresColumnCompressed<LDL<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
369 {
370 enum {value = true};
371 };
372
373 DUNE_REGISTER_SOLVER("ldl",
374 [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
375 -> std::shared_ptr<typename decltype(opTraits)::solver_type>
376 {
377 using OpTraits = decltype(opTraits);
378 using M = typename OpTraits::matrix_type;
379 // works only for sequential operators
380 if constexpr (OpTraits::isParallel){
381 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
382 DUNE_THROW(Dune::InvalidStateException, "LDL works only for sequential operators.");
383 }
384 // check if LDL<M>* is convertible to
385 // InverseOperator*. This allows only the explicit
386 // specialized variants of LDL
387 if constexpr (std::is_convertible_v<LDL<M>*,
388 Dune::InverseOperator<typename OpTraits::domain_type,
389 typename OpTraits::range_type>*> &&
390 std::is_same_v<typename FieldTraits<M>::field_type, double>
391 ){
392 const auto& A = opTraits.getAssembledOpOrThrow(op);
393 const M& mat = A->getmat();
394 int verbose = config.get("verbose", 0);
395 return std::make_shared<LDL<M>>(mat,verbose);
396 }
397 DUNE_THROW(UnsupportedType,
398 "Unsupported Type in LDL (only FieldMatrix<double,...> supported)");
399 });
400
401} // end namespace Dune
402
403
404#endif //HAVE_SUITESPARSE_LDL
405#endif //DUNE_ISTL_LDL_HH
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:466
A vector of blocks with memory management.
Definition: bvector.hh:392
A dense n x m matrix.
Definition: fmatrix.hh:117
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Abstract base class for all solvers.
Definition: solver.hh:101
Use the LDL package to directly solve linear systems – empty default class.
Definition: ldl.hh:56
A generic dynamic dense matrix.
Definition: matrix.hh:561
size_type M() const
Return the number of columns.
Definition: matrix.hh:696
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:188
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:755
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
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
int * getLp()
Get factorization Lp.
Definition: ldl.hh:278
LDLMatrix & getInternalMatrix()
Return the column compress matrix.
Definition: ldl.hh:237
void apply(domain_type &x, range_type &b, InverseOperatorResult &res) override
Apply inverse operator,.
Definition: ldl.hh:152
Dune::ISTL::Impl::BCCSMatrix< T, int > LDLMatrix
The corresponding SuperLU Matrix type.
Definition: ldl.hh:81
double * getLx()
Get factorization Lx.
Definition: ldl.hh:296
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition: ldl.hh:228
void apply(T *x, T *b)
Additional apply method with c-arrays in analogy to superlu.
Definition: ldl.hh:176
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
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res) override
apply inverse operator, with given convergence criteria.
Definition: ldl.hh:166
virtual ~LDL()
Default constructor.
Definition: ldl.hh:145
void free()
Free allocated space.
Definition: ldl.hh:246
int * getLi()
Get factorization Li.
Definition: ldl.hh:287
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: ldl.hh:190
LDL(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: ldl.hh:120
double * getD()
Get factorization diagonal matrix D.
Definition: ldl.hh:269
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category)
Definition: ldl.hh:90
LDL(const Matrix &matrix, const ParameterTree &config)
Constructs the LDL solver.
Definition: ldl.hh:136
ISTL::Impl::BCCSMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A >, int > MatrixInitializer
Type of an associated initializer class.
Definition: ldl.hh:83
const char * name()
Get method name.
Definition: ldl.hh:260
LDL()
Default constructor.
Definition: ldl.hh:141
LDL(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition: ldl.hh:104
Dune namespace.
Definition: alignedallocator.hh:13
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Statistics about the application of an inverse operator.
Definition: solver.hh:50
int iterations
Number of iterations.
Definition: solver.hh:69
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Category
Definition: solvercategory.hh:23
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)