3#ifndef DUNE_ISTL_SUPERLU_HH
4#define DUNE_ISTL_SUPERLU_HH
14#define SUPERLU_NTYPE 1
16#ifdef SUPERLU_POST_2005_VERSION
42#warning Support for SuperLU older than SuperLU 3.0 from August 2005 is deprecated.
55#include "supermatrix.hh"
60#include "istlexception.hh"
79 template<
class Matrix>
83 template<
class M,
class T,
class TM,
class TD,
class TA>
84 class SeqOverlappingSchwarz;
86 template<
class T,
bool tag>
87 struct SeqOverlappingSchwarzAssemblerHelper;
90 struct SuperLUSolveChooser
94 struct SuperLUDenseMatChooser
98 struct SuperLUQueryChooser
102 struct QuerySpaceChooser
107 struct SuperLUDenseMatChooser<float>
109 static void create(SuperMatrix *mat,
int n,
int m,
float *dat,
int n1,
110 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
112 sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
116 static void destroy(SuperMatrix*)
121 struct SuperLUSolveChooser<float>
123 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
124 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
125 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
126 float *rpg,
float *rcond,
float *ferr,
float *berr,
127 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
129 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
130 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
131 memusage, stat, info);
136 struct QuerySpaceChooser<float>
138 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
140 sQuerySpace(L,U,memusage);
149 struct SuperLUDenseMatChooser<double>
151 static void create(SuperMatrix *mat,
int n,
int m,
double *dat,
int n1,
152 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
154 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
158 static void destroy(SuperMatrix * )
162 struct SuperLUSolveChooser<double>
164 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
165 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
166 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
167 double *rpg,
double *rcond,
double *ferr,
double *berr,
168 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
170 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
171 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
172 memusage, stat, info);
177 struct QuerySpaceChooser<double>
179 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
181 dQuerySpace(L,U,memusage);
188 struct SuperLUDenseMatChooser<
std::complex<double> >
190 static void create(SuperMatrix *mat,
int n,
int m, std::complex<double> *dat,
int n1,
191 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
193 zCreate_Dense_Matrix(mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
197 static void destroy(SuperMatrix*)
202 struct SuperLUSolveChooser<
std::complex<double> >
204 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
205 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
206 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
207 double *rpg,
double *rcond,
double *ferr,
double *berr,
208 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
210 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
211 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
212 memusage, stat, info);
217 struct QuerySpaceChooser<
std::complex<double> >
219 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
221 zQuerySpace(L,U,memusage);
228 struct SuperLUDenseMatChooser<
std::complex<float> >
230 static void create(SuperMatrix *mat,
int n,
int m, std::complex<float> *dat,
int n1,
231 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
233 cCreate_Dense_Matrix(mat, n, m,
reinterpret_cast< ::complex*
>(dat), n1, stype, dtype, mtype);
237 static void destroy(SuperMatrix* )
242 struct SuperLUSolveChooser<
std::complex<float> >
244 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
245 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
246 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
247 float *rpg,
float *rcond,
float *ferr,
float *berr,
248 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
250 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
251 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
252 memusage, stat, info);
257 struct QuerySpaceChooser<
std::complex<float> >
259 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
261 cQuerySpace(L,U,memusage);
279 template<
typename T,
typename A,
int n,
int m>
282 BlockVector<FieldVector<T,m>,
283 typename A::template rebind<FieldVector<T,m> >::other>,
284 BlockVector<FieldVector<T,n>,
285 typename A::template rebind<FieldVector<T,n> >::other> >
298 typename A::template rebind<FieldVector<T,m> >::other>
domain_type;
302 typename A::template rebind<FieldVector<T,n> >::other>
range_type;
317 explicit SuperLU(
const Matrix& mat,
bool verbose=
false,
318 bool reusevector=
true);
346 void apply(T* x, T* b);
349 void setMatrix(
const Matrix& mat);
351 typename SuperLUMatrix::size_type nnz()
const
357 void setSubMatrix(
const Matrix& mat,
const S& rowIndexSet);
359 void setVerbosity(
bool v);
367 const char* name() {
return "SuperLU"; }
369 friend class std::mem_fun_ref_t<void,SuperLU>;
370 template<
class M,
class X,
class TM,
class TD,
class T1>
371 friend class SeqOverlappingSchwarz;
372 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
374 SuperLUMatrix& getInternalMatrix() {
return mat; }
380 SuperMatrix L, U, B, X;
381 int *perm_c, *perm_r, *etree;
382 typename GetSuperLUType<T>::float_type *R, *C;
384 superlu_options_t options;
388 bool first, verbose, reusevector;
391 template<
typename T,
typename A,
int n,
int m>
392 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
395 if(mat.
N()+mat.
M()>0)
399 template<
typename T,
typename A,
int n,
int m>
400 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
408 Destroy_SuperNode_Matrix(&L);
409 Destroy_CompCol_Matrix(&U);
412 if(!first && reusevector) {
413 SUPERLU_FREE(B.Store);
414 SUPERLU_FREE(X.Store);
419 template<
typename T,
typename A,
int n,
int m>
420 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
421 ::SuperLU(
const Matrix& mat_,
bool verbose_,
bool reusevector_)
422 : work(0), lwork(0), first(true), verbose(verbose_),
423 reusevector(reusevector_)
428 template<
typename T,
typename A,
int n,
int m>
429 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
430 : work(0), lwork(0),verbose(false),
433 template<
typename T,
typename A,
int n,
int m>
434 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(
bool v)
439 template<
typename T,
typename A,
int n,
int m>
440 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(
const Matrix& mat_)
442 if(mat.
N()+mat.
M()>0) {
452 template<
typename T,
typename A,
int n,
int m>
454 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(
const Matrix& mat_,
457 if(mat.
N()+mat.
M()>0) {
463 mat.setMatrix(mat_,mrs);
467 template<
typename T,
typename A,
int n,
int m>
468 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
472 perm_c =
new int[mat.
M()];
473 perm_r =
new int[mat.
N()];
474 etree =
new int[mat.
M()];
475 R =
new typename GetSuperLUType<T>::float_type[mat.
N()];
476 C =
new typename GetSuperLUType<T>::float_type[mat.
M()];
478 set_default_options(&options);
482 B.Dtype=GetSuperLUType<T>::type;
485 fakeFormat.lda=mat.
N();
488 X.Dtype=GetSuperLUType<T>::type;
493 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
495 mem_usage_t memusage;
499 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
500 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
501 &berr, &memusage, &stat, &info);
504 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
506 if ( info == 0 || info == n+1 ) {
508 if ( options.PivotGrowth )
509 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
510 if ( options.ConditionNumber )
511 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
512 SCformat* Lstore = (SCformat *) L.Store;
513 NCformat* Ustore = (NCformat *) U.Store;
514 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
515 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
516 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
517 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
518 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
521#ifdef HAVE_MEM_USAGE_T_EXPANSIONS
522 std::cout<<memusage.expansions<<std::endl;
524 std::cout<<stat.expansions<<std::endl;
526 }
else if ( info > 0 && lwork == -1 ) {
527 dinfo<<
"** Estimated memory: "<< info - n<<std::endl;
529 if ( options.PrintStat ) StatPrint(&stat);
557 options.Fact = FACTORED;
560 template<
typename T,
typename A,
int n,
int m>
561 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
564 if (mat.
N() != b.
dim())
565 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
566 if (mat.
M() != x.
dim())
567 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
568 if (mat.
M()+mat.
N()==0)
571 SuperMatrix* mB = &B;
572 SuperMatrix* mX = &X;
576 SuperLUDenseMatChooser<T>::create(&B, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
577 SuperLUDenseMatChooser<T>::create(&X, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
580 ((DNformat*)B.Store)->nzval=&b[0];
581 ((DNformat*)X.Store)->nzval=&x[0];
584 SuperLUDenseMatChooser<T>::create(&rB, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
585 SuperLUDenseMatChooser<T>::create(&rX, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
589 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
591 mem_usage_t memusage;
601#ifdef SUPERLU_MIN_VERSION_4_3
602 options.IterRefine=SLU_DOUBLE;
604 options.IterRefine=DOUBLE;
607 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
608 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
609 &memusage, &stat, &info);
629 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
631 if ( info == 0 || info == n+1 ) {
633 if ( options.IterRefine ) {
634 std::cout<<
"Iterative Refinement: steps="
635 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
637 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
638 }
else if ( info > 0 && lwork == -1 ) {
639 std::cout<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
642 if ( options.PrintStat ) StatPrint(&stat);
646 SUPERLU_FREE(rB.Store);
647 SUPERLU_FREE(rX.Store);
651 template<
typename T,
typename A,
int n,
int m>
652 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
655 if(mat.
N()+mat.
M()==0)
658 SuperMatrix* mB = &B;
659 SuperMatrix* mX = &X;
663 SuperLUDenseMatChooser<T>::create(&B, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
664 SuperLUDenseMatChooser<T>::create(&X, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
667 ((DNformat*) B.Store)->nzval=b;
668 ((DNformat*)X.Store)->nzval=x;
671 SuperLUDenseMatChooser<T>::create(&rB, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
672 SuperLUDenseMatChooser<T>::create(&rX, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
677 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
679 mem_usage_t memusage;
684#ifdef SUPERLU_MIN_VERSION_4_3
685 options.IterRefine=SLU_DOUBLE;
687 options.IterRefine=DOUBLE;
690 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
691 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
692 &memusage, &stat, &info);
695 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
697 if ( info == 0 || info == n+1 ) {
699 if ( options.IterRefine ) {
700 dinfo<<
"Iterative Refinement: steps="
701 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
703 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
704 }
else if ( info > 0 && lwork == -1 ) {
705 dinfo<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
707 if ( options.PrintStat ) StatPrint(&stat);
712 SUPERLU_FREE(rB.Store);
713 SUPERLU_FREE(rX.Store);
718 template<
typename T,
typename A,
int n,
int m>
724 template<
typename T,
typename A,
int n,
int m>
725 struct StoresColumnCompressed<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
727 enum { value =
true };
732#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:413
A vector of blocks with memory management.
Definition: bvector.hh:253
A dense n x m matrix.
Definition: fmatrix.hh:67
vector space out of a tensor product of fields.
Definition: fvector.hh:94
derive error class from the base class in common
Definition: istlexception.hh:16
Abstract base class for all solvers.
Definition: solver.hh:79
A generic dynamic dense matrix.
Definition: matrix.hh:25
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
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:302
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:292
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:298
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:337
SuperMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A > > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:294
size_type dim() const
dimension of the vector space
Definition: bvector.hh:223
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:243
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune namespace.
Definition: alignment.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:32
int iterations
Number of iterations.
Definition: solver.hh:50
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18