3 #ifndef DUNE_ISTL_SUPERLU_HH
4 #define DUNE_ISTL_SUPERLU_HH
16 #include <dune/common/fmatrix.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/stdstreams.hh>
35 template<
class M,
class T,
class TM,
class TD,
class TA>
36 class SeqOverlappingSchwarz;
38 template<
class T,
bool tag>
39 struct SeqOverlappingSchwarzAssemblerHelper;
57 #if __has_include("slu_sdefs.h")
61 static void create(SuperMatrix *
mat,
int n,
int m,
float *dat,
int n1,
62 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
64 sCreate_Dense_Matrix(
mat, n, m, dat, n1, stype, dtype, mtype);
68 static void destroy(SuperMatrix*)
73 struct SuperLUSolveChooser<float>
75 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
76 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
77 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
78 float *rpg,
float *rcond,
float *ferr,
float *berr,
79 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
82 sgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
83 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
84 &gLU, memusage, stat, info);
89 struct QuerySpaceChooser<float>
91 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
93 sQuerySpace(L,U,memusage);
99 #if __has_include("slu_ddefs.h")
102 struct SuperLUDenseMatChooser<double>
104 static void create(SuperMatrix *
mat,
int n,
int m,
double *dat,
int n1,
105 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
107 dCreate_Dense_Matrix(
mat, n, m, dat, n1, stype, dtype, mtype);
111 static void destroy(SuperMatrix * )
115 struct SuperLUSolveChooser<double>
117 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
118 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
119 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
120 double *rpg,
double *rcond,
double *ferr,
double *berr,
121 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
124 dgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
125 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
126 &gLU, memusage, stat, info);
131 struct QuerySpaceChooser<double>
133 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
135 dQuerySpace(L,U,memusage);
140 #if __has_include("slu_zdefs.h")
142 struct SuperLUDenseMatChooser<std::complex<double> >
144 static void create(SuperMatrix *
mat,
int n,
int m, std::complex<double> *dat,
int n1,
145 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
147 zCreate_Dense_Matrix(
mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
151 static void destroy(SuperMatrix*)
156 struct SuperLUSolveChooser<std::complex<double> >
158 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
159 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
160 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
161 double *rpg,
double *rcond,
double *ferr,
double *berr,
162 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
165 zgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
166 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
167 &gLU, memusage, stat, info);
172 struct QuerySpaceChooser<std::complex<double> >
174 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
176 zQuerySpace(L,U,memusage);
181 #if __has_include("slu_cdefs.h")
183 struct SuperLUDenseMatChooser<std::complex<float> >
185 static void create(SuperMatrix *
mat,
int n,
int m, std::complex<float> *dat,
int n1,
186 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
188 cCreate_Dense_Matrix(
mat, n, m,
reinterpret_cast< ::complex*
>(dat), n1, stype, dtype, mtype);
192 static void destroy(SuperMatrix* )
197 struct SuperLUSolveChooser<std::complex<float> >
199 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
200 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
201 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
202 float *rpg,
float *rcond,
float *ferr,
float *berr,
203 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
206 cgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
207 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
208 &gLU, memusage, stat, info);
213 struct QuerySpaceChooser<std::complex<float> >
215 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
217 cQuerySpace(L,U,memusage);
225 struct SuperLUVectorChooser
228 template<
typename T,
typename A,
int n,
int m>
229 struct SuperLUVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
232 using domain_type = BlockVector<
234 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >;
236 using range_type = BlockVector<
238 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > >;
241 template<
typename T,
typename A>
242 struct SuperLUVectorChooser<BCRSMatrix<T,A> >
245 using domain_type = BlockVector<T, A>;
247 using range_type = BlockVector<T, A>;
267 typename Impl::SuperLUVectorChooser<M>::domain_type,
268 typename Impl::SuperLUVectorChooser<M>::range_type >
270 using T =
typename M::field_type;
280 using domain_type =
typename Impl::SuperLUVectorChooser<M>::domain_type;
282 using range_type =
typename Impl::SuperLUVectorChooser<M>::range_type;
287 return SolverCategory::Category::sequential;
305 bool reusevector=
true);
319 :
SuperLU(
mat, config.
get<bool>(
"verbose", false), config.
get<bool>(
"reuseVector", true))
348 void apply(T* x, T* b);
353 typename SuperLUMatrix::size_type
nnz()
const
355 return mat.nonzeroes();
369 const char*
name() {
return "SuperLU"; }
371 template<
class Mat,
class X,
class TM,
class TD,
class T1>
381 SuperMatrix L, U, B, X;
382 int *perm_c, *perm_r, *etree;
385 superlu_options_t options;
389 bool first, verbose, reusevector;
409 Destroy_SuperNode_Matrix(&L);
410 Destroy_CompCol_Matrix(&U);
413 if(!first && reusevector) {
414 SUPERLU_FREE(B.Store);
415 SUPERLU_FREE(X.Store);
423 : work(0), lwork(0), first(true), verbose(verbose_),
424 reusevector(reusevector_)
431 : work(0), lwork(0),verbose(false),
464 mat.setMatrix(mat_,mrs);
473 perm_c =
new int[
mat.
M()];
474 perm_r =
new int[
mat.
N()];
475 etree =
new int[
mat.
M()];
479 set_default_options(&options);
486 fakeFormat.lda=
mat.
N();
496 mem_usage_t memusage;
501 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
502 &berr, &memusage, &stat, &info);
505 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
507 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
509 if ( info == 0 || info == nSuperLUCol+1 ) {
511 if ( options.PivotGrowth )
512 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
513 if ( options.ConditionNumber )
514 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
515 SCformat* Lstore = (SCformat *) L.Store;
516 NCformat* Ustore = (NCformat *) U.Store;
517 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
518 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
519 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - nSuperLUCol<<std::endl;
521 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
523 std::cout<<stat.expansions<<std::endl;
525 }
else if ( info > 0 && lwork == -1 ) {
526 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<std::endl;
528 if ( options.PrintStat ) StatPrint(&stat);
556 options.Fact = FACTORED;
563 if (
mat.
N() != b.dim())
564 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
565 if (
mat.
M() != x.dim())
566 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
568 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
570 SuperMatrix* mB = &B;
571 SuperMatrix* mX = &X;
579 ((DNformat*)B.Store)->nzval=&b[0];
580 ((DNformat*)X.Store)->nzval=&x[0];
590 mem_usage_t memusage;
600 options.IterRefine=SLU_DOUBLE;
603 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
604 &memusage, &stat, &info);
624 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
626 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
628 if ( info == 0 || info == nSuperLUCol+1 ) {
630 if ( options.IterRefine ) {
631 std::cout<<
"Iterative Refinement: steps="
632 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
634 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
635 }
else if ( info > 0 && lwork == -1 ) {
636 std::cout<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
639 if ( options.PrintStat ) StatPrint(&stat);
643 SUPERLU_FREE(rB.Store);
644 SUPERLU_FREE(rX.Store);
653 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
655 SuperMatrix* mB = &B;
656 SuperMatrix* mX = &X;
664 ((DNformat*) B.Store)->nzval=b;
665 ((DNformat*)X.Store)->nzval=x;
676 mem_usage_t memusage;
681 options.IterRefine=SLU_DOUBLE;
684 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
685 &memusage, &stat, &info);
688 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
690 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
692 if ( info == 0 || info == nSuperLUCol+1 ) {
694 if ( options.IterRefine ) {
695 dinfo<<
"Iterative Refinement: steps="
696 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
698 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
699 }
else if ( info > 0 && lwork == -1 ) {
700 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
702 if ( options.PrintStat ) StatPrint(&stat);
707 SUPERLU_FREE(rB.Store);
708 SUPERLU_FREE(rX.Store);
713 template<
typename T,
typename A>
719 template<
typename T,
typename A>
728 template<
int k>
struct isValidBlock<
Dune::FieldVector<std::complex<double>,k>> : std::true_type{};
729 template<
typename TL,
typename M>
730 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
731 typename Dune::TypeListElement<2, TL>::type>>
733 std::enable_if_t<
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
735 int verbose = config.get(
"verbose", 0);
736 return std::make_shared<Dune::SuperLU<M>>(
mat,verbose);
740 template<
typename TL,
typename M>
741 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
742 typename Dune::TypeListElement<2, TL>::type>>
744 std::enable_if_t<!
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
747 "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
757 #undef FIRSTCOL_OF_SNODE
762 #undef SUPERLU_MALLOC
782 #undef DROP_SECONDARY
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
void setSubMatrix(const Matrix &mat, const S &rowIndexSet)
Definition: superlu.hh:455
void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: superlu.hh:561
DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator())
void setVerbosity(bool v)
Definition: superlu.hh:435
void free()
free allocated space.
Definition: superlu.hh:401
~SuperLU()
Definition: superlu.hh:394
SuperLU()
Empty default constructor.
Definition: superlu.hh:430
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:441
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
derive error class from the base class in common
Definition: istlexception.hh:17
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
Definition: overlappingschwarz.hh:692
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
Statistics about the application of an inverse operator.
Definition: solver.hh:46
int iterations
Number of iterations.
Definition: solver.hh:65
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Abstract base class for all solvers.
Definition: solver.hh:97
Category
Definition: solvercategory.hh:21
Definition: solverregistry.hh:75
Definition: solvertype.hh:14
@ value
Whether this is a direct solver.
Definition: solvertype.hh:22
Definition: solvertype.hh:28
@ value
whether the solver internally uses column compressed storage
Definition: solvertype.hh:34
Definition: superlu.hh:43
Definition: superlu.hh:47
Definition: superlu.hh:51
Definition: superlu.hh:55
SuperLu Solver.
Definition: superlu.hh:269
void apply(domain_type &x, range_type &b, [[maybe_unused]] double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:340
SuperLUMatrix::size_type nnz() const
Definition: superlu.hh:353
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:282
M matrix_type
Definition: superlu.hh:274
SuperMatrixInitializer< Matrix > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:278
M Matrix
The matrix type.
Definition: superlu.hh:273
typename Impl::SuperLUVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: superlu.hh:280
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:285
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:276
SuperLU(const Matrix &mat, const ParameterTree &config)
Constructs the SuperLU solver.
Definition: superlu.hh:318
const char * name()
Definition: superlu.hh:369
Definition: superlu.hh:725
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: superlu.hh:732
Definition: superlu.hh:726
Definition: supermatrix.hh:130
Definition: supermatrix.hh:177