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

Made superlu work for float, too. 

[[Imported from SVN: r1403]]
parent 354e23f4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/superlu.hh
+ 79
− 21
View file @ 31f75b7b
......@@ -12,8 +12,10 @@
#endif
#ifdef SUPERLU_POST_2005_VERSION
#include "slu_ddefs.h"
//#include "slu_sdefs.h"
#else
#include "dsp_defs.h"
//#include "fsp_defs.h"
#endif
#include "solvers.hh"
#include "supermatrix.hh"
......@@ -104,7 +106,8 @@ namespace Dune
/**
* \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
*/
void apply (domain_type& x, range_type& b, double reduction, InverseOperatorResult& res)
void apply (domain_type& x, range_type& b, typename Dune::FieldTraits<T>::real_type reduction,
InverseOperatorResult& res)
{
apply(x,b,res);
}
......@@ -235,33 +238,33 @@ namespace Dune
perm_c = new int[mat.M()];
perm_r = new int[mat.N()];
etree = new int[mat.M()];
R = new double[mat.N()];
C = new double[mat.M()];
R = new T[mat.N()];
C = new T[mat.M()];
set_default_options(&options);
// Do the factorization
B.ncol=0;
B.Stype=SLU_DN;
B.Dtype=SLU_D;
B.Dtype= static_cast<Dtype_t>(GetSuperLUType<T>::type);
B.Mtype= SLU_GE;
DNformat fakeFormat;
fakeFormat.lda=mat.N();
B.Store=&fakeFormat;
X.Stype=SLU_DN;
X.Dtype=SLU_D;
X.Dtype=static_cast<Dtype_t>(GetSuperLUType<T>::type);
X.Mtype= SLU_GE;
X.ncol=0;
X.Store=&fakeFormat;
double rpg, rcond, ferr, berr;
T rpg, rcond, ferr, berr;
int info;
mem_usage_t memusage;
SuperLUStat_t stat;
StatInit(&stat);
dgssvx(&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);
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);
if(verbose) {
dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
......@@ -320,6 +323,58 @@ namespace Dune
options.Fact = FACTORED;
}
void createDenseSuperLUMatrix(SuperMatrix* B, int rows, int cols, double* b, int size,
Stype_t stype, Mtype_t mtype)
{
dCreate_Dense_Matrix(B, rows, cols, b, size, stype, SLU_D, mtype);
}
// Unfortunately SuperLU uses a lot of copy and paste in its headers.
// This results in some structs being declares in the headers of the float
// AND double version. To get around this we only include the double version
// and define the functions of the other versions as extern.
extern "C"
{
// single precision versions of SuperLU
void sCreate_Dense_Matrix(SuperMatrix* B, int rows, int cols, float* b, int size,
Stype_t stype, Dtype_t dtype, Mtype_t mtype);
void sgssvx(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 createDenseSuperLUMatrix(SuperMatrix* B, int rows, int cols, float* b, int size,
Stype_t stype, Mtype_t mtype)
{
sCreate_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)
{
dgssvx(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)
{
sgssvx(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)
......@@ -328,15 +383,18 @@ namespace Dune
DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
if(first) {
dCreate_Dense_Matrix(&B, mat.N(), 1, reinterpret_cast<T*>(&b[0]), mat.N(), SLU_DN, SLU_D, SLU_GE);
dCreate_Dense_Matrix(&X, mat.N(), 1, reinterpret_cast<T*>(&x[0]), mat.N(), SLU_DN, SLU_D, SLU_GE);
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]),
mat.N(), SLU_DN, SLU_GE);
first=false;
}else{
((DNformat*) B.Store)->nzval=&b[0];
((DNformat*)X.Store)->nzval=&x[0];
}
double rpg, rcond, ferr, berr;
T rpg, rcond, ferr, berr;
int info;
mem_usage_t memusage;
SuperLUStat_t stat;
......@@ -350,9 +408,9 @@ namespace Dune
*/
options.IterRefine=DOUBLE;
dgssvx(&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);
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);
res.iterations=1;
......@@ -398,15 +456,15 @@ namespace Dune
DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
if(first) {
dCreate_Dense_Matrix(&B, mat.N(), 1, b, mat.N(), SLU_DN, SLU_D, SLU_GE);
dCreate_Dense_Matrix(&X, mat.N(), 1, x, mat.N(), SLU_DN, SLU_D, SLU_GE);
createDenseSuperLUMatrix(&B, mat.N(), 1, b, mat.N(), SLU_DN, SLU_GE);
createDenseSuperLUMatrix(&X, mat.N(), 1, x, mat.N(), SLU_DN, SLU_GE);
first=false;
}else{
((DNformat*) B.Store)->nzval=b;
((DNformat*)X.Store)->nzval=x;
}
double rpg, rcond, ferr, berr;
T rpg, rcond, ferr, berr;
int info;
mem_usage_t memusage;
SuperLUStat_t stat;
......@@ -415,9 +473,9 @@ namespace Dune
options.IterRefine=DOUBLE;
dgssvx(&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);
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);
if(verbose) {
dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
......
dune/istl/supermatrix.hh
+ 32
− 2
View file @ 31f75b7b
......@@ -498,13 +498,43 @@ namespace Dune
}
}
// Unfortunately SuperLU uses a lot of copy and paste in its headers.
// This results in some structs being declares in the headers of the float
// AND double version. To get around this we only include the double version
// and define the functions of the other versions as extern.
extern "C"
{
// single precision versions of SuperMatrix creation
void
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 createCompColSuperMatrix(SuperMatrix *A, int m, int n, int nnz,
double *nzval, int *rowind, int *colptr,
Stype_t stype, Mtype_t mtype)
{
dCreate_CompCol_Matrix(A, m, n, nnz, nzval, rowind, colptr, stype,
SLU_D, mtype);
}
void createCompColSuperMatrix(SuperMatrix *A, int m, int n, int nnz,
float *nzval, int *rowind, int *colptr,
Stype_t stype, Mtype_t mtype)
{
sCreate_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
{
delete[] marker;
marker=0;
dCreate_CompCol_Matrix(&mat->A, mat->N_, mat->M_, mat->colstart[cols],
mat->values, mat->rowindex, mat->colstart, SLU_NC, static_cast<Dtype_t>(GetSuperLUType<T>::type), SLU_GE);
createCompColSuperMatrix(&mat->A, mat->N_, mat->M_, mat->colstart[cols],
mat->values, mat->rowindex, mat->colstart, SLU_NC, SLU_GE);
}
template<class F, class MRS>
......
dune/istl/test/superlutest.cc
+ 22
− 12
View file @ 31f75b7b
......@@ -8,20 +8,12 @@
#include <laplacian.hh>
#include <dune/common/timer.hh>
#include <dune/istl/superlu.hh>
int main(int argc, char** argv)
{
const int BS=1;
std::size_t N=100;
if(argc>1)
N = atoi(argv[1]);
std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
typedef Dune::FieldMatrix<double,BS,BS> MatrixBlock;
template<class T, int BS>
void testSuperLU(std::size_t N){
typedef Dune::FieldMatrix<T,BS,BS> MatrixBlock;
typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
typedef Dune::FieldVector<double,BS> VectorBlock;
typedef Dune::FieldVector<T,BS> VectorBlock;
typedef Dune::BlockVector<VectorBlock> Vector;
typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
......@@ -55,3 +47,21 @@ int main(int argc, char** argv)
solver1.apply(x,b, res);
}
int main(int argc, char** argv)
{
const int BS=1;
std::size_t N=100;
if(argc>1)
N = atoi(argv[1]);
std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
testSuperLU<double,BS>(N);
testSuperLU<float,BS>(N);
//testSuperLU<std::complex<double>,1>(N);
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment