3#ifndef DUNE_ISTL_SUPERLU_HH
4#define DUNE_ISTL_SUPERLU_HH
14#define SUPERLU_NTYPE 1
34#include "supermatrix.hh"
39#include "istlexception.hh"
58 template<
class Matrix>
62 template<
class M,
class T,
class TM,
class TD,
class TA>
63 class SeqOverlappingSchwarz;
65 template<
class T,
bool tag>
66 struct SeqOverlappingSchwarzAssemblerHelper;
69 struct SuperLUSolveChooser
73 struct SuperLUDenseMatChooser
77 struct SuperLUQueryChooser
81 struct QuerySpaceChooser
86 struct SuperLUDenseMatChooser<float>
88 static void create(SuperMatrix *mat,
int n,
int m,
float *dat,
int n1,
89 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
91 sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
95 static void destroy(SuperMatrix*)
100 struct SuperLUSolveChooser<float>
102 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
103 char *equed,
float *R,
float *C, SuperMatrix *L, SuperMatrix *U,
104 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
105 float *rpg,
float *rcond,
float *ferr,
float *berr,
106 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
108#if SUPERLU_MIN_VERSION_5
110 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
111 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
112 &gLU, memusage, stat, info);
114 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
115 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
116 memusage, stat, info);
122 struct QuerySpaceChooser<float>
124 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
126 sQuerySpace(L,U,memusage);
135 struct SuperLUDenseMatChooser<double>
137 static void create(SuperMatrix *mat,
int n,
int m,
double *dat,
int n1,
138 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
140 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
144 static void destroy(SuperMatrix * )
148 struct SuperLUSolveChooser<double>
150 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
151 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
152 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
153 double *rpg,
double *rcond,
double *ferr,
double *berr,
154 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
156#if SUPERLU_MIN_VERSION_5
158 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
159 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
160 &gLU, memusage, stat, info);
162 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
163 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
164 memusage, stat, info);
170 struct QuerySpaceChooser<double>
172 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
174 dQuerySpace(L,U,memusage);
181 struct SuperLUDenseMatChooser<
std::complex<double> >
183 static void create(SuperMatrix *mat,
int n,
int m, std::complex<double> *dat,
int n1,
184 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
186 zCreate_Dense_Matrix(mat, n, m,
reinterpret_cast<doublecomplex*
>(dat), n1, stype, dtype, mtype);
190 static void destroy(SuperMatrix*)
195 struct SuperLUSolveChooser<
std::complex<double> >
197 static void solve(superlu_options_t *options, SuperMatrix *mat,
int *perm_c,
int *perm_r,
int *etree,
198 char *equed,
double *R,
double *C, SuperMatrix *L, SuperMatrix *U,
199 void *work,
int lwork, SuperMatrix *B, SuperMatrix *X,
200 double *rpg,
double *rcond,
double *ferr,
double *berr,
201 mem_usage_t *memusage, SuperLUStat_t *stat,
int *info)
203#if SUPERLU_MIN_VERSION_5
205 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
206 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
207 &gLU, memusage, stat, info);
209 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
210 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
211 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#if SUPERLU_MIN_VERSION_5
252 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
253 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
254 &gLU, memusage, stat, info);
256 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
257 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
258 memusage, stat, info);
264 struct QuerySpaceChooser<
std::complex<float> >
266 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
268 cQuerySpace(L,U,memusage);
286 template<
typename T,
typename A,
int n,
int m>
289 BlockVector<FieldVector<T,m>,
290 typename A::template rebind<FieldVector<T,m> >::other>,
291 BlockVector<FieldVector<T,n>,
292 typename A::template rebind<FieldVector<T,n> >::other> >
305 typename A::template rebind<FieldVector<T,m> >::other>
domain_type;
309 typename A::template rebind<FieldVector<T,n> >::other>
range_type;
324 explicit SuperLU(
const Matrix& mat,
bool verbose=
false,
325 bool reusevector=
true);
353 void apply(T* x, T* b);
356 void setMatrix(
const Matrix& mat);
358 typename SuperLUMatrix::size_type nnz()
const
364 void setSubMatrix(
const Matrix& mat,
const S& rowIndexSet);
366 void setVerbosity(
bool v);
374 const char* name() {
return "SuperLU"; }
376 friend class std::mem_fun_ref_t<void,SuperLU>;
377 template<
class M,
class X,
class TM,
class TD,
class T1>
378 friend class SeqOverlappingSchwarz;
379 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
381 SuperLUMatrix& getInternalMatrix() {
return mat; }
387 SuperMatrix L, U, B, X;
388 int *perm_c, *perm_r, *etree;
389 typename GetSuperLUType<T>::float_type *R, *C;
391 superlu_options_t options;
395 bool first, verbose, reusevector;
398 template<
typename T,
typename A,
int n,
int m>
399 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
402 if(mat.
N()+mat.
M()>0)
406 template<
typename T,
typename A,
int n,
int m>
407 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
415 Destroy_SuperNode_Matrix(&L);
416 Destroy_CompCol_Matrix(&U);
419 if(!first && reusevector) {
420 SUPERLU_FREE(B.Store);
421 SUPERLU_FREE(X.Store);
426 template<
typename T,
typename A,
int n,
int m>
427 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
428 ::SuperLU(
const Matrix& mat_,
bool verbose_,
bool reusevector_)
429 : work(0), lwork(0), first(true), verbose(verbose_),
430 reusevector(reusevector_)
435 template<
typename T,
typename A,
int n,
int m>
436 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
437 : work(0), lwork(0),verbose(false),
440 template<
typename T,
typename A,
int n,
int m>
441 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(
bool v)
446 template<
typename T,
typename A,
int n,
int m>
447 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(
const Matrix& mat_)
449 if(mat.
N()+mat.
M()>0) {
459 template<
typename T,
typename A,
int n,
int m>
461 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(
const Matrix& mat_,
464 if(mat.
N()+mat.
M()>0) {
470 mat.setMatrix(mat_,mrs);
474 template<
typename T,
typename A,
int n,
int m>
475 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
479 perm_c =
new int[mat.
M()];
480 perm_r =
new int[mat.
N()];
481 etree =
new int[mat.
M()];
482 R =
new typename GetSuperLUType<T>::float_type[mat.
N()];
483 C =
new typename GetSuperLUType<T>::float_type[mat.
M()];
485 set_default_options(&options);
489 B.Dtype=GetSuperLUType<T>::type;
492 fakeFormat.lda=mat.
N();
495 X.Dtype=GetSuperLUType<T>::type;
500 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
502 mem_usage_t memusage;
506 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
507 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
508 &berr, &memusage, &stat, &info);
511 dinfo<<
"LU factorization: dgssvx() returns info "<< info<<std::endl;
513 if ( info == 0 || info == n+1 ) {
515 if ( options.PivotGrowth )
516 dinfo<<
"Recip. pivot growth = "<<rpg<<std::endl;
517 if ( options.ConditionNumber )
518 dinfo<<
"Recip. condition number = %e\n"<< rcond<<std::endl;
519 SCformat* Lstore = (SCformat *) L.Store;
520 NCformat* Ustore = (NCformat *) U.Store;
521 dinfo<<
"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
522 dinfo<<
"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
523 dinfo<<
"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
524 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
525 dinfo<<
"L\\U MB "<<memusage.for_lu/1e6<<
" \ttotal MB needed "<<memusage.total_needed/1e6
527 std::cout<<stat.expansions<<std::endl;
529 }
else if ( info > 0 && lwork == -1 ) {
530 dinfo<<
"** Estimated memory: "<< info - n<<std::endl;
532 if ( options.PrintStat ) StatPrint(&stat);
560 options.Fact = FACTORED;
563 template<
typename T,
typename A,
int n,
int m>
564 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
567 if (mat.
N() != b.
dim())
568 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
569 if (mat.
M() != x.
dim())
570 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
571 if (mat.
M()+mat.
N()==0)
574 SuperMatrix* mB = &B;
575 SuperMatrix* mX = &X;
579 SuperLUDenseMatChooser<T>::create(&B, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
580 SuperLUDenseMatChooser<T>::create(&X, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
583 ((DNformat*)B.Store)->nzval=&b[0];
584 ((DNformat*)X.Store)->nzval=&x[0];
587 SuperLUDenseMatChooser<T>::create(&rB, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&b[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
588 SuperLUDenseMatChooser<T>::create(&rX, (
int)mat.
N(), 1,
reinterpret_cast<T*
>(&x[0]), (
int)mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
592 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
594 mem_usage_t memusage;
604#ifdef SUPERLU_MIN_VERSION_4_3
605 options.IterRefine=SLU_DOUBLE;
607 options.IterRefine=DOUBLE;
610 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
611 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
612 &memusage, &stat, &info);
632 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
634 if ( info == 0 || info == n+1 ) {
636 if ( options.IterRefine ) {
637 std::cout<<
"Iterative Refinement: steps="
638 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
640 std::cout<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
641 }
else if ( info > 0 && lwork == -1 ) {
642 std::cout<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
645 if ( options.PrintStat ) StatPrint(&stat);
649 SUPERLU_FREE(rB.Store);
650 SUPERLU_FREE(rX.Store);
654 template<
typename T,
typename A,
int n,
int m>
655 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
658 if(mat.
N()+mat.
M()==0)
661 SuperMatrix* mB = &B;
662 SuperMatrix* mX = &X;
666 SuperLUDenseMatChooser<T>::create(&B, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
667 SuperLUDenseMatChooser<T>::create(&X, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
670 ((DNformat*) B.Store)->nzval=b;
671 ((DNformat*)X.Store)->nzval=x;
674 SuperLUDenseMatChooser<T>::create(&rB, mat.
N(), 1, b, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
675 SuperLUDenseMatChooser<T>::create(&rX, mat.
N(), 1, x, mat.
N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
680 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
682 mem_usage_t memusage;
687#ifdef SUPERLU_MIN_VERSION_4_3
688 options.IterRefine=SLU_DOUBLE;
690 options.IterRefine=DOUBLE;
693 SuperLUSolveChooser<T>::solve(&options, &
static_cast<SuperMatrix&
>(mat), perm_c, perm_r, etree, &equed, R, C,
694 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
695 &memusage, &stat, &info);
698 dinfo<<
"Triangular solve: dgssvx() returns info "<< info<<std::endl;
700 if ( info == 0 || info == n+1 ) {
702 if ( options.IterRefine ) {
703 dinfo<<
"Iterative Refinement: steps="
704 <<stat.RefineSteps<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
706 dinfo<<
" FERR="<<ferr<<
" BERR="<<berr<<std::endl;
707 }
else if ( info > 0 && lwork == -1 ) {
708 dinfo<<
"** Estimated memory: "<< info - n<<
" bytes"<<std::endl;
710 if ( options.PrintStat ) StatPrint(&stat);
715 SUPERLU_FREE(rB.Store);
716 SUPERLU_FREE(rX.Store);
721 template<
typename T,
typename A,
int n,
int m>
727 template<
typename T,
typename A,
int n,
int m>
728 struct StoresColumnCompressed<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
730 enum { value =
true };
735#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:313
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:79
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:309
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:299
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:305
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:344
SuperMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A > > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:301
size_type dim() const
dimension of the vector space
Definition: bvector.hh:283
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:216
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune namespace.
Definition: alignment.hh:11
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 intentionally unused function parameters with.
Definition: unused.hh:18