Skip to content
Snippets Groups Projects
Commit b273a8ee authored by Markus Blatt's avatar Markus Blatt
Browse files

Allow non-double values, too. 

[[Imported from SVN: r1407]]
parent 2d5ac6fd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/superlu.hh
+ 72
− 9
View file @ b273a8ee
......@@ -147,8 +147,10 @@ namespace Dune
SuperLUMatrix mat;
SuperMatrix L, U, B, X;
int *perm_c, *perm_r, *etree;
T *R, *C;
T *bstore;
// get the corresponding superlu type (different from the
// original one for complex<(double|float)>
typedef typename GetSuperLUStorageType<T>::type storage_type;
storage_type *R, *C;
superlu_options_t options;
char equed;
void *work;
......@@ -238,8 +240,11 @@ namespace Dune
perm_c = new int[mat.M()];
perm_r = new int[mat.N()];
etree = new int[mat.M()];
R = new T[mat.N()];
C = new T[mat.M()];
// get the corresponding superlu type (different from the
// original one for complex<(double|float)>
typedef typename GetSuperLUStorageType<T>::type storage_type;
R = new storage_type[mat.N()];
C = new storage_type[mat.M()];
set_default_options(&options);
// Do the factorization
......@@ -264,7 +269,7 @@ namespace Dune
StatInit(&stat);
applySuperLU(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
&L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
&berr, &memusage, &stat, &info);
&berr, &memusage, &stat, &info, T());
if(verbose) {
dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
......@@ -345,6 +350,26 @@ namespace Dune
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);
// single precision versions of SuperLU
void cCreate_Dense_Matrix(SuperMatrix* B, int rows, int cols, complex* b, int size,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
void cgssvx(superlu_options_t *options, SuperMatrix *mat, int *permc, int *permr, 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);
void zCreate_Dense_Matrix(SuperMatrix* B, int rows, int cols, doublecomplex* b, int size,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
void zgssvx(superlu_options_t *options, SuperMatrix *mat, int *permc, int *permr, 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);
}
void createDenseSuperLUMatrix(SuperMatrix* B, int rows, int cols, float* b, int size,
......@@ -353,11 +378,23 @@ namespace Dune
sCreate_Dense_Matrix(B, rows, cols, b, size, stype, SLU_S, mtype);
}
void createDenseSuperLUMatrix(SuperMatrix* B, int rows, int cols, complex* b, int size,
Stype_t stype, Mtype_t mtype)
{
cCreate_Dense_Matrix(B, rows, cols, b, size, stype, SLU_S, mtype);
}
void createDenseSuperLUMatrix(SuperMatrix* B, int rows, int cols, doublecomplex* b, int size,
Stype_t stype, Mtype_t mtype)
{
zCreate_Dense_Matrix(B, rows, cols, b, size, stype, SLU_S, mtype);
}
void applySuperLU(superlu_options_t *options, SuperMatrix *mat, int *permc, int *permr, 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)
mem_usage_t *memusage, SuperLUStat_t *stat, int *info, const double&)
{
dgssvx(options, mat, permc, permr, etree, equed, R, C,
L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
......@@ -369,12 +406,36 @@ namespace Dune
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)
mem_usage_t *memusage, SuperLUStat_t *stat, int *info, const float&)
{
sgssvx(options, mat, permc, permr, etree, equed, R, C,
L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
memusage, stat, info);
}
void applySuperLU(superlu_options_t *options, SuperMatrix *mat, int *permc, int *permr, 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, const doublecomplex&)
{
zgssvx(options, mat, permc, permr, etree, equed, R, C,
L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
memusage, stat, info);
}
void applySuperLU(superlu_options_t *options, SuperMatrix *mat, int *permc, int *permr, 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, const complex&)
{
cgssvx(options, mat, permc, permr, etree, equed, R, C,
L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
memusage, stat, info);
}
template<typename T, typename A, int n, int m>
void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
......@@ -384,6 +445,7 @@ namespace Dune
if(first) {
assert(mat.N()<=static_cast<std::size_t>(std::numeric_limits<int>::max()));
createDenseSuperLUMatrix(&B, mat.N(), 1, reinterpret_cast<T*>(&b[0]),
mat.N(), SLU_DN, SLU_GE);
createDenseSuperLUMatrix(&X, mat.N(), 1, reinterpret_cast<T*>(&x[0]),
......@@ -408,9 +470,10 @@ namespace Dune
*/
options.IterRefine=DOUBLE;
applySuperLU(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
&L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr, &berr,
&memusage, &stat, &info);
&memusage, &stat, &info, T());
res.iterations=1;
......@@ -475,7 +538,7 @@ namespace Dune
applySuperLU(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
&L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr, &berr,
&memusage, &stat, &info);
&memusage, &stat, &info, T());
if(verbose) {
dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
......
dune/istl/supermatrix.hh
+ 72
− 4
View file @ b273a8ee
......@@ -6,6 +6,8 @@
#if HAVE_SUPERLU
#ifdef SUPERLU_POST_2005_VERSION
#include "slu_ddefs.h"
#include "slu_dcomplex.h"
#include "slu_scomplex.h"
#else
#include "dsp_defs.h"
#endif
......@@ -18,6 +20,41 @@
namespace Dune
{
template<typename T>
void copySuperLUValue(T& t, const T& s)
{
t=s;
}
void copySuperLUValue(doublecomplex& t, const std::complex<double>& s)
{
t.r = s.real();
t.i = s.imag();
}
void copySuperLUValue(complex& t, const std::complex<float>& s)
{
t.r = s.real();
t.i = s.imag();
}
template<typename T>
struct GetSuperLUStorageType
{
typedef T type;
};
template<>
struct GetSuperLUStorageType<std::complex<float> >
{
typedef complex type;
};
template<>
struct GetSuperLUStorageType<std::complex<double> >
{
typedef doublecomplex type;
};
/**
* @brief Provides access to an iterator over all matrix rows.
*
......@@ -281,7 +318,10 @@ namespace Dune
void setMatrix(const Matrix& mat);
int N_, M_, Nnz_;
B* values;
// get the corresponding superlu type (different from the
// original one for complex<(double|float)>
typedef typename GetSuperLUStorageType<B>::type storage_type;
storage_type* values;
int* rowindex;
int* colstart;
SuperMatrix A;
......@@ -433,8 +473,11 @@ namespace Dune
void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage() const
{
mat->Nnz_*=n*m;
// get the corresponding superlu type (different from the
// original one for complex<(double|float)>
typedef typename GetSuperLUStorageType<T>::type storage_type;
// initialize data
mat->values=new T[mat->Nnz_];
mat->values=new storage_type[mat->Nnz_];
mat->rowindex=new int[mat->Nnz_];
mat->colstart=new int[cols+1];
}
......@@ -492,7 +535,7 @@ namespace Dune
assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
assert((int)marker[colindex*m+j]<mat->Nnz_);
mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
mat->values[marker[colindex*m+j]]=(*col)[i][j];
copySuperLUValue(mat->values[marker[colindex*m+j]],(*col)[i][j]);
++marker[colindex*m+j]; // index for next entry in column
}
}
......@@ -509,6 +552,15 @@ namespace Dune
sCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
float *nzval, int *rowind, int *colptr,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
void
cCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
complex *nzval, int *rowind, int *colptr,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
void
zCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
doublecomplex *nzval, int *rowind, int *colptr,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
}
......@@ -528,6 +580,22 @@ namespace Dune
SLU_S, mtype);
}
void createCompColSuperMatrix(SuperMatrix *A, int m, int n, int nnz,
complex *nzval, int *rowind, int *colptr,
Stype_t stype, Mtype_t mtype)
{
cCreate_CompCol_Matrix(A, m, n, nnz, nzval, rowind, colptr, stype,
SLU_D, mtype);
}
void createCompColSuperMatrix(SuperMatrix *A, int m, int n, int nnz,
doublecomplex *nzval, int *rowind, int *colptr,
Stype_t stype, Mtype_t mtype)
{
zCreate_CompCol_Matrix(A, m, n, nnz, nzval, rowind, colptr, stype,
SLU_S, mtype);
}
template<class T, class A, int n, int m>
void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::createMatrix() const
{
......@@ -678,7 +746,7 @@ namespace Dune
}
if(Nnz_>0) {
values = new B[Nnz_];
values = new storage_type[Nnz_];
rowindex = new int[Nnz_];
for(int i=0; i<Nnz_; ++i)
......
dune/istl/test/superlutest.cc
+ 1
− 2
View file @ b273a8ee
......@@ -62,6 +62,5 @@ int main(int argc, char** argv)
testSuperLU<float,BS>(N);
//testSuperLU<std::complex<double>,1>(N);
testSuperLU<std::complex<double>,1>(N);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment