3#ifndef DUNE_ISTL_SUPERLU_HH
4#define DUNE_ISTL_SUPERLU_HH
8#include "superlufunctions.hh"
10#include "supermatrix.hh"
15#include "istlexception.hh"
34 template<
class Matrix>
38 template<
class M,
class T,
class TM,
class TD,
class TA>
39 class SeqOverlappingSchwarz;
41 template<
class T,
bool tag>
42 struct SeqOverlappingSchwarzAssemblerHelper;
45 struct SuperLUSolveChooser
49 struct SuperLUDenseMatChooser
53 struct SuperLUQueryChooser
57 struct QuerySpaceChooser
62 struct SuperLUDenseMatChooser<float>
64 static void create(SuperMatrix *mat,
int n,
int m,
float *dat,
int n1,
65 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
67 sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
71 static void destroy(SuperMatrix*)
76 struct SuperLUSolveChooser<float>
78 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
79 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
80 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
81 float *rpg,
float *rcond,
float *ferr,
float *berr,
82 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
84#if SUPERLU_MIN_VERSION_5
86 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
87 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
88 &gLU, memusage, stat, info);
90 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
91 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
92 memusage, stat, info);
98 struct QuerySpaceChooser<float>
100 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
102 sQuerySpace(L,U,memusage);
111 struct SuperLUDenseMatChooser<double>
113 static void create(SuperMatrix *mat,
int n,
int m,
double *dat,
int n1,
114 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
116 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
120 static void destroy(SuperMatrix * )
124 struct SuperLUSolveChooser<double>
126 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
127 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
128 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
129 double *rpg,
double *rcond,
double *ferr,
double *berr,
130 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
132#if SUPERLU_MIN_VERSION_5
134 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
135 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
136 &gLU, memusage, stat, info);
138 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
139 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
140 memusage, stat, info);
146 struct QuerySpaceChooser<double>
148 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
150 dQuerySpace(L,U,memusage);
157 struct SuperLUDenseMatChooser<
std::complex<double> >
159 static void create(SuperMatrix *mat,
int n,
int m, std::complex<double> *dat,
int n1,
160 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
162 zCreate_Dense_Matrix(mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
166 static void destroy(SuperMatrix*)
171 struct SuperLUSolveChooser<
std::complex<double> >
173 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
174 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
175 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
176 double *rpg,
double *rcond,
double *ferr,
double *berr,
177 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
179#if SUPERLU_MIN_VERSION_5
181 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
182 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
183 &gLU, memusage, stat, info);
185 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
186 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
187 memusage, stat, info);
193 struct QuerySpaceChooser<
std::complex<double> >
195 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
197 zQuerySpace(L,U,memusage);
204 struct SuperLUDenseMatChooser<
std::complex<float> >
206 static void create(SuperMatrix *mat,
int n,
int m, std::complex<float> *dat,
int n1,
207 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
209 cCreate_Dense_Matrix(mat, n, m,
reinterpret_cast< ::complex*
>(dat), n1, stype, dtype, mtype);
213 static void destroy(SuperMatrix* )
218 struct SuperLUSolveChooser<
std::complex<float> >
220 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
221 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
222 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
223 float *rpg,
float *rcond,
float *ferr,
float *berr,
224 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
226#if SUPERLU_MIN_VERSION_5
228 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
229 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
230 &gLU, memusage, stat, info);
232 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
233 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
234 memusage, stat, info);
240 struct QuerySpaceChooser<
std::complex<float> >
242 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
244 cQuerySpace(L,U,memusage);
262 template<
typename T,
typename A,
int n,
int m>
265 BlockVector<FieldVector<T,m>,
266 typename A::template rebind<FieldVector<T,m> >::other>,
267 BlockVector<FieldVector<T,n>,
268 typename A::template rebind<FieldVector<T,n> >::other> >
281 typename A::template rebind<FieldVector<T,m> >::other>
domain_type;
285 typename A::template rebind<FieldVector<T,n> >::other>
range_type;
290 return SolverCategory::Category::sequential;
307 explicit SuperLU(
const Matrix& mat,
bool verbose=
false,
308 bool reusevector=
true);
336 void apply(T* x, T* b);
339 void setMatrix(
const Matrix& mat);
341 typename SuperLUMatrix::size_type nnz()
const
347 void setSubMatrix(
const Matrix& mat,
const S& rowIndexSet);
349 void setVerbosity(
bool v);
357 const char* name() {
return "SuperLU"; }
359 template<
class M,
class X,
class TM,
class TD,
class T1>
360 friend class SeqOverlappingSchwarz;
361 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
363 SuperLUMatrix& getInternalMatrix() {
return mat; }
369 SuperMatrix L, U, B, X;
370 int *perm_c, *perm_r, *etree;
371 typename GetSuperLUType<T>::float_type *R, *C;
373 superlu_options_t options;
377 bool first, verbose, reusevector;
380 template<
typename T,
typename A,
int n,
int m>
381 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
384 if(mat.
N()+mat.
M()>0)
388 template<
typename T,
typename A,
int n,
int m>
389 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
397 Destroy_SuperNode_Matrix(&L);
398 Destroy_CompCol_Matrix(&U);
401 if(!first && reusevector) {
402 SUPERLU_FREE(B.Store);
403 SUPERLU_FREE(X.Store);
408 template<
typename T,
typename A,
int n,
int m>
409 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
410 ::SuperLU(
const Matrix& mat_,
bool verbose_,
bool reusevector_)
411 : work(0), lwork(0), first(true), verbose(verbose_),
412 reusevector(reusevector_)
417 template<
typename T,
typename A,
int n,
int m>
418 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
419 : work(0), lwork(0),verbose(false),
422 template<
typename T,
typename A,
int n,
int m>
423 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(
bool v)
428 template<
typename T,
typename A,
int n,
int m>
429 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(
const Matrix& mat_)
431 if(mat.
N()+mat.
M()>0) {
441 template<
typename T,
typename A,
int n,
int m>
443 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(
const Matrix& mat_,
446 if(mat.
N()+mat.
M()>0) {
452 mat.setMatrix(mat_,mrs);
456 template<
typename T,
typename A,
int n,
int m>
457 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
461 perm_c =
new int[mat.
M()];
462 perm_r =
new int[mat.
N()];
463 etree =
new int[mat.
M()];
464 R =
new typename GetSuperLUType<T>::float_type[mat.
N()];
465 C =
new typename GetSuperLUType<T>::float_type[mat.
M()];
467 set_default_options(&options);
471 B.Dtype=GetSuperLUType<T>::type;
474 fakeFormat.lda=mat.
N();
477 X.Dtype=GetSuperLUType<T>::type;
482 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
484 mem_usage_t memusage;
488 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
489 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
490 &berr, &memusage, &stat, &info);
493 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
495 if ( info == 0 || info == n+1 ) {
497 if ( options.PivotGrowth )
498 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
499 if ( options.ConditionNumber )
500 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
501 SCformat* Lstore = (SCformat *) L.Store;
502 NCformat* Ustore = (NCformat *) U.Store;
503 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
504 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
505 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
506 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
507 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
509 std::cout<<stat.expansions<<std::endl;
511 }
else if ( info > 0 && lwork == -1 ) {
512 dinfo<<
"** Estimated memory: "<< info - n<<std::endl;
514 if ( options.PrintStat ) StatPrint(&stat);
542 options.Fact = FACTORED;
545 template<
typename T,
typename A,
int n,
int m>
546 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
549 if (mat.
N() != b.dim())
550 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
551 if (mat.
M() != x.dim())
552 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
553 if (mat.
M()+mat.
N()==0)
556 SuperMatrix* mB = &B;
557 SuperMatrix* mX = &X;
561 SuperLUDenseMatChooser<T>::create(&B, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
562 SuperLUDenseMatChooser<T>::create(&X, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
565 ((DNformat*)B.Store)->nzval=&b[0];
566 ((DNformat*)X.Store)->nzval=&x[0];
569 SuperLUDenseMatChooser<T>::create(&rB, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
570 SuperLUDenseMatChooser<T>::create(&rX, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
574 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
576 mem_usage_t memusage;
586#ifdef SUPERLU_MIN_VERSION_4_3
587 options.IterRefine=SLU_DOUBLE;
589 options.IterRefine=DOUBLE;
592 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
593 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
594 &memusage, &stat, &info);
614 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
616 if ( info == 0 || info == n+1 ) {
618 if ( options.IterRefine ) {
619 std::cout<<
"Iterative Refinement: steps="
620 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
622 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
623 }
else if ( info > 0 && lwork == -1 ) {
624 std::cout<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
627 if ( options.PrintStat ) StatPrint(&stat);
631 SUPERLU_FREE(rB.Store);
632 SUPERLU_FREE(rX.Store);
636 template<
typename T,
typename A,
int n,
int m>
637 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
640 if(mat.
N()+mat.
M()==0)
643 SuperMatrix* mB = &B;
644 SuperMatrix* mX = &X;
648 SuperLUDenseMatChooser<T>::create(&B, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
649 SuperLUDenseMatChooser<T>::create(&X, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
652 ((DNformat*) B.Store)->nzval=b;
653 ((DNformat*)X.Store)->nzval=x;
656 SuperLUDenseMatChooser<T>::create(&rB, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
657 SuperLUDenseMatChooser<T>::create(&rX, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
662 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
664 mem_usage_t memusage;
669#ifdef SUPERLU_MIN_VERSION_4_3
670 options.IterRefine=SLU_DOUBLE;
672 options.IterRefine=DOUBLE;
675 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
676 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
677 &memusage, &stat, &info);
680 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
682 if ( info == 0 || info == n+1 ) {
684 if ( options.IterRefine ) {
685 dinfo<<
"Iterative Refinement: steps="
686 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
688 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
689 }
else if ( info > 0 && lwork == -1 ) {
690 dinfo<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
692 if ( options.PrintStat ) StatPrint(&stat);
697 SUPERLU_FREE(rB.Store);
698 SUPERLU_FREE(rX.Store);
703 template<
typename T,
typename A,
int n,
int m>
709 template<
typename T,
typename A,
int n,
int m>
710 struct StoresColumnCompressed<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
712 enum { value =
true };
717#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:423
A vector of blocks with memory management.
Definition: bvector.hh:317
A dense n x m matrix.
Definition: fmatrix.hh:68
vector space out of a tensor product of fields.
Definition: fvector.hh:93
derive error class from the base class in common
Definition: istlexception.hh:16
Abstract base class for all solvers.
Definition: solver.hh:91
A generic dynamic dense matrix.
Definition: matrix.hh:555
size_type M() const
Return the number of columns.
Definition: matrix.hh:695
size_type N() const
Return the number of rows.
Definition: matrix.hh:690
Dune::BlockVector< FieldVector< T, n >, typename A::template rebind< FieldVector< T, n > >::other > range_type
The type of the range of the solver.
Definition: superlu.hh:285
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:288
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:275
Dune::BlockVector< FieldVector< T, m >, typename A::template rebind< FieldVector< T, m > >::other > domain_type
The type of the domain of the solver.
Definition: superlu.hh:281
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:327
SuperMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A > > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:277
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.
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:58
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune namespace.
Definition: alignedallocator.hh:10
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:41
int iterations
Number of iterations.
Definition: solver.hh:59
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
Category
Definition: solvercategory.hh:21