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>
37 template<
class M,
class T,
class TM,
class TD,
class TA>
38 class SeqOverlappingSchwarz;
40 template<
class T,
bool tag>
41 struct SeqOverlappingSchwarzAssemblerHelper;
48 template<
class Matrix>
57 static constexpr
bool valid = false ;
63 static constexpr
bool valid = true ;
65 template<
typename... A>
68 umfpack_di_defaults(args...);
70 template<
typename... A>
73 umfpack_di_free_numeric(args...);
75 template<
typename... A>
78 umfpack_di_free_symbolic(args...);
80 template<
typename... A>
83 return umfpack_di_load_numeric(args...);
85 template<
typename... A>
88 umfpack_di_numeric(args...);
90 template<
typename... A>
93 umfpack_di_report_info(args...);
95 template<
typename... A>
98 umfpack_di_report_status(args...);
100 template<
typename... A>
103 return umfpack_di_save_numeric(args...);
105 template<
typename... A>
108 umfpack_di_solve(args...);
110 template<
typename... A>
113 umfpack_di_symbolic(args...);
120 static constexpr
bool valid = true ;
122 template<
typename... A>
125 umfpack_zi_defaults(args...);
127 template<
typename... A>
130 umfpack_zi_free_numeric(args...);
132 template<
typename... A>
135 umfpack_zi_free_symbolic(args...);
137 template<
typename... A>
140 return umfpack_zi_load_numeric(args...);
142 template<
typename... A>
143 static void numeric(
const int* cs,
const int* ri,
const double* val, A... args)
145 umfpack_zi_numeric(cs,ri,val,NULL,args...);
147 template<
typename... A>
150 umfpack_zi_report_info(args...);
152 template<
typename... A>
155 umfpack_zi_report_status(args...);
157 template<
typename... A>
160 return umfpack_zi_save_numeric(args...);
162 template<
typename... A>
163 static void solve(
int m,
const int* cs,
const int* ri, std::complex<double>* val,
double* x,
const double* b,A... args)
165 const double* cval =
reinterpret_cast<const double*
>(val);
166 umfpack_zi_solve(m,cs,ri,cval,NULL,x,NULL,b,NULL,args...);
168 template<
typename... A>
169 static void symbolic(
int m,
int n,
const int* cs,
const int* ri,
const double* val, A... args)
171 umfpack_zi_symbolic(m,n,cs,ri,val,NULL,args...);
188 template<
typename T,
typename A,
int n,
int m>
191 BlockVector<FieldVector<T,m>,
192 typename A::template rebind<FieldVector<T,m> >::other>,
193 BlockVector<FieldVector<T,n>,
194 typename A::template rebind<FieldVector<T,n> >::other> >
207 typename A::template rebind<FieldVector<T,m> >::other>
domain_type;
211 typename A::template rebind<FieldVector<T,n> >::other>
range_type;
216 return SolverCategory::Category::sequential;
230 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
231 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
232 Caller::defaults(UMF_Control);
233 setVerbosity(verbose);
248 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
249 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
250 Caller::defaults(UMF_Control);
251 setVerbosity(verbose);
257 UMFPack() : matrixIsLoaded_(false), verbosity_(0)
260 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
261 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
262 Caller::defaults(UMF_Control);
278 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
279 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
280 Caller::defaults(UMF_Control);
281 setVerbosity(verbose);
282 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
283 if ((errcode == UMFPACK_ERROR_out_of_memory) || (errcode == UMFPACK_ERROR_file_IO))
285 matrixIsLoaded_ =
false;
287 saveDecomposition(file);
291 matrixIsLoaded_ =
true;
292 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
305 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
306 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
307 Caller::defaults(UMF_Control);
308 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
309 if (errcode == UMFPACK_ERROR_out_of_memory)
310 DUNE_THROW(Dune::Exception,
"ran out of memory while loading UMFPack decomposition");
311 if (errcode == UMFPACK_ERROR_file_IO)
312 DUNE_THROW(Dune::Exception,
"IO error while loading UMFPack decomposition");
313 matrixIsLoaded_ =
true;
314 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
315 setVerbosity(verbose);
320 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
329 if (umfpackMatrix_.N() != b.dim())
330 DUNE_THROW(
Dune::ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
331 if (umfpackMatrix_.M() != x.dim())
332 DUNE_THROW(
Dune::ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
334 double UMF_Apply_Info[UMFPACK_INFO];
335 Caller::solve(UMFPACK_A,
336 umfpackMatrix_.getColStart(),
337 umfpackMatrix_.getRowIndex(),
338 umfpackMatrix_.getValues(),
339 reinterpret_cast<double*
>(&x[0]),
340 reinterpret_cast<double*
>(&b[0]),
348 res.
elapsed = UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME];
350 printOnApply(UMF_Apply_Info);
358 DUNE_UNUSED_PARAMETER(reduction);
369 double UMF_Apply_Info[UMFPACK_INFO];
370 Caller::solve(UMFPACK_A,
371 umfpackMatrix_.getColStart(),
372 umfpackMatrix_.getRowIndex(),
373 umfpackMatrix_.getValues(),
379 printOnApply(UMF_Apply_Info);
395 if (option >= UMFPACK_CONTROL)
396 DUNE_THROW(RangeError,
"Requested non-existing UMFPack option");
398 UMF_Control[option] = value;
406 int errcode = Caller::save_numeric(UMF_Numeric,
const_cast<char*
>(file));
407 if (errcode != UMFPACK_OK)
408 DUNE_THROW(Dune::Exception,
"IO ERROR while trying to save UMFPack decomposition");
414 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
416 umfpackMatrix_ = matrix;
423 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
425 umfpackMatrix_.setMatrix(_mat,rowIndexSet);
441 UMF_Control[UMFPACK_PRL] = 1;
443 UMF_Control[UMFPACK_PRL] = 2;
445 UMF_Control[UMFPACK_PRL] = 4;
463 return umfpackMatrix_;
472 if (!matrixIsLoaded_)
474 Caller::free_symbolic(&UMF_Symbolic);
475 umfpackMatrix_.free();
477 Caller::free_numeric(&UMF_Numeric);
478 matrixIsLoaded_ =
false;
481 const char*
name() {
return "UMFPACK"; }
486 template<
class M,
class X,
class TM,
class TD,
class T1>
493 double UMF_Decomposition_Info[UMFPACK_INFO];
494 Caller::symbolic(
static_cast<int>(umfpackMatrix_.N()),
495 static_cast<int>(umfpackMatrix_.N()),
496 umfpackMatrix_.getColStart(),
497 umfpackMatrix_.getRowIndex(),
498 reinterpret_cast<double*
>(umfpackMatrix_.getValues()),
501 UMF_Decomposition_Info);
502 Caller::numeric(umfpackMatrix_.getColStart(),
503 umfpackMatrix_.getRowIndex(),
504 reinterpret_cast<double*
>(umfpackMatrix_.getValues()),
508 UMF_Decomposition_Info);
509 Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
512 std::cout <<
"[UMFPack Decomposition]" << std::endl;
513 std::cout <<
"Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME] <<
" (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] <<
")" << std::endl;
514 std::cout <<
"Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
515 std::cout <<
"Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT] <<
" bytes" << std::endl;
516 std::cout <<
"Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
517 std::cout <<
"Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ] <<
" U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
521 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
525 void printOnApply(
double* UMF_Info)
527 Caller::report_status(UMF_Control,UMF_Info[UMFPACK_STATUS]);
530 std::cout <<
"[UMFPack Solve]" << std::endl;
531 std::cout <<
"Wallclock Time: " << UMF_Info[UMFPACK_SOLVE_WALLTIME] <<
" (CPU Time: " << UMF_Info[UMFPACK_SOLVE_TIME] <<
")" << std::endl;
532 std::cout <<
"Flops Taken: " << UMF_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
533 std::cout <<
"Iterative Refinement steps taken: " << UMF_Info[UMFPACK_IR_TAKEN] << std::endl;
534 std::cout <<
"Error Estimate: " << UMF_Info[UMFPACK_OMEGA1] <<
" resp. " << UMF_Info[UMFPACK_OMEGA2] << std::endl;
538 UMFPackMatrix umfpackMatrix_;
539 bool matrixIsLoaded_;
543 double UMF_Control[UMFPACK_CONTROL];
546 template<
typename T,
typename A,
int n,
int m>
552 template<
typename T,
typename A,
int n,
int m>
559 #endif // HAVE_SUITESPARSE_UMFPACK
561 #endif //DUNE_ISTL_UMFPACK_HH