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*)
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);
100 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
102 sQuerySpace(L,U,memusage);
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 * )
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);
148 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
150 dQuerySpace(L,U,memusage);
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*)
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);
195 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
197 zQuerySpace(L,U,memusage);
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* )
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);
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);
329 DUNE_UNUSED_PARAMETER(reduction);
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>
363 SuperLUMatrix& getInternalMatrix() {
return mat; }
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>
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> >
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>
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_)
441 template<
typename T,
typename A,
int n,
int m>
443 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(
const Matrix& mat_,
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()];
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>
637 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
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>
712 enum { value =
true };
717 #undef FIRSTCOL_OF_SNODE 722 #undef SUPERLU_MALLOC 742 #undef DROP_SECONDARY 747 #endif // HAVE_SUPERLU 748 #endif // DUNE_SUPERLU_HH Definition: allocator.hh:7
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:422
SuperLUMatrix::size_type nnz() const
Definition: superlu.hh:341
This file implements a vector space as a tensor product of a given vector space. The number of compon...
const char * name()
Definition: superlu.hh:357
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:288
Implementations of the inverse operator interface.
Definition: colcompmatrix.hh:160
Definition: solvertype.hh:27
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: matrixutils.hh:25
Implementation of the BCRSMatrix class.
Definition: superlu.hh:53
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
Sequential overlapping Schwarz preconditioner.
Definition: colcompmatrix.hh:157
Definition: superlu.hh:35
derive error class from the base class in common
Definition: istlexception.hh:16
Definition: supermatrix.hh:176
bool converged
True if convergence criterion has been met.
Definition: solver.hh:65
Abstract base class for all solvers.
Definition: solver.hh:91
Definition: supermatrix.hh:129
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
Category
Definition: solvercategory.hh:21
Definition: supermatrix.hh:287
Definition: solvertype.hh:13
A vector of blocks with memory management.
Definition: bvector.hh:316
Definition: superlu.hh:49
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:275
Definition: superlu.hh:57
Statistics about the application of an inverse operator.
Definition: solver.hh:40
Definition: superlu.hh:45
int iterations
Number of iterations.
Definition: solver.hh:59
Templates characterizing the type of a solver.
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:327