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;
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)
81 #if SUPERLU_MIN_VERSION_5
83 sgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
84 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
85 &gLU, memusage, stat, info);
87 sgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
88 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
89 memusage, stat, info);
95 struct QuerySpaceChooser<float>
97 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
99 sQuerySpace(L,U,memusage);
108 struct SuperLUDenseMatChooser<double>
110 static void create(SuperMatrix *
mat,
int n,
int m,
double *dat,
int n1,
111 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
113 dCreate_Dense_Matrix(
mat, n, m, dat, n1, stype, dtype, mtype);
117 static void destroy(SuperMatrix * )
121 struct SuperLUSolveChooser<double>
123 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
124 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
125 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
126 double *rpg,
double *rcond,
double *ferr,
double *berr,
127 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
129 #if SUPERLU_MIN_VERSION_5
131 dgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
132 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
133 &gLU, memusage, stat, info);
135 dgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
136 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
137 memusage, stat, info);
143 struct QuerySpaceChooser<double>
145 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
147 dQuerySpace(L,U,memusage);
154 struct SuperLUDenseMatChooser<std::complex<double> >
156 static void create(SuperMatrix *
mat,
int n,
int m, std::complex<double> *dat,
int n1,
157 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
159 zCreate_Dense_Matrix(
mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
163 static void destroy(SuperMatrix*)
168 struct SuperLUSolveChooser<std::complex<double> >
170 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
171 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
172 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
173 double *rpg,
double *rcond,
double *ferr,
double *berr,
174 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
176 #if SUPERLU_MIN_VERSION_5
178 zgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
179 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
180 &gLU, memusage, stat, info);
182 zgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
183 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
184 memusage, stat, info);
190 struct QuerySpaceChooser<std::complex<double> >
192 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
194 zQuerySpace(L,U,memusage);
201 struct SuperLUDenseMatChooser<std::complex<float> >
203 static void create(SuperMatrix *
mat,
int n,
int m, std::complex<float> *dat,
int n1,
204 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
206 cCreate_Dense_Matrix(
mat, n, m,
reinterpret_cast< ::complex*
>(dat), n1, stype, dtype, mtype);
210 static void destroy(SuperMatrix* )
215 struct SuperLUSolveChooser<std::complex<float> >
217 static void solve(superlu_options_t *options, SuperMatrix *
mat,
int *perm_c,
int *perm_r,
int *etree,
218 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
219 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
220 float *rpg,
float *rcond,
float *ferr,
float *berr,
221 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
223 #if SUPERLU_MIN_VERSION_5
225 cgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
226 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
227 &gLU, memusage, stat, info);
229 cgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
230 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
231 memusage, stat, info);
237 struct QuerySpaceChooser<std::complex<float> >
239 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
241 cQuerySpace(L,U,memusage);
249 struct SuperLUVectorChooser
252 template<
typename T,
typename A,
int n,
int m>
253 struct SuperLUVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
256 using domain_type = BlockVector<
258 typename A::template rebind<FieldVector<T,m> >::other>;
260 using range_type = BlockVector<
262 typename A::template rebind<FieldVector<T,n> >::other>;
265 template<
typename T,
typename A>
266 struct SuperLUVectorChooser<BCRSMatrix<T,A> >
269 using domain_type = BlockVector<T, A>;
271 using range_type = BlockVector<T, A>;
291 typename Impl::SuperLUVectorChooser<M>::domain_type,
292 typename Impl::SuperLUVectorChooser<M>::range_type >
294 using T =
typename M::field_type;
304 using domain_type =
typename Impl::SuperLUVectorChooser<M>::domain_type;
306 using range_type =
typename Impl::SuperLUVectorChooser<M>::range_type;
311 return SolverCategory::Category::sequential;
329 bool reusevector=
true);
350 DUNE_UNUSED_PARAMETER(reduction);
357 void apply(T* x, T* b);
362 typename SuperLUMatrix::size_type
nnz()
const
378 const char*
name() {
return "SuperLU"; }
380 template<
class Mat,
class X,
class TM,
class TD,
class T1>
390 SuperMatrix L, U, B, X;
391 int *perm_c, *perm_r, *etree;
394 superlu_options_t options;
398 bool first, verbose, reusevector;
418 Destroy_SuperNode_Matrix(&L);
419 Destroy_CompCol_Matrix(&U);
422 if(!first && reusevector) {
423 SUPERLU_FREE(B.Store);
424 SUPERLU_FREE(X.Store);
432 : work(0), lwork(0), first(true), verbose(verbose_),
433 reusevector(reusevector_)
440 : work(0), lwork(0),verbose(false),
473 mat.setMatrix(mat_,mrs);
482 perm_c =
new int[
mat.
M()];
483 perm_r =
new int[
mat.
N()];
484 etree =
new int[
mat.
M()];
488 set_default_options(&options);
495 fakeFormat.lda=
mat.
N();
505 mem_usage_t memusage;
510 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
511 &berr, &memusage, &stat, &info);
514 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
516 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
518 if ( info == 0 || info == nSuperLUCol+1 ) {
520 if ( options.PivotGrowth )
521 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
522 if ( options.ConditionNumber )
523 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
524 SCformat* Lstore = (SCformat *) L.Store;
525 NCformat* Ustore = (NCformat *) U.Store;
526 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
527 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
528 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - nSuperLUCol<<std::endl;
530 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
532 std::cout<<stat.expansions<<std::endl;
534 }
else if ( info > 0 && lwork == -1 ) {
535 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<std::endl;
537 if ( options.PrintStat ) StatPrint(&stat);
565 options.Fact = FACTORED;
572 if (
mat.
N() != b.dim())
573 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
574 if (
mat.
M() != x.dim())
575 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
577 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
579 SuperMatrix* mB = &B;
580 SuperMatrix* mX = &X;
588 ((DNformat*)B.Store)->nzval=&b[0];
589 ((DNformat*)X.Store)->nzval=&x[0];
599 mem_usage_t memusage;
609 #ifdef SUPERLU_MIN_VERSION_4_3
610 options.IterRefine=SLU_DOUBLE;
612 options.IterRefine=DOUBLE;
616 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
617 &memusage, &stat, &info);
637 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
639 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
641 if ( info == 0 || info == nSuperLUCol+1 ) {
643 if ( options.IterRefine ) {
644 std::cout<<
"Iterative Refinement: steps="
645 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
647 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
648 }
else if ( info > 0 && lwork == -1 ) {
649 std::cout<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
652 if ( options.PrintStat ) StatPrint(&stat);
656 SUPERLU_FREE(rB.Store);
657 SUPERLU_FREE(rX.Store);
666 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
668 SuperMatrix* mB = &B;
669 SuperMatrix* mX = &X;
677 ((DNformat*) B.Store)->nzval=b;
678 ((DNformat*)X.Store)->nzval=x;
689 mem_usage_t memusage;
694 #ifdef SUPERLU_MIN_VERSION_4_3
695 options.IterRefine=SLU_DOUBLE;
697 options.IterRefine=DOUBLE;
701 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
702 &memusage, &stat, &info);
705 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
707 auto nSuperLUCol =
static_cast<SuperMatrix&
>(
mat).ncol;
709 if ( info == 0 || info == nSuperLUCol+1 ) {
711 if ( options.IterRefine ) {
712 dinfo<<
"Iterative Refinement: steps="
713 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
715 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
716 }
else if ( info > 0 && lwork == -1 ) {
717 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
719 if ( options.PrintStat ) StatPrint(&stat);
724 SUPERLU_FREE(rB.Store);
725 SUPERLU_FREE(rX.Store);
730 template<
typename T,
typename A>
736 template<
typename T,
typename A>
745 template<
int k>
struct isValidBlock<
Dune::FieldVector<std::complex<double>,k>> : std::true_type{};
746 template<
typename TL,
typename M>
747 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
748 typename Dune::TypeListElement<2, TL>::type>>
750 std::enable_if_t<
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
752 int verbose = config.get(
"verbose", 0);
753 return std::make_shared<Dune::SuperLU<M>>(
mat,verbose);
757 template<
typename TL,
typename M>
758 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
759 typename Dune::TypeListElement<2, TL>::type>>
761 std::enable_if_t<!
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
764 "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
774 #undef FIRSTCOL_OF_SNODE
779 #undef SUPERLU_MALLOC
799 #undef DROP_SECONDARY
804 #endif // HAVE_SUPERLU
805 #endif // DUNE_SUPERLU_HH