3 #ifndef DUNE_ISTL_SUPERLU_HH 4 #define DUNE_ISTL_SUPERLU_HH 14 #define SUPERLU_NTYPE 1 18 #include "slu_sdefs.h" 22 #include "slu_ddefs.h" 26 #include "slu_cdefs.h" 30 #include "slu_zdefs.h" 40 #include <dune/common/fmatrix.hh> 41 #include <dune/common/fvector.hh> 42 #include <dune/common/stdstreams.hh> 58 template<
class Matrix>
62 template<
class M,
class T,
class TM,
class TD,
class TA>
65 template<
class T,
bool tag>
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*)
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);
124 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
126 sQuerySpace(L,U,memusage);
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);
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);
172 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
174 dQuerySpace(L,U,memusage);
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*)
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);
219 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
221 zQuerySpace(L,U,memusage);
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* )
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);
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);
346 DUNE_UNUSED_PARAMETER(reduction);
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>
381 SuperLUMatrix& getInternalMatrix() {
return mat; }
387 SuperMatrix L, U, B, X;
388 int *perm_c, *perm_r, *etree;
391 superlu_options_t options;
395 bool first, verbose, reusevector;
398 template<
typename T,
typename A,
int n,
int m>
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> >
429 : work(0), lwork(0), first(true), verbose(verbose_),
430 reusevector(reusevector_)
435 template<
typename T,
typename A,
int n,
int m>
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_)
459 template<
typename T,
typename A,
int n,
int m>
461 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(
const Matrix& mat_,
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()];
485 set_default_options(&options);
492 fakeFormat.lda=
mat.N();
502 mem_usage_t memusage;
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;
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> >
568 DUNE_THROW(
ISTLError,
"Size of right-hand-side vector b does not match the number of matrix rows!");
570 DUNE_THROW(
ISTLError,
"Size of solution vector x does not match the number of matrix columns!");
572 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
574 SuperMatrix* mB = &B;
575 SuperMatrix* mX = &X;
583 ((DNformat*)B.Store)->nzval=&b[0];
584 ((DNformat*)X.Store)->nzval=&x[0];
594 mem_usage_t memusage;
604 #ifdef SUPERLU_MIN_VERSION_4_3 605 options.IterRefine=SLU_DOUBLE;
607 options.IterRefine=DOUBLE;
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> >
659 DUNE_THROW(
ISTLError,
"Matrix of SuperLU is null!");
661 SuperMatrix* mB = &B;
662 SuperMatrix* mX = &X;
670 ((DNformat*) B.Store)->nzval=b;
671 ((DNformat*)X.Store)->nzval=x;
682 mem_usage_t memusage;
687 #ifdef SUPERLU_MIN_VERSION_4_3 688 options.IterRefine=SLU_DOUBLE;
690 options.IterRefine=DOUBLE;
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>
730 enum { value =
true };
735 #undef FIRSTCOL_OF_SNODE 740 #undef SUPERLU_MALLOC 760 #undef DROP_SECONDARY 765 #endif // HAVE_SUPERLU 766 #endif // DUNE_SUPERLU_HH Definition: superlu.hh:69
static void create(SuperMatrix *mat, int n, int m, double *dat, int n1, Stype_t stype, Dtype_t dtype, Mtype_t mtype)
Definition: superlu.hh:137
static void querySpace(SuperMatrix *L, SuperMatrix *U, mem_usage_t *memusage)
Definition: superlu.hh:172
Definition: solvertype.hh:27
A vector of blocks with memory management.
Definition: bvector.hh:312
Definition: superlu.hh:77
Definition: superlu.hh:59
const char * name()
Definition: superlu.hh:374
Implementations of the inverse operator interface.
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:344
Templates characterizing the type of a solver.
Definition: colcompmatrix.hh:160
static void destroy(SuperMatrix *)
Definition: superlu.hh:144
Definition: supermatrix.hh:147
Abstract base class for all solvers.
Definition: solver.hh:79
Definition: matrixutils.hh:25
size_type dim() const
dimension of the vector space
Definition: bvector.hh:283
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: basearray.hh:19
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:299
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Definition: supermatrix.hh:194
Definition: solvertype.hh:13
Implementation of the BCRSMatrix class.
Definition: superlu.hh:73
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
int iterations
Number of iterations.
Definition: solver.hh:50
Definition: supermatrix.hh:305
static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree, char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix *B, SuperMatrix *X, double *rpg, double *rcond, double *ferr, double *berr, mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
Definition: superlu.hh:150
Definition: superlu.hh:81
derive error class from the base class in common
Definition: istlexception.hh:16
Sequential overlapping Schwarz preconditioner.
Definition: colcompmatrix.hh:157
SuperLUMatrix::size_type nnz() const
Definition: superlu.hh:358
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
Statistics about the application of an inverse operator.
Definition: solver.hh:31