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 std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >;
260 using range_type = BlockVector<
262 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > >;
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
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:464
void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: superlu.hh:570
void setVerbosity(bool v)
Definition: superlu.hh:444
void free()
free allocated space.
Definition: superlu.hh:410
~SuperLU()
Definition: superlu.hh:403
SuperLU()
Empty default constructor.
Definition: superlu.hh:439
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:450
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:7
DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator())
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:426
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:752
Definition: overlappingschwarz.hh:691
derive error class from the base class in common
Definition: istlexception.hh:16
A generic dynamic dense matrix.
Definition: matrix.hh:558
size_type M() const
Return the number of columns.
Definition: matrix.hh:696
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
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: solverfactory.hh:124
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:293
SuperLUMatrix::size_type nnz() const
Definition: superlu.hh:362
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:348
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:306
M matrix_type
Definition: superlu.hh:298
SuperMatrixInitializer< Matrix > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:302
M Matrix
The matrix type.
Definition: superlu.hh:297
typename Impl::SuperLUVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: superlu.hh:304
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:309
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:300
const char * name()
Definition: superlu.hh:378
Definition: superlu.hh:742
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:749
Definition: superlu.hh:743
Definition: supermatrix.hh:130
Definition: supermatrix.hh:177