5#ifndef DUNE_ISTL_SUPERLU_HH
6#define DUNE_ISTL_SUPERLU_HH
10#include "superlufunctions.hh"
12#include "supermatrix.hh"
15#include "bcrsmatrix.hh"
17#include "istlexception.hh"
22#include <dune/istl/solverfactory.hh>
37 template<
class M,
class T,
class TM,
class TD,
class TA>
38 class SeqOverlappingSchwarz;
40 template<
class T,
bool tag>
41 struct SeqOverlappingSchwarzAssemblerHelper;
44 struct SuperLUSolveChooser
48 struct SuperLUDenseMatChooser
52 struct SuperLUQueryChooser
56 struct QuerySpaceChooser
59#if __has_include("slu_sdefs.h")
61 struct SuperLUDenseMatChooser<float>
63 static void create(SuperMatrix *mat,
int n,
int m,
float *dat,
int n1,
64 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
66 sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
70 static void destroy(SuperMatrix*)
75 struct SuperLUSolveChooser<float>
77 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
78 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
79 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
80 float *rpg,
float *rcond,
float *ferr,
float *berr,
81 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
84 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
85 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
86 &gLU, memusage, stat, info);
91 struct QuerySpaceChooser<float>
93 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
95 sQuerySpace(L,U,memusage);
101#if __has_include("slu_ddefs.h")
104 struct SuperLUDenseMatChooser<double>
106 static void create(SuperMatrix *mat,
int n,
int m,
double *dat,
int n1,
107 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
109 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
113 static void destroy(SuperMatrix * )
117 struct SuperLUSolveChooser<double>
119 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
120 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
121 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
122 double *rpg,
double *rcond,
double *ferr,
double *berr,
123 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
126 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
127 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
128 &gLU, memusage, stat, info);
133 struct QuerySpaceChooser<double>
135 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
137 dQuerySpace(L,U,memusage);
142#if __has_include("slu_zdefs.h")
144 struct SuperLUDenseMatChooser<
std::complex<double> >
146 static void create(SuperMatrix *mat,
int n,
int m, std::complex<double> *dat,
int n1,
147 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
149 zCreate_Dense_Matrix(mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
153 static void destroy(SuperMatrix*)
158 struct SuperLUSolveChooser<
std::complex<double> >
160 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
161 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
162 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
163 double *rpg,
double *rcond,
double *ferr,
double *berr,
164 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
167 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
168 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
169 &gLU, memusage, stat, info);
174 struct QuerySpaceChooser<
std::complex<double> >
176 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
178 zQuerySpace(L,U,memusage);
183#if __has_include("slu_cdefs.h")
185 struct SuperLUDenseMatChooser<
std::complex<float> >
187 static void create(SuperMatrix *mat,
int n,
int m, std::complex<float> *dat,
int n1,
188 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
190 cCreate_Dense_Matrix(mat, n, m,
reinterpret_cast< ::complex*
>(dat), n1, stype, dtype, mtype);
194 static void destroy(SuperMatrix* )
199 struct SuperLUSolveChooser<
std::complex<float> >
201 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
202 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
203 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
204 float *rpg,
float *rcond,
float *ferr,
float *berr,
205 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
208 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
209 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
210 &gLU, memusage, stat, info);
215 struct QuerySpaceChooser<
std::complex<float> >
217 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
219 cQuerySpace(L,U,memusage);
227 struct SuperLUVectorChooser
230 template<
typename T,
typename A,
int n,
int m>
231 struct SuperLUVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
234 using domain_type = BlockVector<
236 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >;
238 using range_type = BlockVector<
240 typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > >;
243 template<
typename T,
typename A>
244 struct SuperLUVectorChooser<BCRSMatrix<T,A> >
247 using domain_type = BlockVector<T, A>;
249 using range_type = BlockVector<T, A>;
269 typename Impl::SuperLUVectorChooser<M>::domain_type,
270 typename Impl::SuperLUVectorChooser<M>::range_type >
272 using T =
typename M::field_type;
276 using matrix_type = M;
282 using domain_type =
typename Impl::SuperLUVectorChooser<M>::domain_type;
284 using range_type =
typename Impl::SuperLUVectorChooser<M>::range_type;
289 return SolverCategory::Category::sequential;
307 bool reusevector=
true);
321 :
SuperLU(mat, config.
get<bool>(
"verbose", false), config.
get<bool>(
"reuseVector", true))
355 typename SuperLUMatrix::size_type nnz()
const
357 return mat.nonzeroes();
361 void setSubMatrix(
const Matrix& mat,
const S& rowIndexSet);
363 void setVerbosity(
bool v);
371 const char* name() {
return "SuperLU"; }
373 template<
class Mat,
class X,
class TM,
class TD,
class T1>
375 friend struct SeqOverlappingSchwarzAssemblerHelper<
SuperLU<
Matrix>,true>;
383 SuperMatrix L, U, B, X;
384 int *perm_c, *perm_r, *etree;
385 typename GetSuperLUType<T>::float_type *R, *C;
387 superlu_options_t options;
391 bool first, verbose, reusevector;
398 if(mat.
N()+mat.
M()>0)
411 Destroy_SuperNode_Matrix(&L);
412 Destroy_CompCol_Matrix(&U);
415 if(!first && reusevector) {
416 SUPERLU_FREE(B.Store);
417 SUPERLU_FREE(X.Store);
425 : work(0), lwork(0), first(true), verbose(verbose_),
426 reusevector(reusevector_)
433 : work(0), lwork(0),verbose(false),
445 if(mat.
N()+mat.
M()>0) {
460 if(mat.
N()+mat.
M()>0) {
466 mat.setMatrix(mat_,mrs);
471 void SuperLU<M>::decompose()
475 perm_c =
new int[mat.
M()];
476 perm_r =
new int[mat.
N()];
477 etree =
new int[mat.
M()];
478 R =
new typename GetSuperLUType<T>::float_type[mat.
N()];
479 C =
new typename GetSuperLUType<T>::float_type[mat.
M()];
481 set_default_options(&options);
485 B.Dtype=GetSuperLUType<T>::type;
488 fakeFormat.lda=mat.
N();
491 X.Dtype=GetSuperLUType<T>::type;
496 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
498 mem_usage_t memusage;
502 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
503 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
504 &berr, &memusage, &stat, &info);
507 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
509 auto nSuperLUCol =
static_cast<SuperMatrix&
>(mat).ncol;
511 if ( info == 0 || info == nSuperLUCol+1 ) {
513 if ( options.PivotGrowth )
514 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
515 if ( options.ConditionNumber )
516 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
517 SCformat* Lstore = (SCformat *) L.Store;
518 NCformat* Ustore = (NCformat *) U.Store;
519 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
520 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
521 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - nSuperLUCol<<std::endl;
522 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
523 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
525 std::cout<<stat.expansions<<std::endl;
527 }
else if ( info > 0 && lwork == -1 ) {
528 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<std::endl;
530 if ( options.PrintStat ) StatPrint(&stat);
558 options.Fact = FACTORED;
565 if (mat.
N() != b.dim())
566 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
567 if (mat.
M() != x.dim())
568 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
569 if (mat.
M()+mat.
N()==0)
572 SuperMatrix* mB = &B;
573 SuperMatrix* mX = &X;
577 SuperLUDenseMatChooser<T>::create(&B, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
578 SuperLUDenseMatChooser<T>::create(&X, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
581 ((DNformat*)B.Store)->nzval=&b[0];
582 ((DNformat*)X.Store)->nzval=&x[0];
585 SuperLUDenseMatChooser<T>::create(&rB, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
586 SuperLUDenseMatChooser<T>::create(&rX, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
590 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
592 mem_usage_t memusage;
602 options.IterRefine=SLU_DOUBLE;
604 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
605 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
606 &memusage, &stat, &info);
626 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
628 auto nSuperLUCol =
static_cast<SuperMatrix&
>(mat).ncol;
630 if ( info == 0 || info == nSuperLUCol+1 ) {
632 if ( options.IterRefine ) {
633 std::cout<<
"Iterative Refinement: steps="
634 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
636 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
637 }
else if ( info > 0 && lwork == -1 ) {
638 std::cout<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
641 if ( options.PrintStat ) StatPrint(&stat);
645 SUPERLU_FREE(rB.Store);
646 SUPERLU_FREE(rX.Store);
654 if(mat.
N()+mat.
M()==0)
657 SuperMatrix* mB = &B;
658 SuperMatrix* mX = &X;
662 SuperLUDenseMatChooser<T>::create(&B, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
663 SuperLUDenseMatChooser<T>::create(&X, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
666 ((DNformat*) B.Store)->nzval=b;
667 ((DNformat*)X.Store)->nzval=x;
670 SuperLUDenseMatChooser<T>::create(&rB, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
671 SuperLUDenseMatChooser<T>::create(&rX, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
676 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
678 mem_usage_t memusage;
683 options.IterRefine=SLU_DOUBLE;
685 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
686 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
687 &memusage, &stat, &info);
690 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
692 auto nSuperLUCol =
static_cast<SuperMatrix&
>(mat).ncol;
694 if ( info == 0 || info == nSuperLUCol+1 ) {
696 if ( options.IterRefine ) {
697 dinfo<<
"Iterative Refinement: steps="
698 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
700 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
701 }
else if ( info > 0 && lwork == -1 ) {
702 dinfo<<
"** Estimated memory: "<< info - nSuperLUCol<<
" bytes"<<std::endl;
704 if ( options.PrintStat ) StatPrint(&stat);
709 SUPERLU_FREE(rB.Store);
710 SUPERLU_FREE(rX.Store);
715 template<
typename T,
typename A>
721 template<
typename T,
typename A>
722 struct StoresColumnCompressed<SuperLU<BCRSMatrix<T,A> > >
724 enum { value =
true };
728 DUNE_REGISTER_SOLVER(
"superlu",
730 -> std::shared_ptr<
typename decltype(opTraits)::solver_type>
732 using OpTraits =
decltype(opTraits);
733 using M =
typename OpTraits::matrix_type;
735 if constexpr (OpTraits::isParallel){
736 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
739 if constexpr (OpTraits::isAssembled &&
741 std::is_same_v<typename FieldTraits<M>::real_type,
double>){
745 if constexpr (std::is_convertible_v<SuperLU<M>*,
747 typename OpTraits::range_type>*>
749 const auto& A = opTraits.getAssembledOpOrThrow(op);
750 const M& mat = A->getmat();
751 int verbose = config.
get(
"verbose", 0);
752 return std::make_shared<Dune::SuperLU<M>>(mat,verbose);
756 "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
761#undef FIRSTCOL_OF_SNODE
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:466
derive error class from the base class in common
Definition: istlexception.hh:19
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Abstract base class for all solvers.
Definition: solver.hh:101
Impl::SuperLUVectorChooser< ISTLM >::range_type range_type
Definition: solver.hh:107
Impl::SuperLUVectorChooser< ISTLM >::domain_type domain_type
Definition: solver.hh:104
A generic dynamic dense matrix.
Definition: matrix.hh:561
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
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:188
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:755
SuperLu Solver.
Definition: superlu.hh:271
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res) override
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:342
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:284
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:287
SuperMatrixInitializer< Matrix > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:280
M Matrix
The matrix type.
Definition: superlu.hh:275
typename Impl::SuperLUVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: superlu.hh:282
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:278
SuperLU(const Matrix &mat, const ParameterTree &config)
Constructs the SuperLU solver.
Definition: superlu.hh:320
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:218
SuperLU(const Matrix &mat, bool verbose=false, bool reusevector=true)
Constructs the SuperLU solver.
Definition: superlu.hh:424
void apply(T *x, T *b)
Apply SuperLu to C arrays.
Definition: superlu.hh:652
void apply(domain_type &x, range_type &b, InverseOperatorResult &res) override
Apply inverse operator,.
Definition: superlu.hh:563
void free()
free allocated space.
Definition: superlu.hh:403
SuperLU()
Empty default constructor.
Definition: superlu.hh:432
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:443
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:141
Dune namespace.
Definition: alignedallocator.hh:13
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
Templates characterizing the type of a solver.
Standard Dune debug streams.
Statistics about the application of an inverse operator.
Definition: solver.hh:50
int iterations
Number of iterations.
Definition: solver.hh:69
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Category
Definition: solvercategory.hh:23