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>
34 template<
class Matrix>
38 template<
class M,
class T,
class TM,
class TD,
class TA>
41 template<
class T,
bool tag>
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;
308 bool reusevector=
true);
329 DUNE_UNUSED_PARAMETER(reduction);
336 void apply(T* x, T* b);
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>
369 SuperMatrix L, U, B, X;
370 int *perm_c, *perm_r, *etree;
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> >
388 template<
typename T,
typename A,
int n,
int m>
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>
411 : work(0), lwork(0), first(true), verbose(verbose_),
412 reusevector(reusevector_)
417 template<
typename T,
typename A,
int n,
int m>
419 : work(0), lwork(0),verbose(false),
422 template<
typename T,
typename A,
int n,
int m>
428 template<
typename T,
typename A,
int n,
int m>
441 template<
typename T,
typename A,
int n,
int m>
452 mat.setMatrix(mat_,mrs);
456 template<
typename T,
typename A,
int n,
int m>
461 perm_c =
new int[
mat.
M()];
462 perm_r =
new int[
mat.
N()];
463 etree =
new int[
mat.
M()];
467 set_default_options(&options);
474 fakeFormat.lda=
mat.
N();
484 mem_usage_t memusage;
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;
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!");
554 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
556 SuperMatrix* mB = &B;
557 SuperMatrix* mX = &X;
565 ((DNformat*)B.Store)->nzval=&b[0];
566 ((DNformat*)X.Store)->nzval=&x[0];
576 mem_usage_t memusage;
586 #ifdef SUPERLU_MIN_VERSION_4_3
587 options.IterRefine=SLU_DOUBLE;
589 options.IterRefine=DOUBLE;
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>
641 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
643 SuperMatrix* mB = &B;
644 SuperMatrix* mX = &X;
652 ((DNformat*) B.Store)->nzval=b;
653 ((DNformat*)X.Store)->nzval=x;
664 mem_usage_t memusage;
669 #ifdef SUPERLU_MIN_VERSION_4_3
670 options.IterRefine=SLU_DOUBLE;
672 options.IterRefine=DOUBLE;
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>
717 #undef FIRSTCOL_OF_SNODE
722 #undef SUPERLU_MALLOC
742 #undef DROP_SECONDARY
747 #endif // HAVE_SUPERLU
748 #endif // DUNE_SUPERLU_HH