Dune Core Modules (2.9.1)

umfpack.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) 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_UMFPACK_HH
6#define DUNE_ISTL_UMFPACK_HH
7
8#if HAVE_SUITESPARSE_UMFPACK || defined DOXYGEN
9
10#include<complex>
11#include<type_traits>
12
13#include<umfpack.h>
14
18#include<dune/istl/bccsmatrixinitializer.hh>
22#include <dune/istl/solverfactory.hh>
23
24
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
45 // wrapper class for C-Function Calls in the backend. Choose the right function namespace
46 // depending on the template parameter used.
47 template<typename T>
48 struct UMFPackMethodChooser
49 {
50 static constexpr bool valid = false ;
51 };
52
53 template<>
54 struct UMFPackMethodChooser<double>
55 {
56 static constexpr bool valid = true ;
57
58 template<typename... A>
59 static void defaults(A... args)
60 {
61 umfpack_dl_defaults(args...);
62 }
63 template<typename... A>
64 static void free_numeric(A... args)
65 {
66 umfpack_dl_free_numeric(args...);
67 }
68 template<typename... A>
69 static void free_symbolic(A... args)
70 {
71 umfpack_dl_free_symbolic(args...);
72 }
73 template<typename... A>
74 static int load_numeric(A... args)
75 {
76 return umfpack_dl_load_numeric(args...);
77 }
78 template<typename... A>
79 static void numeric(A... args)
80 {
81 umfpack_dl_numeric(args...);
82 }
83 template<typename... A>
84 static void report_info(A... args)
85 {
86 umfpack_dl_report_info(args...);
87 }
88 template<typename... A>
89 static void report_status(A... args)
90 {
91 umfpack_dl_report_status(args...);
92 }
93 template<typename... A>
94 static int save_numeric(A... args)
95 {
96 return umfpack_dl_save_numeric(args...);
97 }
98 template<typename... A>
99 static void solve(A... args)
100 {
101 umfpack_dl_solve(args...);
102 }
103 template<typename... A>
104 static void symbolic(A... args)
105 {
106 umfpack_dl_symbolic(args...);
107 }
108 };
109
110 template<>
111 struct UMFPackMethodChooser<std::complex<double> >
112 {
113 static constexpr bool valid = true ;
114 using size_type = SuiteSparse_long;
115
116 template<typename... A>
117 static void defaults(A... args)
118 {
119 umfpack_zl_defaults(args...);
120 }
121 template<typename... A>
122 static void free_numeric(A... args)
123 {
124 umfpack_zl_free_numeric(args...);
125 }
126 template<typename... A>
127 static void free_symbolic(A... args)
128 {
129 umfpack_zl_free_symbolic(args...);
130 }
131 template<typename... A>
132 static int load_numeric(A... args)
133 {
134 return umfpack_zl_load_numeric(args...);
135 }
136 template<typename... A>
137 static void numeric(const size_type* cs, const size_type* ri, const double* val, A... args)
138 {
139 umfpack_zl_numeric(cs,ri,val,NULL,args...);
140 }
141 template<typename... A>
142 static void report_info(A... args)
143 {
144 umfpack_zl_report_info(args...);
145 }
146 template<typename... A>
147 static void report_status(A... args)
148 {
149 umfpack_zl_report_status(args...);
150 }
151 template<typename... A>
152 static int save_numeric(A... args)
153 {
154 return umfpack_zl_save_numeric(args...);
155 }
156 template<typename... A>
157 static void solve(size_type m, const size_type* cs, const size_type* ri, std::complex<double>* val, double* x, const double* b,A... args)
158 {
159 const double* cval = reinterpret_cast<const double*>(val);
160 umfpack_zl_solve(m,cs,ri,cval,NULL,x,NULL,b,NULL,args...);
161 }
162 template<typename... A>
163 static void symbolic(size_type m, size_type n, const size_type* cs, const size_type* ri, const double* val, A... args)
164 {
165 umfpack_zl_symbolic(m,n,cs,ri,val,NULL,args...);
166 }
167 };
168
169 namespace Impl
170 {
171 template<class M>
172 struct UMFPackVectorChooser
173 {};
174
175 template<typename T, typename A, int n, int m>
176 struct UMFPackVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
177 {
179 using domain_type = BlockVector<
180 FieldVector<T,m>,
181 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >;
183 using range_type = BlockVector<
184 FieldVector<T,n>,
185 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > >;
186 };
187
188 template<typename T, typename A>
189 struct UMFPackVectorChooser<BCRSMatrix<T,A> >
190 {
192 using domain_type = BlockVector<T, A>;
194 using range_type = BlockVector<T, A>;
195 };
196 }
197
211 template<typename M>
213 : public InverseOperator<
214 typename Impl::UMFPackVectorChooser<M>::domain_type,
215 typename Impl::UMFPackVectorChooser<M>::range_type >
216 {
217 using T = typename M::field_type;
218
219
220 public:
221 using size_type = SuiteSparse_long;
222
224 using Matrix = M;
225 using matrix_type = M;
227 using UMFPackMatrix = ISTL::Impl::BCCSMatrix<typename Matrix::field_type, size_type>;
229 using MatrixInitializer = ISTL::Impl::BCCSMatrixInitializer<M, size_type>;
231 using domain_type = typename Impl::UMFPackVectorChooser<M>::domain_type;
233 using range_type = typename Impl::UMFPackVectorChooser<M>::range_type;
234
237 {
238 return SolverCategory::Category::sequential;
239 }
240
249 UMFPack(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false)
250 {
251 //check whether T is a supported type
252 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
253 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
254 Caller::defaults(UMF_Control);
255 setVerbosity(verbose);
256 setMatrix(matrix);
257 }
258
267 UMFPack(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false)
268 {
269 //check whether T is a supported type
270 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
271 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
272 Caller::defaults(UMF_Control);
273 setVerbosity(verbose);
274 setMatrix(matrix);
275 }
276
286 UMFPack(const Matrix& mat_, const ParameterTree& config)
287 : UMFPack(mat_, config.get<int>("verbose", 0))
288 {}
289
292 UMFPack() : matrixIsLoaded_(false), verbosity_(0)
293 {
294 //check whether T is a supported type
295 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
296 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
297 Caller::defaults(UMF_Control);
298 }
299
310 UMFPack(const Matrix& mat_, const char* file, int verbose=0)
311 {
312 //check whether T is a supported type
313 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
314 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
315 Caller::defaults(UMF_Control);
316 setVerbosity(verbose);
317 int errcode = Caller::load_numeric(&UMF_Numeric, const_cast<char*>(file));
318 if ((errcode == UMFPACK_ERROR_out_of_memory) || (errcode == UMFPACK_ERROR_file_IO))
319 {
320 matrixIsLoaded_ = false;
321 setMatrix(mat_);
322 saveDecomposition(file);
323 }
324 else
325 {
326 matrixIsLoaded_ = true;
327 std::cout << "UMFPack decomposition successfully loaded from " << file << std::endl;
328 }
329 }
330
337 UMFPack(const char* file, int verbose=0)
338 {
339 //check whether T is a supported type
340 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
341 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
342 Caller::defaults(UMF_Control);
343 int errcode = Caller::load_numeric(&UMF_Numeric, const_cast<char*>(file));
344 if (errcode == UMFPACK_ERROR_out_of_memory)
345 DUNE_THROW(Dune::Exception, "ran out of memory while loading UMFPack decomposition");
346 if (errcode == UMFPACK_ERROR_file_IO)
347 DUNE_THROW(Dune::Exception, "IO error while loading UMFPack decomposition");
348 matrixIsLoaded_ = true;
349 std::cout << "UMFPack decomposition successfully loaded from " << file << std::endl;
350 setVerbosity(verbose);
351 }
352
353 virtual ~UMFPack()
354 {
355 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
356 free();
357 }
358
363 {
364 if (umfpackMatrix_.N() != b.dim())
365 DUNE_THROW(Dune::ISTLError, "Size of right-hand-side vector b does not match the number of matrix rows!");
366 if (umfpackMatrix_.M() != x.dim())
367 DUNE_THROW(Dune::ISTLError, "Size of solution vector x does not match the number of matrix columns!");
368 if (b.size() == 0)
369 return;
370
371 double UMF_Apply_Info[UMFPACK_INFO];
372 Caller::solve(UMFPACK_A,
373 umfpackMatrix_.getColStart(),
374 umfpackMatrix_.getRowIndex(),
375 umfpackMatrix_.getValues(),
376 reinterpret_cast<double*>(&x[0]),
377 reinterpret_cast<double*>(&b[0]),
378 UMF_Numeric,
379 UMF_Control,
380 UMF_Apply_Info);
381
382 //this is a direct solver
383 res.iterations = 1;
384 res.converged = true;
385 res.elapsed = UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME];
386
387 printOnApply(UMF_Apply_Info);
388 }
389
393 virtual void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
394 {
395 apply(x,b,res);
396 }
397
403 void apply(T* x, T* b)
404 {
405 double UMF_Apply_Info[UMFPACK_INFO];
406 Caller::solve(UMFPACK_A,
407 umfpackMatrix_.getColStart(),
408 umfpackMatrix_.getRowIndex(),
409 umfpackMatrix_.getValues(),
410 x,
411 b,
412 UMF_Numeric,
413 UMF_Control,
414 UMF_Apply_Info);
415 printOnApply(UMF_Apply_Info);
416 }
417
429 void setOption(unsigned int option, double value)
430 {
431 if (option >= UMFPACK_CONTROL)
432 DUNE_THROW(RangeError, "Requested non-existing UMFPack option");
433
434 UMF_Control[option] = value;
435 }
436
440 void saveDecomposition(const char* file)
441 {
442 int errcode = Caller::save_numeric(UMF_Numeric, const_cast<char*>(file));
443 if (errcode != UMFPACK_OK)
444 DUNE_THROW(Dune::Exception,"IO ERROR while trying to save UMFPack decomposition");
445 }
446
448 void setMatrix(const Matrix& matrix)
449 {
450 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
451 free();
452 if (matrix.N() == 0 or matrix.M() == 0)
453 return;
454
455 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
456 umfpackMatrix_.free();
457 umfpackMatrix_.setSize(MatrixDimension<Matrix>::rowdim(matrix),
458 MatrixDimension<Matrix>::coldim(matrix));
459 ISTL::Impl::BCCSMatrixInitializer<Matrix, size_type> initializer(umfpackMatrix_);
460
461 copyToBCCSMatrix(initializer, matrix);
462
463 decompose();
464 }
465
466 template<class S>
467 void setSubMatrix(const Matrix& _mat, const S& rowIndexSet)
468 {
469 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
470 free();
471
472 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
473 umfpackMatrix_.free();
474
475 umfpackMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(_mat) / _mat.N(),
476 rowIndexSet.size()*MatrixDimension<Matrix>::coldim(_mat) / _mat.M());
477 ISTL::Impl::BCCSMatrixInitializer<Matrix, SuiteSparse_long> initializer(umfpackMatrix_);
478
479 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(_mat,rowIndexSet));
480
481 decompose();
482 }
483
491 void setVerbosity(int v)
492 {
493 verbosity_ = v;
494 // set the verbosity level in UMFPack
495 if (verbosity_ == 0)
496 UMF_Control[UMFPACK_PRL] = 1;
497 if (verbosity_ == 1)
498 UMF_Control[UMFPACK_PRL] = 2;
499 if (verbosity_ == 2)
500 UMF_Control[UMFPACK_PRL] = 4;
501 }
502
508 {
509 return UMF_Numeric;
510 }
511
517 {
518 return umfpackMatrix_;
519 }
520
525 void free()
526 {
527 if (!matrixIsLoaded_)
528 {
529 Caller::free_symbolic(&UMF_Symbolic);
530 umfpackMatrix_.free();
531 }
532 Caller::free_numeric(&UMF_Numeric);
533 matrixIsLoaded_ = false;
534 }
535
536 const char* name() { return "UMFPACK"; }
537
538 private:
539 typedef typename Dune::UMFPackMethodChooser<T> Caller;
540
541 template<class Mat,class X, class TM, class TD, class T1>
542 friend class SeqOverlappingSchwarz;
543 friend struct SeqOverlappingSchwarzAssemblerHelper<UMFPack<Matrix>,true>;
544
546 void decompose()
547 {
548 double UMF_Decomposition_Info[UMFPACK_INFO];
549 Caller::symbolic(static_cast<SuiteSparse_long>(umfpackMatrix_.N()),
550 static_cast<SuiteSparse_long>(umfpackMatrix_.N()),
551 umfpackMatrix_.getColStart(),
552 umfpackMatrix_.getRowIndex(),
553 reinterpret_cast<double*>(umfpackMatrix_.getValues()),
554 &UMF_Symbolic,
555 UMF_Control,
556 UMF_Decomposition_Info);
557 Caller::numeric(umfpackMatrix_.getColStart(),
558 umfpackMatrix_.getRowIndex(),
559 reinterpret_cast<double*>(umfpackMatrix_.getValues()),
560 UMF_Symbolic,
561 &UMF_Numeric,
562 UMF_Control,
563 UMF_Decomposition_Info);
564 Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
565 if (verbosity_ == 1)
566 {
567 std::cout << "[UMFPack Decomposition]" << std::endl;
568 std::cout << "Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME] << " (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] << ")" << std::endl;
569 std::cout << "Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
570 std::cout << "Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT] << " bytes" << std::endl;
571 std::cout << "Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
572 std::cout << "Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ] << " U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
573 }
574 if (verbosity_ == 2)
575 {
576 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
577 }
578 }
579
580 void printOnApply(double* UMF_Info)
581 {
582 Caller::report_status(UMF_Control,UMF_Info[UMFPACK_STATUS]);
583 if (verbosity_ > 0)
584 {
585 std::cout << "[UMFPack Solve]" << std::endl;
586 std::cout << "Wallclock Time: " << UMF_Info[UMFPACK_SOLVE_WALLTIME] << " (CPU Time: " << UMF_Info[UMFPACK_SOLVE_TIME] << ")" << std::endl;
587 std::cout << "Flops Taken: " << UMF_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
588 std::cout << "Iterative Refinement steps taken: " << UMF_Info[UMFPACK_IR_TAKEN] << std::endl;
589 std::cout << "Error Estimate: " << UMF_Info[UMFPACK_OMEGA1] << " resp. " << UMF_Info[UMFPACK_OMEGA2] << std::endl;
590 }
591 }
592
593 UMFPackMatrix umfpackMatrix_;
594 bool matrixIsLoaded_;
595 int verbosity_;
596 void *UMF_Symbolic;
597 void *UMF_Numeric;
598 double UMF_Control[UMFPACK_CONTROL];
599 };
600
601 template<typename T, typename A, int n, int m>
602 struct IsDirectSolver<UMFPack<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
603 {
604 enum { value=true};
605 };
606
607 template<typename T, typename A>
608 struct StoresColumnCompressed<UMFPack<BCRSMatrix<T,A> > >
609 {
610 enum { value = true };
611 };
612
613 struct UMFPackCreator {
614 template<class F,class=void> struct isValidBlock : std::false_type{};
615 template<class B> struct isValidBlock<B, std::enable_if_t<std::is_same<typename FieldTraits<B>::real_type,double>::value>> : std::true_type {};
616
617 template<typename TL, typename M>
618 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
619 typename Dune::TypeListElement<2, TL>::type>>
620 operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
621 std::enable_if_t<
622 isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
623 {
624 int verbose = config.get("verbose", 0);
625 return std::make_shared<Dune::UMFPack<M>>(mat,verbose);
626 }
627
628 // second version with SFINAE to validate the template parameters of UMFPack
629 template<typename TL, typename M>
630 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
631 typename Dune::TypeListElement<2, TL>::type>>
632 operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
633 std::enable_if_t<
634 !isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
635 {
636 DUNE_THROW(UnsupportedType,
637 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
638 }
639 };
640 DUNE_REGISTER_DIRECT_SOLVER("umfpack",Dune::UMFPackCreator());
641} // end namespace Dune
642
643#endif // HAVE_SUITESPARSE_UMFPACK
644
645#endif //DUNE_ISTL_UMFPACK_HH
Implementation of the BCRSMatrix class.
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
derive error class from the base class in common
Definition: istlexception.hh:19
Abstract base class for all solvers.
Definition: solver.hh:99
A generic dynamic dense matrix.
Definition: matrix.hh:561
size_type M() const
Return the number of columns.
Definition: matrix.hh:700
size_type N() const
Return the number of rows.
Definition: matrix.hh:695
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:185
Default exception class for range errors.
Definition: exceptions.hh:254
The UMFPack direct sparse solver.
Definition: umfpack.hh:216
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
void free()
free allocated space.
Definition: umfpack.hh:525
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: umfpack.hh:362
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: umfpack.hh:236
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: umfpack.hh:448
UMFPack(const Matrix &mat_, const ParameterTree &config)
Construct a solver object from a matrix.
Definition: umfpack.hh:286
UMFPack(const Matrix &mat_, const char *file, int verbose=0)
Try loading a decomposition from file and do a decomposition if unsuccessful.
Definition: umfpack.hh:310
typename Impl::UMFPackVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: umfpack.hh:233
UMFPack()
default constructor
Definition: umfpack.hh:292
void apply(T *x, T *b)
additional apply method with c-arrays in analogy to superlu
Definition: umfpack.hh:403
void setVerbosity(int v)
sets the verbosity level for the UMFPack solver
Definition: umfpack.hh:491
UMFPack(const char *file, int verbose=0)
try loading a decomposition from file
Definition: umfpack.hh:337
void saveDecomposition(const char *file)
saves a decomposition to a file
Definition: umfpack.hh:440
UMFPackMatrix & getInternalMatrix()
Return the column compress matrix from UMFPack.
Definition: umfpack.hh:516
typename Impl::UMFPackVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: umfpack.hh:231
UMFPack(const Matrix &matrix, int verbose=0)
Construct a solver object from a matrix.
Definition: umfpack.hh:249
ISTL::Impl::BCCSMatrixInitializer< M, size_type > MatrixInitializer
Type of an associated initializer class.
Definition: umfpack.hh:229
virtual void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: umfpack.hh:393
ISTL::Impl::BCCSMatrix< typename Matrix::field_type, size_type > UMFPackMatrix
The corresponding UMFPack matrix type.
Definition: umfpack.hh:227
void setOption(unsigned int option, double value)
Set UMFPack-specific options.
Definition: umfpack.hh:429
M Matrix
The matrix type.
Definition: umfpack.hh:224
UMFPack(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: umfpack.hh:267
void * getFactorization()
Return the matrix factorization.
Definition: umfpack.hh:507
Dune namespace.
Definition: alignedallocator.hh:13
STL namespace.
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Statistics about the application of an inverse operator.
Definition: solver.hh:48
double elapsed
Elapsed time in seconds.
Definition: solver.hh:82
int iterations
Number of iterations.
Definition: solver.hh:67
bool converged
True if convergence criterion has been met.
Definition: solver.hh:73
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)