Go to the documentation of this file.
3 #ifndef DUNE_ISTL_UMFPACK_HH
4 #define DUNE_ISTL_UMFPACK_HH
6 #if HAVE_SUITESPARSE_UMFPACK || defined DOXYGEN
13 #include<dune/common/exceptions.hh>
14 #include<dune/common/fmatrix.hh>
15 #include<dune/common/fvector.hh>
16 #include<dune/common/unused.hh>
38 template<
class M,
class T,
class TM,
class TD,
class TA>
39 class SeqOverlappingSchwarz;
41 template<
class T,
bool tag>
42 struct SeqOverlappingSchwarzAssemblerHelper;
49 static constexpr
bool valid = false ;
55 static constexpr
bool valid = true ;
57 template<
typename... A>
60 umfpack_dl_defaults(args...);
62 template<
typename... A>
65 umfpack_dl_free_numeric(args...);
67 template<
typename... A>
70 umfpack_dl_free_symbolic(args...);
72 template<
typename... A>
75 return umfpack_dl_load_numeric(args...);
77 template<
typename... A>
80 umfpack_dl_numeric(args...);
82 template<
typename... A>
85 umfpack_dl_report_info(args...);
87 template<
typename... A>
90 umfpack_dl_report_status(args...);
92 template<
typename... A>
95 return umfpack_dl_save_numeric(args...);
97 template<
typename... A>
100 umfpack_dl_solve(args...);
102 template<
typename... A>
105 umfpack_dl_symbolic(args...);
112 static constexpr
bool valid = true ;
114 template<
typename... A>
117 umfpack_zl_defaults(args...);
119 template<
typename... A>
122 umfpack_zl_free_numeric(args...);
124 template<
typename... A>
127 umfpack_zl_free_symbolic(args...);
129 template<
typename... A>
132 return umfpack_zl_load_numeric(args...);
134 template<
typename... A>
135 static void numeric(
const long int* cs,
const long int* ri,
const double* val, A... args)
137 umfpack_zl_numeric(cs,ri,val,NULL,args...);
139 template<
typename... A>
142 umfpack_zl_report_info(args...);
144 template<
typename... A>
147 umfpack_zl_report_status(args...);
149 template<
typename... A>
152 return umfpack_zl_save_numeric(args...);
154 template<
typename... A>
155 static void solve(
long int m,
const long int* cs,
const long int* ri, std::complex<double>* val,
double* x,
const double* b,A... args)
157 const double* cval =
reinterpret_cast<const double*
>(val);
158 umfpack_zl_solve(m,cs,ri,cval,NULL,x,NULL,b,NULL,args...);
160 template<
typename... A>
161 static void symbolic(
long int m,
long int n,
const long int* cs,
const long int* ri,
const double* val, A... args)
163 umfpack_zl_symbolic(m,n,cs,ri,val,NULL,args...);
170 struct UMFPackVectorChooser
173 template<
typename T,
typename A,
int n,
int m>
174 struct UMFPackVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
177 using domain_type = BlockVector<
179 typename A::template rebind<FieldVector<T,m> >::other>;
181 using range_type = BlockVector<
183 typename A::template rebind<FieldVector<T,n> >::other>;
186 template<
typename T,
typename A>
187 struct UMFPackVectorChooser<BCRSMatrix<T,A> >
190 using domain_type = BlockVector<T, A>;
192 using range_type = BlockVector<T, A>;
212 typename Impl::UMFPackVectorChooser<M>::domain_type,
213 typename Impl::UMFPackVectorChooser<M>::range_type >
215 using T =
typename M::field_type;
226 using domain_type =
typename Impl::UMFPackVectorChooser<M>::domain_type;
228 using range_type =
typename Impl::UMFPackVectorChooser<M>::range_type;
233 return SolverCategory::Category::sequential;
247 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
248 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
249 Caller::defaults(UMF_Control);
265 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
266 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
267 Caller::defaults(UMF_Control);
274 UMFPack() : matrixIsLoaded_(false), verbosity_(0)
277 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
278 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
279 Caller::defaults(UMF_Control);
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);
299 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
300 if ((errcode == UMFPACK_ERROR_out_of_memory) || (errcode == UMFPACK_ERROR_file_IO))
302 matrixIsLoaded_ =
false;
308 matrixIsLoaded_ =
true;
309 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
322 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
323 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
324 Caller::defaults(UMF_Control);
325 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
326 if (errcode == UMFPACK_ERROR_out_of_memory)
327 DUNE_THROW(Dune::Exception,
"ran out of memory while loading UMFPack decomposition");
328 if (errcode == UMFPACK_ERROR_file_IO)
329 DUNE_THROW(Dune::Exception,
"IO error while loading UMFPack decomposition");
330 matrixIsLoaded_ =
true;
331 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
337 if ((umfpackMatrix_.
N() + umfpackMatrix_.
M() > 0) || matrixIsLoaded_)
346 if (umfpackMatrix_.
N() != b.dim())
347 DUNE_THROW(
Dune::ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
348 if (umfpackMatrix_.
M() != x.dim())
349 DUNE_THROW(
Dune::ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
353 double UMF_Apply_Info[UMFPACK_INFO];
354 Caller::solve(UMFPACK_A,
358 reinterpret_cast<double*
>(&x[0]),
359 reinterpret_cast<double*
>(&b[0]),
367 res.
elapsed = UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME];
369 printOnApply(UMF_Apply_Info);
377 DUNE_UNUSED_PARAMETER(reduction);
388 double UMF_Apply_Info[UMFPACK_INFO];
389 Caller::solve(UMFPACK_A,
398 printOnApply(UMF_Apply_Info);
414 if (option >= UMFPACK_CONTROL)
415 DUNE_THROW(RangeError,
"Requested non-existing UMFPack option");
417 UMF_Control[option] = value;
425 int errcode = Caller::save_numeric(UMF_Numeric,
const_cast<char*
>(file));
426 if (errcode != UMFPACK_OK)
427 DUNE_THROW(Dune::Exception,
"IO ERROR while trying to save UMFPack decomposition");
433 if ((umfpackMatrix_.
N() + umfpackMatrix_.
M() > 0) || matrixIsLoaded_)
435 if (matrix.N() == 0 or matrix.M() == 0)
437 umfpackMatrix_ = matrix;
444 if ((umfpackMatrix_.
N() + umfpackMatrix_.
M() > 0) || matrixIsLoaded_)
446 umfpackMatrix_.
setMatrix(_mat,rowIndexSet);
462 UMF_Control[UMFPACK_PRL] = 1;
464 UMF_Control[UMFPACK_PRL] = 2;
466 UMF_Control[UMFPACK_PRL] = 4;
484 return umfpackMatrix_;
493 if (!matrixIsLoaded_)
495 Caller::free_symbolic(&UMF_Symbolic);
496 umfpackMatrix_.
free();
498 Caller::free_numeric(&UMF_Numeric);
499 matrixIsLoaded_ =
false;
502 const char*
name() {
return "UMFPACK"; }
507 template<
class Mat,
class X,
class TM,
class TD,
class T1>
514 double UMF_Decomposition_Info[UMFPACK_INFO];
515 Caller::symbolic(
static_cast<int>(umfpackMatrix_.
N()),
516 static_cast<int>(umfpackMatrix_.
N()),
519 reinterpret_cast<double*
>(umfpackMatrix_.
getValues()),
522 UMF_Decomposition_Info);
525 reinterpret_cast<double*
>(umfpackMatrix_.
getValues()),
529 UMF_Decomposition_Info);
530 Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
533 std::cout <<
"[UMFPack Decomposition]" << std::endl;
534 std::cout <<
"Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME] <<
" (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] <<
")" << std::endl;
535 std::cout <<
"Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
536 std::cout <<
"Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT] <<
" bytes" << std::endl;
537 std::cout <<
"Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
538 std::cout <<
"Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ] <<
" U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
542 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
546 void printOnApply(
double* UMF_Info)
548 Caller::report_status(UMF_Control,UMF_Info[UMFPACK_STATUS]);
551 std::cout <<
"[UMFPack Solve]" << std::endl;
552 std::cout <<
"Wallclock Time: " << UMF_Info[UMFPACK_SOLVE_WALLTIME] <<
" (CPU Time: " << UMF_Info[UMFPACK_SOLVE_TIME] <<
")" << std::endl;
553 std::cout <<
"Flops Taken: " << UMF_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
554 std::cout <<
"Iterative Refinement steps taken: " << UMF_Info[UMFPACK_IR_TAKEN] << std::endl;
555 std::cout <<
"Error Estimate: " << UMF_Info[UMFPACK_OMEGA1] <<
" resp. " << UMF_Info[UMFPACK_OMEGA2] << std::endl;
560 bool matrixIsLoaded_;
564 double UMF_Control[UMFPACK_CONTROL];
567 template<
typename T,
typename A,
int n,
int m>
573 template<
typename T,
typename A>
581 template<
class B>
struct isValidBlock<B, std::enable_if_t<std::is_same<typename FieldTraits<B>::real_type,double>::value>> : std::true_type {};
583 template<
typename TL,
typename M>
584 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
585 typename Dune::TypeListElement<2, TL>::type>>
588 isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
590 int verbose = config.get(
"verbose", 0);
591 return std::make_shared<Dune::UMFPack<M>>(
mat,verbose);
595 template<
typename TL,
typename M>
596 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
597 typename Dune::TypeListElement<2, TL>::type>>
600 !
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
603 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
609 #endif // HAVE_SUITESPARSE_UMFPACK
611 #endif //DUNE_ISTL_UMFPACK_HH
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition: umfpack.hh:431
virtual void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: umfpack.hh:375
Index * getRowIndex() const
Definition: colcompmatrix.hh:221
void setVerbosity(int v)
sets the verbosity level for the UMFPack solver
Definition: umfpack.hh:457
std::shared_ptr< Dune::InverseOperator< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL, const M &mat, const Dune::ParameterTree &config, std::enable_if_t< isValidBlock< typename Dune::TypeListElement< 1, TL >::type::block_type >::value, int >=0) const
Definition: umfpack.hh:586
M Matrix
The matrix type.
Definition: umfpack.hh:219
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:292
virtual void free()
free allocated space.
static int load_numeric(A... args)
Definition: umfpack.hh:73
double elapsed
Elapsed time in seconds.
Definition: solver.hh:80
static void free_numeric(A... args)
Definition: umfpack.hh:63
DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator())
static void free_symbolic(A... args)
Definition: umfpack.hh:68
void apply(T *x, T *b)
additional apply method with c-arrays in analogy to superlu
Definition: umfpack.hh:386
typename Impl::UMFPackVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: umfpack.hh:226
static void symbolic(long int m, long int n, const long int *cs, const long int *ri, const double *val, A... args)
Definition: umfpack.hh:161
static int save_numeric(A... args)
Definition: umfpack.hh:150
Sequential overlapping Schwarz preconditioner.
Definition: colcompmatrix.hh:150
Definition: solverfactory.hh:124
void setSubMatrix(const Matrix &_mat, const S &rowIndexSet)
Definition: umfpack.hh:442
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: umfpack.hh:344
Definition: matrixutils.hh:26
Abstract base class for all solvers.
Definition: solver.hh:97
static void report_info(A... args)
Definition: umfpack.hh:140
static void report_status(A... args)
Definition: umfpack.hh:88
UMFPack(const Matrix &matrix, int verbose=0)
Construct a solver object from a matrix.
Definition: umfpack.hh:244
typename Impl::UMFPackVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: umfpack.hh:228
Definition: umfpack.hh:580
static constexpr bool valid
Definition: umfpack.hh:49
void setOption(unsigned int option, double value)
Set UMFPack-specific options.
Definition: umfpack.hh:412
@ value
Whether this is a direct solver.
Definition: solvertype.hh:22
virtual ~UMFPack()
Definition: umfpack.hh:335
Dune::ColCompMatrix< Matrix, long int > UMFPackMatrix
The corresponding UMFPack matrix type.
Definition: umfpack.hh:222
const char * name()
Definition: umfpack.hh:502
Index * getColStart() const
Definition: colcompmatrix.hh:226
static void free_numeric(A... args)
Definition: umfpack.hh:120
void saveDecomposition(const char *file)
saves a decomposition to a file
Definition: umfpack.hh:423
ColCompMatrixInitializer< M, long int > MatrixInitializer
Type of an associated initializer class.
Definition: umfpack.hh:224
UMFPack()
default constructor
Definition: umfpack.hh:274
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
static void numeric(A... args)
Definition: umfpack.hh:78
Definition: solvertype.hh:27
Statistics about the application of an inverse operator.
Definition: solver.hh:45
void free()
free allocated space.
Definition: umfpack.hh:491
static void solve(long int m, const long int *cs, const long int *ri, std::complex< double > *val, double *x, const double *b, A... args)
Definition: umfpack.hh:155
int iterations
Number of iterations.
Definition: solver.hh:65
static void report_info(A... args)
Definition: umfpack.hh:83
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: umfpack.hh:231
Implementation of the BCRSMatrix class.
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:424
static void defaults(A... args)
Definition: umfpack.hh:58
derive error class from the base class in common
Definition: istlexception.hh:16
virtual void setMatrix(const Matrix &mat, const std::set< std::size_t > &mrs)
Initialize data from a given set of matrix rows and columns.
UMFPack(const char *file, int verbose=0)
try loading a decomposition from file
Definition: umfpack.hh:319
Definition: umfpack.hh:579
@ value
whether the solver internally uses column compressed storage
Definition: solvertype.hh:34
Definition: colcompmatrix.hh:153
Implementations of the inverse operator interface.
Definition: allocator.hh:7
Definition: solvertype.hh:13
void * getFactorization()
Return the matrix factorization.
Definition: umfpack.hh:473
static void solve(A... args)
Definition: umfpack.hh:98
Inititializer for the ColCompMatrix as needed by OverlappingSchwarz.
Definition: colcompmatrix.hh:147
static int save_numeric(A... args)
Definition: umfpack.hh:93
size_type N() const
Get the number of rows.
Definition: colcompmatrix.hh:192
Matrix & mat
Definition: matrixmatrix.hh:345
static void report_status(A... args)
Definition: umfpack.hh:145
UMFPack(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: umfpack.hh:262
Definition: umfpack.hh:47
The UMFPack direct sparse solver.
Definition: umfpack.hh:210
Templates characterizing the type of a solver.
Category
Definition: solvercategory.hh:21
static int load_numeric(A... args)
Definition: umfpack.hh:130
UMFPackMatrix & getInternalMatrix()
Return the column compress matrix from UMFPack.
Definition: umfpack.hh:482
static void symbolic(A... args)
Definition: umfpack.hh:103
M matrix_type
Definition: umfpack.hh:220
static void defaults(A... args)
Definition: umfpack.hh:115
static void free_symbolic(A... args)
Definition: umfpack.hh:125
B * getValues() const
Definition: colcompmatrix.hh:216
static void numeric(const long int *cs, const long int *ri, const double *val, A... args)
Definition: umfpack.hh:135
size_type M() const
Get the number of columns.
Definition: colcompmatrix.hh:211