3#ifndef DUNE_ISTL_SUPERLU_HH
4#define DUNE_ISTL_SUPERLU_HH
8#include "superlufunctions.hh"
10#include "supermatrix.hh"
15#include "istlexception.hh"
20#include <dune/istl/solverfactory.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;
42 struct SuperLUSolveChooser
46 struct SuperLUDenseMatChooser
50 struct SuperLUQueryChooser
54 struct QuerySpaceChooser
57#if __has_include("slu_sdefs.h")
59 struct SuperLUDenseMatChooser<float>
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;
274 using matrix_type = M;
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();
359 void setSubMatrix(
const Matrix& mat,
const S& rowIndexSet);
361 void setVerbosity(
bool v);
369 const char* name() {
return "SuperLU"; }
371 template<
class Mat,
class X,
class TM,
class TD,
class T1>
372 friend class SeqOverlappingSchwarz;
373 friend struct SeqOverlappingSchwarzAssemblerHelper<
SuperLU<
Matrix>,true>;
381 SuperMatrix L, U, B, X;
382 int *perm_c, *perm_r, *etree;
383 typename GetSuperLUType<T>::float_type *R, *C;
385 superlu_options_t options;
389 bool first, verbose, reusevector;
396 if(mat.
N()+mat.
M()>0)
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),
443 if(mat.
N()+mat.
M()>0) {
458 if(mat.
N()+mat.
M()>0) {
464 mat.setMatrix(mat_,mrs);
469 void SuperLU<M>::decompose()
473 perm_c =
new int[mat.
M()];
474 perm_r =
new int[mat.
N()];
475 etree =
new int[mat.
M()];
476 R =
new typename GetSuperLUType<T>::float_type[mat.
N()];
477 C =
new typename GetSuperLUType<T>::float_type[mat.
M()];
479 set_default_options(&options);
483 B.Dtype=GetSuperLUType<T>::type;
486 fakeFormat.lda=mat.
N();
489 X.Dtype=GetSuperLUType<T>::type;
494 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
496 mem_usage_t memusage;
500 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
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;
520 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
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!");
567 if (mat.
M()+mat.
N()==0)
570 SuperMatrix* mB = &B;
571 SuperMatrix* mX = &X;
575 SuperLUDenseMatChooser<T>::create(&B, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
576 SuperLUDenseMatChooser<T>::create(&X, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
579 ((DNformat*)B.Store)->nzval=&b[0];
580 ((DNformat*)X.Store)->nzval=&x[0];
583 SuperLUDenseMatChooser<T>::create(&rB, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
584 SuperLUDenseMatChooser<T>::create(&rX, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
588 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
590 mem_usage_t memusage;
600 options.IterRefine=SLU_DOUBLE;
602 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
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);
652 if(mat.
N()+mat.
M()==0)
655 SuperMatrix* mB = &B;
656 SuperMatrix* mX = &X;
660 SuperLUDenseMatChooser<T>::create(&B, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
661 SuperLUDenseMatChooser<T>::create(&X, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
664 ((DNformat*) B.Store)->nzval=b;
665 ((DNformat*)X.Store)->nzval=x;
668 SuperLUDenseMatChooser<T>::create(&rB, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
669 SuperLUDenseMatChooser<T>::create(&rX, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
674 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
676 mem_usage_t memusage;
681 options.IterRefine=SLU_DOUBLE;
683 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
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>
720 struct StoresColumnCompressed<SuperLU<BCRSMatrix<T,A> > >
722 enum { value =
true };
725 struct SuperLUCreator {
726 template<
class>
struct isValidBlock : std::false_type{};
727 template<
int k>
struct isValidBlock<
Dune::FieldVector<double,k>> : std::true_type{};
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)");
750 template<>
struct SuperLUCreator::isValidBlock<double> : std::true_type{};
751 template<>
struct SuperLUCreator::isValidBlock<
std::complex<double>> : std::true_type{};
753 DUNE_REGISTER_DIRECT_SOLVER(
"superlu", SuperLUCreator());
757#undef FIRSTCOL_OF_SNODE
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
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
Abstract base class for all solvers.
Definition: solver.hh:97
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
Hierarchical structure of string parameters.
Definition: parametertree.hh:35
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:183
SuperLu Solver.
Definition: superlu.hh:269
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:340
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:282
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
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:216
void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: superlu.hh:561
void free()
free allocated space.
Definition: superlu.hh:401
SuperLU()
Empty default constructor.
Definition: superlu.hh:430
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:441
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune namespace.
Definition: alignedallocator.hh:11
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
Standard Dune debug streams.
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
Category
Definition: solvercategory.hh:21