fdurastante
/
amg4psblas
mirror of https://github.com/sfilippone/amg4psblas.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

mld2p4-2:

Browse Source 
mlprec/impl/mld_cslu_interface.c
 mlprec/impl/mld_dslu_interface.c
 mlprec/impl/mld_sslu_interface.c
 mlprec/impl/mld_zslu_interface.c
 mlprec/mld_c_slu_solver.F90
 mlprec/mld_d_slu_solver.F90
 mlprec/mld_s_slu_solver.F90
 mlprec/mld_z_slu_solver.F90


Fixed interface to SuperLU solver; SuperLU now working.
stopcriterion
Salvatore Filippone 12 years ago
parent 6b8ff2b597
commit 91d3f5a043
8 changed files with 105 additions and 125 deletions
Show Stats Download Patch File Download Diff File

39
mlprec/impl/mld_cslu_interface.c
Unescape Escape View File

@ -38,7 +38,7 @@
*
* File: mld_cslu_interface.c
*
* Functions: mld_cslu_fact_, mld_cslu_solve_, mld_cslu_free_.
* Functions: mld_cslu_fact, mld_cslu_solve, mld_cslu_free.
*
* This file is an interface to the SuperLU routines for sparse factorization and
* solve. It was obtained by modifying the c_fortran_cgssv.c file from the SuperLU
@ -110,15 +110,13 @@ typedef struct {
int
mld_cslu_fact(int n, int nnz,
int mld_cslu_fact(int n, int nnz,
#ifdef HAVE_SLU_
complex *values,
complex *values,
#else
void *values,
void *values,
#endif
int *rowptr, int *colind, void **f_factors)
int *rowptr, int *colind, void **f_factors)
{
/*
* This routine can be called from Fortran.
@ -140,7 +138,6 @@ mld_cslu_fact(int n, int nnz,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
float drop_tol = 0.0;
mem_usage_t mem_usage;
superlu_options_t options;
SuperLUStat_t stat;
@ -180,7 +177,7 @@ mld_cslu_fact(int n, int nnz,
panel_size = sp_ienv(1);
relax = sp_ienv(2);
cgstrf(&options, &AC, drop_tol, relax, panel_size,
cgstrf(&options, &AC, relax, panel_size,
etree, NULL, 0, perm_c, perm_r, L, U, &stat, &info);
if ( info == 0 ) {
@ -191,17 +188,15 @@ mld_cslu_fact(int n, int nnz,
printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
#endif
} else {
printf("cgstrf() error returns INFO= %d\n", info);
if ( info <= n ) { /* factorization completes */
cQuerySpace(L, U, &mem_usage);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
}
}
@ -226,15 +221,13 @@ mld_cslu_fact(int n, int nnz,
}
int
mld_cslu_solve(int itrans, int n, int nrhs,
int mld_cslu_solve(int itrans, int n, int nrhs,
#ifdef HAVE_SLU_
complex *b,
complex *b,
#else
void *b,
void *b,
#endif
int ldb,void *f_factors)
int ldb,void *f_factors)
{
/*
* This routine can be called from Fortran.
@ -296,9 +289,7 @@ mld_cslu_solve(int itrans, int n, int nrhs,
}
int
mld_cslu_free(void *f_factors)
int mld_cslu_free(void *f_factors)
{
/*
* This routine can be called from Fortran.

39
mlprec/impl/mld_dslu_interface.c
Unescape Escape View File

@ -38,7 +38,7 @@
*
* File: mld_slu_interface.c
*
* Functions: mld_dslu_fact_, mld_dslu_solve_, mld_dslu_free_.
* Functions: mld_dslu_fact, mld_dslu_solve, mld_dslu_free.
*
* This file is an interface to the SuperLU routines for sparse factorization and
* solve. It was obtained by modifying the c_fortran_dgssv.c file from the SuperLU
@ -110,10 +110,8 @@ typedef struct {
int
mld_dslu_fact(int n, int nnz, double *values,
int *rowptr, int *colind, void **f_factors)
int mld_dslu_fact(int n, int nnz, double *values,
int *rowptr, int *colind, void **f_factors)
{
/*
* This routine can be called from Fortran.
@ -135,7 +133,6 @@ mld_dslu_fact(int n, int nnz, double *values,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
double drop_tol = 0.0;
mem_usage_t mem_usage;
superlu_options_t options;
SuperLUStat_t stat;
@ -166,8 +163,9 @@ mld_dslu_fact(int n, int nnz, double *values,
* permc_spec = 2: minimum degree on structure of A'+A
* permc_spec = 3: approximate minimum degree for unsymmetric matrices
*/
options.ColPerm=2;
permc_spec = options.ColPerm;
/* options.ColPerm=2; */
/* permc_spec = options.ColPerm; */
permc_spec=0;
get_perm_c(permc_spec, &A, perm_c);
sp_preorder(&options, &A, perm_c, etree, &AC);
@ -175,7 +173,7 @@ mld_dslu_fact(int n, int nnz, double *values,
panel_size = sp_ienv(1);
relax = sp_ienv(2);
dgstrf(&options, &AC, drop_tol, relax, panel_size,
dgstrf(&options, &AC, relax, panel_size,
etree, NULL, 0, perm_c, perm_r, L, U, &stat, &info);
if ( info == 0 ) {
@ -186,17 +184,15 @@ mld_dslu_fact(int n, int nnz, double *values,
printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
#endif
} else {
printf("dgstrf() error returns INFO= %d\n", info);
if ( info <= n ) { /* factorization completes */
dQuerySpace(L, U, &mem_usage);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
}
}
@ -221,10 +217,8 @@ mld_dslu_fact(int n, int nnz, double *values,
}
int
mld_dslu_solve(int itrans, int n, int nrhs, double *b, int ldb,
void *f_factors)
int mld_dslu_solve(int itrans, int n, int nrhs, double *b, int ldb,
void *f_factors)
{
/*
* This routine can be called from Fortran.
@ -242,7 +236,6 @@ mld_dslu_solve(int itrans, int n, int nrhs, double *b, int ldb,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
double drop_tol = 0.0;
SuperLUStat_t stat;
factors_t *LUfactors;
@ -267,7 +260,7 @@ mld_dslu_solve(int itrans, int n, int nrhs, double *b, int ldb,
dCreate_Dense_Matrix(&B, n, nrhs, b, ldb, SLU_DN, SLU_D, SLU_GE);
/* Solve the system A*X=B, overwriting B with X. */
dgstrs (trans, L, U, perm_c, perm_r, &B, &stat, &info);
dgstrs(trans, L, U, perm_c, perm_r, &B, &stat, &info);
Destroy_SuperMatrix_Store(&B);
StatFree(&stat);
@ -279,9 +272,7 @@ mld_dslu_solve(int itrans, int n, int nrhs, double *b, int ldb,
}
int
mld_dslu_free(void *f_factors)
int mld_dslu_free(void *f_factors)
{
/*
* This routine can be called from Fortran.

32
mlprec/impl/mld_sslu_interface.c
Unescape Escape View File

@ -38,7 +38,7 @@
*
* File: mld_slu_interface.c
*
* Functions: mld_sslu_fact_, mld_sslu_solve_, mld_sslu_free_.
* Functions: mld_sslu_fact, mld_sslu_solve, mld_sslu_free.
*
* This file is an interface to the SuperLU routines for sparse factorization and
* solve. It was obtained by modifying the c_fortran_dgssv.c file from the SuperLU
@ -110,10 +110,8 @@ typedef struct {
int
mld_sslu_fact(int n, int nnz, float *values,
int *rowptr, int *colind, void **f_factors)
int mld_sslu_fact(int n, int nnz, float *values,
int *rowptr, int *colind, void **f_factors)
{
/*
* This routine can be called from Fortran.
@ -135,7 +133,6 @@ mld_sslu_fact(int n, int nnz, float *values,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
float drop_tol = 0.0;
mem_usage_t mem_usage;
superlu_options_t options;
SuperLUStat_t stat;
@ -175,7 +172,7 @@ mld_sslu_fact(int n, int nnz, float *values,
panel_size = sp_ienv(1);
relax = sp_ienv(2);
sgstrf(&options, &AC, drop_tol, relax, panel_size,
sgstrf(&options, &AC, relax, panel_size,
etree, NULL, 0, perm_c, perm_r, L, U, &stat, &info);
if ( info == 0 ) {
@ -186,17 +183,15 @@ mld_sslu_fact(int n, int nnz, float *values,
printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
#endif
} else {
printf("sgstrf() error returns INFO= %d\n", info);
if ( info <= n ) { /* factorization completes */
sQuerySpace(L, U, &mem_usage);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
}
}
@ -221,10 +216,8 @@ mld_sslu_fact(int n, int nnz, float *values,
}
int
mld_sslu_solve(int itrans, int n, int nrhs, float *b, int ldb,
void *f_factors)
int mld_sslu_solve(int itrans, int n, int nrhs, float *b, int ldb,
void *f_factors)
{
/*
* This routine can be called from Fortran.
@ -242,7 +235,6 @@ mld_sslu_solve(int itrans, int n, int nrhs, float *b, int ldb,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
float drop_tol = 0.0;
SuperLUStat_t stat;
factors_t *LUfactors;
@ -279,9 +271,7 @@ mld_sslu_solve(int itrans, int n, int nrhs, float *b, int ldb,
}
int
mld_sslu_free(void *f_factors)
int mld_sslu_free(void *f_factors)
{
/*
* This routine can be called from Fortran.

32
mlprec/impl/mld_zslu_interface.c
Unescape Escape View File

@ -38,7 +38,7 @@
*
* File: mld_zslu_interface.c
*
* Functions: mld_zslu_fact_, mld_zslu_solve_, mld_zslu_free_.
* Functions: mld_zslu_fact, mld_zslu_solve, mld_zslu_free.
*
* This file is an interface to the SuperLU routines for sparse factorization and
* solve. It was obtained by modifying the c_fortran_zgssv.c file from the SuperLU
@ -109,8 +109,7 @@ typedef struct {
int
mld_zslu_fact(int n, int nnz,
int mld_zslu_fact(int n, int nnz,
#ifdef HAVE_SLU_
doublecomplex *values,
#else
@ -139,7 +138,6 @@ mld_zslu_fact(int n, int nnz,
NCformat *Ustore;
int i, panel_size, permc_spec, relax;
trans_t trans;
double drop_tol = 0.0;
mem_usage_t mem_usage;
superlu_options_t options;
SuperLUStat_t stat;
@ -179,7 +177,7 @@ mld_zslu_fact(int n, int nnz,
panel_size = sp_ienv(1);
relax = sp_ienv(2);
zgstrf(&options, &AC, drop_tol, relax, panel_size,
zgstrf(&options, &AC, relax, panel_size,
etree, NULL, 0, perm_c, perm_r, L, U, &stat, &info);
if ( info == 0 ) {
@ -190,17 +188,15 @@ mld_zslu_fact(int n, int nnz,
printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
#endif
} else {
printf("zgstrf() error returns INFO= %d\n", info);
if ( info <= n ) { /* factorization completes */
zQuerySpace(L, U, &mem_usage);
printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
mem_usage.expansions);
printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
}
}
@ -225,15 +221,13 @@ mld_zslu_fact(int n, int nnz,
}
int
mld_zslu_solve(int itrans, int n, int nrhs,
int mld_zslu_solve(int itrans, int n, int nrhs,
#ifdef HAVE_SLU_
doublecomplex *b,
doublecomplex *b,
#else
void *b,
void *b,
#endif
int ldb,void *f_factors)
int ldb,void *f_factors)
{
/*
* This routine can be called from Fortran.
@ -295,9 +289,7 @@ mld_zslu_solve(int itrans, int n, int nrhs,
}
int
mld_zslu_free(void *f_factors)
int mld_zslu_free(void *f_factors)
{
/*
* This routine can be called from Fortran.

22
mlprec/mld_c_slu_solver.F90
Unescape Escape View File

@ -95,12 +95,12 @@ module mld_c_slu_solver
end interface
interface
function mld_cslu_solve(itrans,n,x, b, ldb, lufactors)&
function mld_cslu_solve(itrans,n,nrhs,b,ldb,lufactors)&
& bind(c,name='mld_cslu_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
complex(c_float_complex) :: x(*), b(ldb,*)
integer(c_int), value :: itrans,n,nrhs,ldb
complex(c_float_complex) :: b(ldb,*)
type(c_ptr), value :: lufactors
end function mld_cslu_solve
end interface
@ -162,23 +162,27 @@ contains
end if
endif
ww(1:n_row) = x(1:n_row)
select case(trans_)
case('N')
info = mld_cslu_solve(0,n_row,ww,x,n_row,sv%lufactors)
info = mld_cslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
case('T')
info = mld_cslu_solve(1,n_row,ww,x,n_row,sv%lufactors)
info = mld_cslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
case('C')
info = mld_cslu_solve(2,n_row,ww,x,n_row,sv%lufactors)
info = mld_cslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
call psb_errpush(psb_err_internal_error_, &
& name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == psb_success_) &
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in subsolve')
goto 9999
endif

22
mlprec/mld_d_slu_solver.F90
Unescape Escape View File

@ -95,12 +95,12 @@ module mld_d_slu_solver
end interface
interface
function mld_dslu_solve(itrans,n,x, b, ldb, lufactors)&
function mld_dslu_solve(itrans,n,nrhs,b,ldb,lufactors)&
& bind(c,name='mld_dslu_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
real(c_double) :: x(*), b(ldb,*)
integer(c_int), value :: itrans,n,nrhs,ldb
real(c_double) :: b(ldb,*)
type(c_ptr), value :: lufactors
end function mld_dslu_solve
end interface
@ -162,23 +162,27 @@ contains
end if
endif
ww(1:n_row) = x(1:n_row)
select case(trans_)
case('N')
info = mld_dslu_solve(0,n_row,ww,x,n_row,sv%lufactors)
info = mld_dslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
case('T')
info = mld_dslu_solve(1,n_row,ww,x,n_row,sv%lufactors)
info = mld_dslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
case('C')
info = mld_dslu_solve(2,n_row,ww,x,n_row,sv%lufactors)
info = mld_dslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
call psb_errpush(psb_err_internal_error_, &
& name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == psb_success_) &
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in subsolve')
goto 9999
endif

22
mlprec/mld_s_slu_solver.F90
Unescape Escape View File

@ -95,12 +95,12 @@ module mld_s_slu_solver
end interface
interface
function mld_sslu_solve(itrans,n,x, b, ldb, lufactors)&
function mld_sslu_solve(itrans,n,nrhs,b,ldb,lufactors)&
& bind(c,name='mld_sslu_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
real(c_float) :: x(*), b(ldb,*)
integer(c_int), value :: itrans,n,nrhs,ldb
real(c_float) :: b(ldb,*)
type(c_ptr), value :: lufactors
end function mld_sslu_solve
end interface
@ -162,23 +162,27 @@ contains
end if
endif
ww(1:n_row) = x(1:n_row)
select case(trans_)
case('N')
info = mld_sslu_solve(0,n_row,ww,x,n_row,sv%lufactors)
info = mld_sslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
case('T')
info = mld_sslu_solve(1,n_row,ww,x,n_row,sv%lufactors)
info = mld_sslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
case('C')
info = mld_sslu_solve(2,n_row,ww,x,n_row,sv%lufactors)
info = mld_sslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
call psb_errpush(psb_err_internal_error_, &
& name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == psb_success_) &
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in subsolve')
goto 9999
endif

22
mlprec/mld_z_slu_solver.F90
Unescape Escape View File

@ -95,12 +95,12 @@ module mld_z_slu_solver
end interface
interface
function mld_zslu_solve(itrans,n,x, b, ldb, lufactors)&
function mld_zslu_solve(itrans,n,nrhs,b,ldb,lufactors)&
& bind(c,name='mld_zslu_solve') result(info)
use iso_c_binding
integer(c_int) :: info
integer(c_int), value :: itrans,n,ldb
complex(c_double_complex) :: x(*), b(ldb,*)
integer(c_int), value :: itrans,n,nrhs,ldb
complex(c_double_complex) :: b(ldb,*)
type(c_ptr), value :: lufactors
end function mld_zslu_solve
end interface
@ -162,23 +162,27 @@ contains
end if
endif
ww(1:n_row) = x(1:n_row)
select case(trans_)
case('N')
info = mld_zslu_solve(0,n_row,ww,x,n_row,sv%lufactors)
info = mld_zslu_solve(0,n_row,1,ww,n_row,sv%lufactors)
case('T')
info = mld_zslu_solve(1,n_row,ww,x,n_row,sv%lufactors)
info = mld_zslu_solve(1,n_row,1,ww,n_row,sv%lufactors)
case('C')
info = mld_zslu_solve(2,n_row,ww,x,n_row,sv%lufactors)
info = mld_zslu_solve(2,n_row,1,ww,n_row,sv%lufactors)
case default
call psb_errpush(psb_err_internal_error_,name,a_err='Invalid TRANS in ILU subsolve')
call psb_errpush(psb_err_internal_error_, &
& name,a_err='Invalid TRANS in ILU subsolve')
goto 9999
end select
if (info == psb_success_) call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info == psb_success_) &
& call psb_geaxpby(alpha,ww,beta,y,desc_data,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Error in subsolve')
call psb_errpush(psb_err_internal_error_,&
& name,a_err='Error in subsolve')
goto 9999
endif

Write Preview
Loading…
Cancel
Save