│ │ │
│ │ │ -
│ │ │ -
│ │ │ -
5#ifndef DUNE_ISTL_UMFPACK_HH
│ │ │ -
6#define DUNE_ISTL_UMFPACK_HH
│ │ │ -
│ │ │ -
8#if HAVE_SUITESPARSE_UMFPACK || defined DOXYGEN
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
15#include<dune/common/exceptions.hh>
│ │ │ -
16#include<dune/common/fmatrix.hh>
│ │ │ -
17#include<dune/common/fvector.hh>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
43 template<
class M,
class T,
class TM,
class TD,
class TA>
│ │ │ -
│ │ │ -
│ │ │ -
46 template<
class T,
bool tag>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
54 static constexpr bool valid = false ;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
60 static constexpr bool valid = true ;
│ │ │ -
│ │ │ -
62 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
65 umfpack_dl_defaults(args...);
│ │ │ -
│ │ │ -
│ │ │ -
67 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
70 umfpack_dl_free_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
72 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
75 umfpack_dl_free_symbolic(args...);
│ │ │ +
3#ifndef DUNE_ISTL_COMMON_REGISTRY_HH
│ │ │ +
4#define DUNE_ISTL_COMMON_REGISTRY_HH
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
14#include <dune/common/typelist.hh>
│ │ │ +
15#include <dune/common/hybridutilities.hh>
│ │ │ +
16#include <dune/common/parameterizedobject.hh>
│ │ │ +
│ │ │ +
│ │ │ +
18#define DUNE_REGISTRY_PUT(Tag, id, ...) \
│ │ │ +
│ │ │ +
│ │ │ +
21 struct Registry<Tag, DUNE_GET_COUNTER(Tag)> \
│ │ │ +
│ │ │ +
23 static auto getCreator() \
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
27 static std::string name() { return id; } \
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
35 template<
class Tag, std::
size_t index>
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
40 template<
template<
class>
class Base,
class V,
class Tag,
typename... Args>
│ │ │ +
41 auto registryGet(Tag , std::string name, Args... args)
│ │ │ +
│ │ │ +
│ │ │ +
44 std::shared_ptr<Base<V> > result;
│ │ │ +
45 Dune::Hybrid::forEach(std::make_index_sequence<count>{},
│ │ │ +
│ │ │ +
47 using Reg = Registry<Tag, index>;
│ │ │ +
48 if(!result && Reg::name() == name) {
│ │ │ +
49 result = Reg::getCreator()(Dune::MetaType<V>{}, args...);
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
60 template<
class V,
class Type,
class Tag,
class... Args>
│ │ │ +
61 int addRegistryToFactory(Dune::ParameterizedObjectFactory<Type(Args...), std::string>& factory,
│ │ │ +
│ │ │ +
│ │ │ +
64 Dune::Hybrid::forEach(std::make_index_sequence<count>{},
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
69 using Reg = Registry<Tag, index>;
│ │ │ +
70 auto genericcreator = Reg::getCreator();
│ │ │ +
71 factory.define(Reg::name(), [genericcreator](Args... args){
│ │ │ +
72 return genericcreator(V{}, args...);
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │
│ │ │ -
│ │ │ -
77 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
80 return umfpack_dl_load_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
82 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
85 umfpack_dl_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
87 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
90 umfpack_dl_report_info(args...);
│ │ │ -
│ │ │ -
│ │ │ -
92 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
95 umfpack_dl_report_status(args...);
│ │ │ -
│ │ │ -
│ │ │ -
97 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
100 return umfpack_dl_save_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
102 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
105 umfpack_dl_solve(args...);
│ │ │ -
│ │ │ -
│ │ │ -
107 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
110 umfpack_dl_symbolic(args...);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
117 static constexpr bool valid = true ;
│ │ │ -
│ │ │ -
│ │ │ -
120 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
123 umfpack_zl_defaults(args...);
│ │ │ -
│ │ │ -
│ │ │ -
125 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
128 umfpack_zl_free_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
130 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
133 umfpack_zl_free_symbolic(args...);
│ │ │ -
│ │ │ -
│ │ │ -
135 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
138 return umfpack_zl_load_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
140 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
143 umfpack_zl_numeric(cs,ri,val,NULL,args...);
│ │ │ -
│ │ │ -
│ │ │ -
145 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
148 umfpack_zl_report_info(args...);
│ │ │ -
│ │ │ -
│ │ │ -
150 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
153 umfpack_zl_report_status(args...);
│ │ │ -
│ │ │ -
│ │ │ -
155 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
158 return umfpack_zl_save_numeric(args...);
│ │ │ -
│ │ │ -
│ │ │ -
160 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
163 const double* cval =
reinterpret_cast<const double*
>(val);
│ │ │ -
164 umfpack_zl_solve(m,cs,ri,cval,NULL,x,NULL,b,NULL,args...);
│ │ │ -
│ │ │ -
│ │ │ -
166 template<
typename... A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
169 umfpack_zl_symbolic(m,n,cs,ri,val,NULL,args...);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
175 template<
class M,
class =
void>
│ │ │ -
176 struct UMFPackVectorChooser;
│ │ │ -
│ │ │ -
179 template<
class M>
using UMFPackDomainType =
typename UMFPackVectorChooser<M>::domain_type;
│ │ │ -
│ │ │ -
182 template<
class M>
using UMFPackRangeType =
typename UMFPackVectorChooser<M>::range_type;
│ │ │ -
│ │ │ -
│ │ │ -
185 struct UMFPackVectorChooser<M,
│ │ │ -
186 std::enable_if_t<(std::is_same<M,double>::value) || (std::is_same<M,std::complex<double> >::value)>>
│ │ │ -
│ │ │ -
188 using domain_type = M;
│ │ │ -
189 using range_type = M;
│ │ │ -
│ │ │ -
│ │ │ -
192 template<
typename T,
int n,
int m>
│ │ │ -
193 struct UMFPackVectorChooser<FieldMatrix<T,n,m>,
│ │ │ -
194 std::enable_if_t<(std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value)>>
│ │ │ -
│ │ │ -
197 using domain_type = FieldVector<T,m>;
│ │ │ -
199 using range_type = FieldVector<T,n>;
│ │ │ -
│ │ │ -
│ │ │ -
202 template<
typename T,
typename A>
│ │ │ -
203 struct UMFPackVectorChooser<BCRSMatrix<T,A>,
│ │ │ -
204 std::void_t<UMFPackDomainType<T>, UMFPackRangeType<T>>>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
210 using domain_type = BlockVector<UMFPackDomainType<T>,
typename std::allocator_traits<A>::template rebind_alloc<UMFPackDomainType<T>>>;
│ │ │ -
212 using range_type = BlockVector<UMFPackRangeType<T>,
typename std::allocator_traits<A>::template rebind_alloc<UMFPackRangeType<T>>>;
│ │ │ -
│ │ │ -
│ │ │ -
215 template<
typename T,
typename A>
│ │ │ -
216 struct UMFPackVectorChooser<Matrix<T,A>,
│ │ │ -
217 std::void_t<UMFPackDomainType<T>, UMFPackRangeType<T>>>
│ │ │ -
218 :
public UMFPackVectorChooser<BCRSMatrix<T,A>, std::void_t<UMFPackDomainType<T>, UMFPackRangeType<T>>>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
222 template<
typename FirstBlock,
typename... Blocks>
│ │ │ -
223 struct UMFPackVectorChooser<MultiTypeBlockVector<FirstBlock, Blocks...>,
│ │ │ -
224 std::void_t<UMFPackDomainType<FirstBlock>, UMFPackRangeType<FirstBlock>, UMFPackDomainType<Blocks>...>>
│ │ │ -
│ │ │ -
227 using domain_type = MultiTypeBlockVector<UMFPackDomainType<FirstBlock>, UMFPackDomainType<Blocks>...>;
│ │ │ -
229 using range_type = UMFPackRangeType<FirstBlock>;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
233 template<
typename FirstRow,
typename... Rows>
│ │ │ -
234 struct UMFPackVectorChooser<MultiTypeBlockMatrix<FirstRow, Rows...>,
│ │ │ -
235 std::void_t<UMFPackDomainType<FirstRow>, UMFPackRangeType<FirstRow>, UMFPackRangeType<Rows>...>>
│ │ │ -
│ │ │ -
238 using domain_type = UMFPackDomainType<FirstRow>;
│ │ │ -
240 using range_type = MultiTypeBlockVector< UMFPackRangeType<FirstRow>, UMFPackRangeType<Rows>... >;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
266 using T =
typename M::field_type;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
275 using UMFPackMatrix = ISTL::Impl::BCCSMatrix<typename Matrix::field_type, size_type>;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
300 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
│ │ │ -
301 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
302 Caller::defaults(UMF_Control);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
318 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
│ │ │ -
319 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
320 Caller::defaults(UMF_Control);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
335 :
UMFPack(mat_, config.
get<int>(
"verbose", 0))
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
340 UMFPack() : matrixIsLoaded_(false), verbosity_(0)
│ │ │ -
│ │ │ -
│ │ │ -
343 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
│ │ │ -
344 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
345 Caller::defaults(UMF_Control);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
361 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
│ │ │ -
362 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
363 Caller::defaults(UMF_Control);
│ │ │ -
│ │ │ -
365 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
│ │ │ -
366 if ((errcode == UMFPACK_ERROR_out_of_memory) || (errcode == UMFPACK_ERROR_file_IO))
│ │ │ -
│ │ │ -
368 matrixIsLoaded_ =
false;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
374 matrixIsLoaded_ =
true;
│ │ │ -
375 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
388 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
│ │ │ -
389 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
390 Caller::defaults(UMF_Control);
│ │ │ -
391 int errcode = Caller::load_numeric(&UMF_Numeric,
const_cast<char*
>(file));
│ │ │ -
392 if (errcode == UMFPACK_ERROR_out_of_memory)
│ │ │ -
393 DUNE_THROW(Dune::Exception,
"ran out of memory while loading UMFPack decomposition");
│ │ │ -
394 if (errcode == UMFPACK_ERROR_file_IO)
│ │ │ -
395 DUNE_THROW(Dune::Exception,
"IO error while loading UMFPack decomposition");
│ │ │ -
396 matrixIsLoaded_ =
true;
│ │ │ -
397 std::cout <<
"UMFPack decomposition successfully loaded from " << file << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
403 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
412 if (umfpackMatrix_.N() != b.dim())
│ │ │ -
413 DUNE_THROW(
Dune::ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
│ │ │ -
414 if (umfpackMatrix_.M() != x.dim())
│ │ │ -
415 DUNE_THROW(
Dune::ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
421 std::vector<T> xFlat(x.dim()), bFlat(b.dim());
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
425 xFlat[ offset ] = entry;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
430 bFlat[ offset ] = entry;
│ │ │ -
│ │ │ -
│ │ │ -
433 double UMF_Apply_Info[UMFPACK_INFO];
│ │ │ -
434 Caller::solve(UMFPACK_A,
│ │ │ -
435 umfpackMatrix_.getColStart(),
│ │ │ -
436 umfpackMatrix_.getRowIndex(),
│ │ │ -
437 umfpackMatrix_.getValues(),
│ │ │ -
438 reinterpret_cast<double*
>(&xFlat[0]),
│ │ │ -
439 reinterpret_cast<double*
>(&bFlat[0]),
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
447 entry = xFlat[offset];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
453 res.
elapsed = UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME];
│ │ │ -
│ │ │ -
455 printOnApply(UMF_Apply_Info);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
475 double UMF_Apply_Info[UMFPACK_INFO];
│ │ │ -
476 Caller::solve(UMFPACK_A,
│ │ │ -
477 umfpackMatrix_.getColStart(),
│ │ │ -
478 umfpackMatrix_.getRowIndex(),
│ │ │ -
479 umfpackMatrix_.getValues(),
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
485 printOnApply(UMF_Apply_Info);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
501 if (option >= UMFPACK_CONTROL)
│ │ │ -
502 DUNE_THROW(RangeError,
"Requested non-existing UMFPack option");
│ │ │ -
│ │ │ -
504 UMF_Control[option] = value;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
512 int errcode = Caller::save_numeric(UMF_Numeric,
const_cast<char*
>(file));
│ │ │ -
513 if (errcode != UMFPACK_OK)
│ │ │ -
514 DUNE_THROW(Dune::Exception,
"IO ERROR while trying to save UMFPack decomposition");
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
526 template<
class BitVector = Impl::NoBitVector>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
529 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
│ │ │ -
│ │ │ -
531 if (matrix.N() == 0 or matrix.M() == 0)
│ │ │ -
│ │ │ -
│ │ │ -
534 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
│ │ │ -
535 umfpackMatrix_.free();
│ │ │ -
│ │ │ -
537 constexpr bool useBitVector = not std::is_same_v<BitVector,Impl::NoBitVector>;
│ │ │ -
│ │ │ -
│ │ │ -
540 std::vector<bool> flatBitVector;
│ │ │ -
│ │ │ -
542 std::vector<size_type> subIndices;
│ │ │ -
│ │ │ -
544 [[maybe_unused]]
int numberOfIgnoredDofs = 0;
│ │ │ -
│ │ │ -
│ │ │ -
547 if constexpr ( useBitVector )
│ │ │ -
│ │ │ -
│ │ │ -
550 flatBitVector.resize(flatSize);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
554 flatBitVector[ offset ] = entry;
│ │ │ -
│ │ │ -
│ │ │ -
557 numberOfIgnoredDofs++;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
565 if constexpr ( useBitVector )
│ │ │ -
566 if ( flatBitVector[row] or flatBitVector[
col] )
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
572 if constexpr ( useBitVector )
│ │ │ -
│ │ │ -
│ │ │ -
575 subIndices.resize(flatRows,std::numeric_limits<std::size_t>::max());
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
579 if ( not flatBitVector[ i ] )
│ │ │ -
580 subIndices[ i ] = subIndexCounter++;
│ │ │ -
│ │ │ -
│ │ │ -
583 flatRows -= numberOfIgnoredDofs;
│ │ │ -
584 flatCols -= numberOfIgnoredDofs;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
588 umfpackMatrix_.setSize(flatRows,flatCols);
│ │ │ -
589 umfpackMatrix_.Nnz_ = nonZeros;
│ │ │ -
│ │ │ -
│ │ │ -
592 umfpackMatrix_.colstart =
new size_type[flatCols+1];
│ │ │ -
593 umfpackMatrix_.rowindex =
new size_type[nonZeros];
│ │ │ -
594 umfpackMatrix_.values =
new T[nonZeros];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
598 umfpackMatrix_.colstart[i] = 0;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
606 if constexpr ( useBitVector )
│ │ │ -
607 if ( flatBitVector[flatRowIndex] or flatBitVector[flatColIndex] )
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
612 auto colIdx = useBitVector ? subIndices[flatColIndex] : flatColIndex;
│ │ │ -
│ │ │ -
614 umfpackMatrix_.colstart[colIdx+1]++;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
620 umfpackMatrix_.colstart[i+1] += umfpackMatrix_.colstart[i];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
624 std::vector<size_type> colPosition(flatCols,0);
│ │ │ -
│ │ │ -
│ │ │ -
627 flatMatrixForEach(matrix, [&](
auto&& entry,
auto&& flatRowIndex,
auto&& flatColIndex)
│ │ │ -
│ │ │ -
│ │ │ -
630 if constexpr ( useBitVector )
│ │ │ -
631 if ( flatBitVector[flatRowIndex] or flatBitVector[flatColIndex] )
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
636 auto rowIdx = useBitVector ? subIndices[flatRowIndex] : flatRowIndex;
│ │ │ -
637 auto colIdx = useBitVector ? subIndices[flatColIndex] : flatColIndex;
│ │ │ -
│ │ │ -
│ │ │ -
640 auto colStart = umfpackMatrix_.colstart[colIdx];
│ │ │ -
│ │ │ -
642 auto colPos = colPosition[colIdx];
│ │ │ -
│ │ │ -
644 umfpackMatrix_.rowindex[ colStart + colPos ] = rowIdx;
│ │ │ -
645 umfpackMatrix_.values[ colStart + colPos ] = entry;
│ │ │ -
│ │ │ -
647 colPosition[colIdx]++;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
659 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
│ │ │ -
│ │ │ -
│ │ │ -
662 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
│ │ │ -
663 umfpackMatrix_.free();
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
667 ISTL::Impl::BCCSMatrixInitializer<Matrix, SuiteSparse_long> initializer(umfpackMatrix_);
│ │ │ -
│ │ │ -
669 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<
Matrix,std::set<typename Matrix::size_type> >(_mat,rowIndexSet));
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
686 UMF_Control[UMFPACK_PRL] = 1;
│ │ │ -
│ │ │ -
688 UMF_Control[UMFPACK_PRL] = 2;
│ │ │ -
│ │ │ -
690 UMF_Control[UMFPACK_PRL] = 4;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
708 return umfpackMatrix_;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
717 if (!matrixIsLoaded_)
│ │ │ -
│ │ │ -
719 Caller::free_symbolic(&UMF_Symbolic);
│ │ │ -
720 umfpackMatrix_.free();
│ │ │ -
│ │ │ -
722 Caller::free_numeric(&UMF_Numeric);
│ │ │ -
723 matrixIsLoaded_ =
false;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
726 const char*
name() {
return "UMFPACK"; }
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
731 template<
class Mat,
class X,
class TM,
class TD,
class T1>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
738 double UMF_Decomposition_Info[UMFPACK_INFO];
│ │ │ -
739 Caller::symbolic(
static_cast<SuiteSparse_long
>(umfpackMatrix_.N()),
│ │ │ -
740 static_cast<SuiteSparse_long
>(umfpackMatrix_.N()),
│ │ │ -
741 umfpackMatrix_.getColStart(),
│ │ │ -
742 umfpackMatrix_.getRowIndex(),
│ │ │ -
743 reinterpret_cast<double*
>(umfpackMatrix_.getValues()),
│ │ │ -
│ │ │ -
│ │ │ -
746 UMF_Decomposition_Info);
│ │ │ -
747 Caller::numeric(umfpackMatrix_.getColStart(),
│ │ │ -
748 umfpackMatrix_.getRowIndex(),
│ │ │ -
749 reinterpret_cast<double*
>(umfpackMatrix_.getValues()),
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
753 UMF_Decomposition_Info);
│ │ │ -
754 Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
│ │ │ -
│ │ │ -
│ │ │ -
757 std::cout <<
"[UMFPack Decomposition]" << std::endl;
│ │ │ -
758 std::cout <<
"Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME] <<
" (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] <<
")" << std::endl;
│ │ │ -
759 std::cout <<
"Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
│ │ │ -
760 std::cout <<
"Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT] <<
" bytes" << std::endl;
│ │ │ -
761 std::cout <<
"Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
│ │ │ -
762 std::cout <<
"Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ] <<
" U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
766 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
770 void printOnApply(
double* UMF_Info)
│ │ │ -
│ │ │ -
772 Caller::report_status(UMF_Control,UMF_Info[UMFPACK_STATUS]);
│ │ │ -
│ │ │ -
│ │ │ -
775 std::cout <<
"[UMFPack Solve]" << std::endl;
│ │ │ -
776 std::cout <<
"Wallclock Time: " << UMF_Info[UMFPACK_SOLVE_WALLTIME] <<
" (CPU Time: " << UMF_Info[UMFPACK_SOLVE_TIME] <<
")" << std::endl;
│ │ │ -
777 std::cout <<
"Flops Taken: " << UMF_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
│ │ │ -
778 std::cout <<
"Iterative Refinement steps taken: " << UMF_Info[UMFPACK_IR_TAKEN] << std::endl;
│ │ │ -
779 std::cout <<
"Error Estimate: " << UMF_Info[UMFPACK_OMEGA1] <<
" resp. " << UMF_Info[UMFPACK_OMEGA2] << std::endl;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
784 bool matrixIsLoaded_;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
788 double UMF_Control[UMFPACK_CONTROL];
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
797 template<
typename T,
typename A>
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
804 template<
class OpTraits,
class=
void>
struct isValidBlock : std::false_type{};
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
807 std::is_same_v<Impl::UMFPackDomainType<typename OpTraits::matrix_type>, typename OpTraits::domain_type>
│ │ │ -
808 && std::is_same_v<Impl::UMFPackRangeType<typename OpTraits::matrix_type>, typename OpTraits::range_type>
│ │ │ -
809 && std::is_same_v<typename FieldTraits<typename OpTraits::domain_type::field_type>::real_type, double>
│ │ │ -
810 && std::is_same_v<typename FieldTraits<typename OpTraits::range_type::field_type>::real_type, double>
│ │ │ -
811 >> : std::true_type {};
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
814 [](
auto opTraits,
const auto& op,
const Dune::ParameterTree& config)
│ │ │ -
815 -> std::shared_ptr<
typename decltype(opTraits)::solver_type>
│ │ │ -
│ │ │ -
817 using OpTraits =
decltype(opTraits);
│ │ │ -
│ │ │ -
819 if constexpr (OpTraits::isParallel){
│ │ │ -
820 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
│ │ │ -
821 DUNE_THROW(Dune::InvalidStateException,
"UMFPack works only for sequential operators.");
│ │ │ -
│ │ │ -
823 if constexpr (OpTraits::isAssembled){
│ │ │ -
824 using M =
typename OpTraits::matrix_type;
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
829 const auto& A = opTraits.getAssembledOpOrThrow(op);
│ │ │ -
830 const M& mat = A->getmat();
│ │ │ -
831 int verbose = config.get(
"verbose", 0);
│ │ │ -
832 return std::make_shared<Dune::UMFPack<M>>(mat,verbose);
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
836 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
Templates characterizing the type of a solver.
│ │ │ -
Implementations of the inverse operator interface.
│ │ │ -
#define DUNE_REGISTER_SOLVER(name,...)
Definition solverregistry.hh:16
│ │ │ -
│ │ │ -
│ │ │ -
│ │ │ -
A dynamic dense block matrix class.
│ │ │ -
│ │ │ -
Implementation of the BCRSMatrix class.
│ │ │ -
│ │ │ -
void setSubMatrix(const Matrix &_mat, const std::set< typename Matrix::size_type > &rowIndexSet)
Definition umfpack.hh:657
│ │ │ -
void free()
free allocated space.
Definition umfpack.hh:715
│ │ │ -
friend class SeqOverlappingSchwarz
Definition umfpack.hh:732
│ │ │ -
SuiteSparse_long size_type
Definition umfpack.hh:269
│ │ │ -
static void symbolic(size_type m, size_type n, const size_type *cs, const size_type *ri, const double *val, A... args)
Definition umfpack.hh:167
│ │ │ -
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)
Definition umfpack.hh:161
│ │ │ -
M matrix_type
Definition umfpack.hh:273
│ │ │ -
static void numeric(const size_type *cs, const size_type *ri, const double *val, A... args)
Definition umfpack.hh:141
│ │ │ -
static void report_info(A... args)
Definition umfpack.hh:146
│ │ │ -
UMFPack(const Matrix &mat_, const ParameterTree &config)
Construct a solver object from a matrix.
Definition umfpack.hh:334
│ │ │ -
static int load_numeric(A... args)
Definition umfpack.hh:136
│ │ │ -
static int load_numeric(A... args)
Definition umfpack.hh:78
│ │ │ -
static void report_status(A... args)
Definition umfpack.hh:93
│ │ │ -
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:358
│ │ │ -
Impl::UMFPackRangeType< M > range_type
The type of the range of the solver.
Definition umfpack.hh:281
│ │ │ -
UMFPack()
default constructor
Definition umfpack.hh:340
│ │ │ -
static void symbolic(A... args)
Definition umfpack.hh:108
│ │ │ -
static void report_info(A... args)
Definition umfpack.hh:88
│ │ │ -
static void free_symbolic(A... args)
Definition umfpack.hh:73
│ │ │ -
static int save_numeric(A... args)
Definition umfpack.hh:156
│ │ │ -
static void free_numeric(A... args)
Definition umfpack.hh:126
│ │ │ -
static int save_numeric(A... args)
Definition umfpack.hh:98
│ │ │ -
static void report_status(A... args)
Definition umfpack.hh:151
│ │ │ -
Impl::UMFPackDomainType< M > domain_type
The type of the domain of the solver.
Definition umfpack.hh:279
│ │ │ -
void apply(T *x, T *b)
additional apply method with c-arrays in analogy to superlu
Definition umfpack.hh:473
│ │ │ -
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category).
Definition umfpack.hh:284
│ │ │ -
static void defaults(A... args)
Definition umfpack.hh:121
│ │ │ -
static void free_numeric(A... args)
Definition umfpack.hh:68
│ │ │ -
void setVerbosity(int v)
sets the verbosity level for the UMFPack solver
Definition umfpack.hh:681
│ │ │ -
UMFPack(const char *file, int verbose=0)
try loading a decomposition from file
Definition umfpack.hh:385
│ │ │ -
static void numeric(A... args)
Definition umfpack.hh:83
│ │ │ -
static constexpr bool valid
Definition umfpack.hh:60
│ │ │ -
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res) override
apply inverse operator, with given convergence criteria.
Definition umfpack.hh:461
│ │ │ -
static constexpr bool valid
Definition umfpack.hh:54
│ │ │ -
virtual ~UMFPack()
Definition umfpack.hh:401
│ │ │ -
const char * name()
Definition umfpack.hh:726
│ │ │ -
void setMatrix(const Matrix &matrix, const BitVector &bitVector={})
Initialize data from given matrix.
Definition umfpack.hh:527
│ │ │ -
void saveDecomposition(const char *file)
saves a decomposition to a file
Definition umfpack.hh:510
│ │ │ -
UMFPackMatrix & getInternalMatrix()
Return the column compress matrix from UMFPack.
Definition umfpack.hh:706
│ │ │ -
SuiteSparse_long size_type
Definition umfpack.hh:118
│ │ │ -
UMFPack(const Matrix &matrix, int verbose=0)
Construct a solver object from a matrix.
Definition umfpack.hh:297
│ │ │ -
ISTL::Impl::BCCSMatrixInitializer< M, size_type > MatrixInitializer
Type of an associated initializer class.
Definition umfpack.hh:277
│ │ │ -
ISTL::Impl::BCCSMatrix< typename Matrix::field_type, size_type > UMFPackMatrix
The corresponding (scalar) UMFPack matrix type.
Definition umfpack.hh:275
│ │ │ -
void setOption(unsigned int option, double value)
Set UMFPack-specific options.
Definition umfpack.hh:499
│ │ │ -
static void free_symbolic(A... args)
Definition umfpack.hh:131
│ │ │ -
M Matrix
The matrix type.
Definition umfpack.hh:272
│ │ │ -
static void defaults(A... args)
Definition umfpack.hh:63
│ │ │ -
void apply(domain_type &x, range_type &b, InverseOperatorResult &res) override
Apply inverse operator,.
Definition umfpack.hh:410
│ │ │ -
static void solve(A... args)
Definition umfpack.hh:103
│ │ │ -
static constexpr bool valid
Definition umfpack.hh:117
│ │ │ -
UMFPack(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition umfpack.hh:315
│ │ │ -
void * getFactorization()
Return the matrix factorization.
Definition umfpack.hh:697
│ │ │ -
@ value
Definition umfpack.hh:794
│ │ │ -
@ value
Definition umfpack.hh:800
│ │ │ -
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
│ │ │ +
#define DUNE_GET_COUNTER(Tag)
Definition counter.hh:17
│ │ │
Definition allocator.hh:11
│ │ │ -
std::pair< std::size_t, std::size_t > flatMatrixForEach(Matrix &&matrix, F &&f, std::size_t rowOffset=0, std::size_t colOffset=0)
Traverse a blocked matrix and call a functor at each scalar entry.
Definition foreach.hh:132
│ │ │ -
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
│ │ │ -
MatrixSparsityPatternGatherScatter< M, I >::ColIter MatrixSparsityPatternGatherScatter< M, I >::col
Definition matrixredistribute.hh:591
│ │ │ -
std::size_t flatVectorForEach(Vector &&vector, F &&f, std::size_t offset=0)
Traverse a blocked vector and call a functor at each scalar entry.
Definition foreach.hh:95
│ │ │ -
Definition umfpack.hh:803
│ │ │ -
static auto coldim(const M &)
Definition matrixutils.hh:219
│ │ │ -
static auto rowdim(const M &)
Definition matrixutils.hh:214
│ │ │ -
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:467
│ │ │ -
derive error class from the base class in common
Definition istlexception.hh:19
│ │ │ -
Sequential overlapping Schwarz preconditioner.
Definition overlappingschwarz.hh:755
│ │ │ -
Definition overlappingschwarz.hh:694
│ │ │ -
Statistics about the application of an inverse operator.
Definition solver.hh:50
│ │ │ -
double elapsed
Elapsed time in seconds.
Definition solver.hh:84
│ │ │ -
int iterations
Number of iterations.
Definition solver.hh:69
│ │ │ -
bool converged
True if convergence criterion has been met.
Definition solver.hh:75
│ │ │ -
Abstract base class for all solvers.
Definition solver.hh:101
│ │ │ -
Category
Definition solvercategory.hh:23
│ │ │ -
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25
│ │ │ -
Definition solverregistry.hh:97
│ │ │ -
Definition solvertype.hh:16
│ │ │ -
Definition solvertype.hh:30
│ │ │ -
│ │ │ -
The UMFPack direct sparse solver.
Definition umfpack.hh:265
│ │ │ -
Definition umfpack.hh:804
│ │ │
1.15.0
│ │ │
│ │ │