3#ifndef DUNE_ISTL_SUPERLU_HH
4#define DUNE_ISTL_SUPERLU_HH
14#define SUPERLU_NTYPE 1
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>
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,
114 sgssvx(options,
mat, perm_c, perm_r, etree, equed, R, C,
115 L, U, work, lwork, B, X,
rpg,
rcond,
ferr,
berr,
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> >
325 bool reusevector=
true);
353 void apply(T* x, T* b);
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"; }
377 template<
class M,
class X,
class TM,
class TD,
class T1>
388 int *perm_c, *perm_r, *etree;
395 bool first, verbose, reusevector;
398 template<
typename T,
typename A,
int n,
int m>
399 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
406 template<
typename T,
typename A,
int n,
int m>
419 if(!first && reusevector) {
426 template<
typename T,
typename A,
int n,
int m>
429 : work(0), lwork(0), first(
true), verbose(verbose_),
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>
446 template<
typename T,
typename A,
int n,
int m>
459 template<
typename T,
typename A,
int n,
int m>
474 template<
typename T,
typename A,
int n,
int m>
479 perm_c =
new int[
mat.
M()];
480 perm_r =
new int[
mat.
N()];
481 etree =
new int[
mat.
M()];
511 dinfo<<
"LU factorization: dgssvx() returns info "<<
info<<std::endl;
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;
521 dinfo<<
"No of nonzeros in factor L = "<<
Lstore->nnz<<std::endl;
522 dinfo<<
"No of nonzeros in factor U = "<<
Ustore->nnz<<std::endl;
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!");
604#ifdef SUPERLU_MIN_VERSION_4_3
607 options.IterRefine=
DOUBLE;
611 &L, &U, work, lwork,
mB,
mX, &
rpg, &
rcond, &
ferr, &
berr,
632 dinfo<<
"Triangular solve: dgssvx() returns info "<<
info<<std::endl;
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;
654 template<
typename T,
typename A,
int n,
int m>
687#ifdef SUPERLU_MIN_VERSION_4_3
690 options.IterRefine=
DOUBLE;
694 &L, &U, work, lwork,
mB,
mX, &
rpg, &
rcond, &
ferr, &
berr,
698 dinfo<<
"Triangular solve: dgssvx() returns info "<<
info<<std::endl;
702 if ( options.IterRefine ) {
703 dinfo<<
"Iterative Refinement: steps="
704 <<
stat.RefineSteps<<
" FERR="<<
ferr<<
" BERR="<<
berr<<std::endl;
707 }
else if (
info > 0 && lwork == -1 ) {
708 dinfo<<
"** Estimated memory: "<<
info - n<<
" bytes"<<std::endl;
721 template<
typename T,
typename A,
int n,
int m>
727 template<
typename T,
typename A,
int n,
int m>
735#undef FIRSTCOL_OF_SNODE
Implementations of the inverse operator interface.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Templates characterizing the type of a solver.
Implementation of the BCRSMatrix class.
Matrix & mat
Definition matrixmatrix.hh:343
Definition basearray.hh:19
Statistics about compression achieved in implicit mode.
Definition bcrsmatrix.hh:81
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:412
size_type dim() const
dimension of the vector space
Definition bvector.hh:279
A vector of blocks with memory management.
Definition bvector.hh:309
Sequential overlapping Schwarz preconditioner.
Definition overlappingschwarz.hh:742
Definition overlappingschwarz.hh:683
derive error class from the base class in common
Definition istlexception.hh:16
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
Definition matrixutils.hh:25
Statistics about the application of an inverse operator.
Definition solver.hh:32
Abstract base class for all solvers.
Definition solver.hh:79
Definition solvertype.hh:14
@ value
Whether this is a direct solver.
Definition solvertype.hh:22
Definition solvertype.hh:28
@ value
whether the solver internally uses column compressed storage
Definition solvertype.hh:34
static void destroy(SuperMatrix *)
Definition superlu.hh:95
static void create(SuperMatrix *mat, int n, int m, float *dat, int n1, Stype_t stype, Dtype_t dtype, Mtype_t mtype)
Definition superlu.hh:88
static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree, char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix *B, SuperMatrix *X, float *rpg, float *rcond, float *ferr, float *berr, mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
Definition superlu.hh:102
static void querySpace(SuperMatrix *L, SuperMatrix *U, mem_usage_t *memusage)
Definition superlu.hh:124
const char * name()
Definition superlu.hh:374
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::BCRSMatrix< FieldMatrix< T, n, m >, A > Matrix
The matrix type.
Definition superlu.hh:296
SuperLUMatrix::size_type nnz() const
Definition superlu.hh:358
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
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > matrix_type
Definition superlu.hh:297
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
Definition supermatrix.hh:148