From 8fed44deeae4f491b24f7d973928330123079cb3 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 12 Jul 2013 14:51:05 +0000 Subject: [PATCH] mld2p4-2: README mlprec/impl/mld_cmlprec_aply.f90 mlprec/impl/mld_cprecaply.f90 mlprec/impl/mld_cslud_interface.c mlprec/impl/mld_dmlprec_aply.f90 mlprec/impl/mld_dprecaply.f90 mlprec/impl/mld_dslud_interface.c mlprec/impl/mld_smlprec_aply.f90 mlprec/impl/mld_sprecaply.f90 mlprec/impl/mld_sslud_interface.c mlprec/impl/mld_zmlprec_aply.f90 mlprec/impl/mld_zprecaply.f90 mlprec/impl/mld_zslud_interface.c mlprec/impl/smoother/mld_c_as_smoother_apply.f90 mlprec/impl/smoother/mld_c_base_smoother_apply.f90 mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 mlprec/impl/smoother/mld_d_as_smoother_apply.f90 mlprec/impl/smoother/mld_d_base_smoother_apply.f90 mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 mlprec/impl/smoother/mld_s_as_smoother_apply.f90 mlprec/impl/smoother/mld_s_base_smoother_apply.f90 mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 mlprec/impl/smoother/mld_z_as_smoother_apply.f90 mlprec/impl/smoother/mld_z_base_smoother_apply.f90 mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 mlprec/impl/solver/mld_c_base_solver_apply.f90 mlprec/impl/solver/mld_c_diag_solver_apply.f90 mlprec/impl/solver/mld_c_id_solver_apply.f90 mlprec/impl/solver/mld_c_ilu_solver_apply.f90 mlprec/impl/solver/mld_d_base_solver_apply.f90 mlprec/impl/solver/mld_d_diag_solver_apply.f90 mlprec/impl/solver/mld_d_id_solver_apply.f90 mlprec/impl/solver/mld_d_ilu_solver_apply.f90 mlprec/impl/solver/mld_s_base_solver_apply.f90 mlprec/impl/solver/mld_s_diag_solver_apply.f90 mlprec/impl/solver/mld_s_id_solver_apply.f90 mlprec/impl/solver/mld_s_ilu_solver_apply.f90 mlprec/impl/solver/mld_z_base_solver_apply.f90 mlprec/impl/solver/mld_z_diag_solver_apply.f90 mlprec/impl/solver/mld_z_id_solver_apply.f90 mlprec/impl/solver/mld_z_ilu_solver_apply.f90 mlprec/mld_c_as_smoother.f90 mlprec/mld_c_base_smoother_mod.f90 mlprec/mld_c_base_solver_mod.f90 mlprec/mld_c_diag_solver.f90 mlprec/mld_c_id_solver.f90 mlprec/mld_c_ilu_solver.f90 mlprec/mld_c_jac_smoother.f90 mlprec/mld_c_prec_type.f90 mlprec/mld_c_slu_solver.F90 mlprec/mld_c_sludist_solver.F90 mlprec/mld_c_umf_solver.F90 mlprec/mld_d_as_smoother.f90 mlprec/mld_d_base_smoother_mod.f90 mlprec/mld_d_base_solver_mod.f90 mlprec/mld_d_diag_solver.f90 mlprec/mld_d_id_solver.f90 mlprec/mld_d_ilu_solver.f90 mlprec/mld_d_jac_smoother.f90 mlprec/mld_d_prec_type.f90 mlprec/mld_d_slu_solver.F90 mlprec/mld_d_sludist_solver.F90 mlprec/mld_d_umf_solver.F90 mlprec/mld_s_as_smoother.f90 mlprec/mld_s_base_smoother_mod.f90 mlprec/mld_s_base_solver_mod.f90 mlprec/mld_s_diag_solver.f90 mlprec/mld_s_id_solver.f90 mlprec/mld_s_ilu_solver.f90 mlprec/mld_s_jac_smoother.f90 mlprec/mld_s_prec_type.f90 mlprec/mld_s_slu_solver.F90 mlprec/mld_s_sludist_solver.F90 mlprec/mld_s_umf_solver.F90 mlprec/mld_z_as_smoother.f90 mlprec/mld_z_base_smoother_mod.f90 mlprec/mld_z_base_solver_mod.f90 mlprec/mld_z_diag_solver.f90 mlprec/mld_z_id_solver.f90 mlprec/mld_z_ilu_solver.f90 mlprec/mld_z_jac_smoother.f90 mlprec/mld_z_prec_type.f90 mlprec/mld_z_slu_solver.F90 mlprec/mld_z_sludist_solver.F90 mlprec/mld_z_umf_solver.F90 Fix SuperLU_Dist. SuperLU does not work completely yet. Unify INTENT(INOUT) on solver_apply. --- README | 6 +- mlprec/impl/mld_cmlprec_aply.f90 | 4 +- mlprec/impl/mld_cprecaply.f90 | 4 +- mlprec/impl/mld_cslud_interface.c | 180 ++++++------------ mlprec/impl/mld_dmlprec_aply.f90 | 4 +- mlprec/impl/mld_dprecaply.f90 | 4 +- mlprec/impl/mld_dslud_interface.c | 13 +- mlprec/impl/mld_smlprec_aply.f90 | 4 +- mlprec/impl/mld_sprecaply.f90 | 4 +- mlprec/impl/mld_sslud_interface.c | 178 ++++++----------- mlprec/impl/mld_zmlprec_aply.f90 | 4 +- mlprec/impl/mld_zprecaply.f90 | 4 +- mlprec/impl/mld_zslud_interface.c | 155 +++++---------- .../impl/smoother/mld_c_as_smoother_apply.f90 | 2 +- .../smoother/mld_c_base_smoother_apply.f90 | 2 +- .../smoother/mld_c_jac_smoother_apply.f90 | 2 +- .../impl/smoother/mld_d_as_smoother_apply.f90 | 2 +- .../smoother/mld_d_base_smoother_apply.f90 | 2 +- .../smoother/mld_d_jac_smoother_apply.f90 | 2 +- .../impl/smoother/mld_s_as_smoother_apply.f90 | 2 +- .../smoother/mld_s_base_smoother_apply.f90 | 2 +- .../smoother/mld_s_jac_smoother_apply.f90 | 2 +- .../impl/smoother/mld_z_as_smoother_apply.f90 | 2 +- .../smoother/mld_z_base_smoother_apply.f90 | 2 +- .../smoother/mld_z_jac_smoother_apply.f90 | 2 +- .../impl/solver/mld_c_base_solver_apply.f90 | 2 +- .../impl/solver/mld_c_diag_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_c_id_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_c_ilu_solver_apply.f90 | 2 +- .../impl/solver/mld_d_base_solver_apply.f90 | 2 +- .../impl/solver/mld_d_diag_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_d_id_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_d_ilu_solver_apply.f90 | 2 +- .../impl/solver/mld_s_base_solver_apply.f90 | 2 +- .../impl/solver/mld_s_diag_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_s_id_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_s_ilu_solver_apply.f90 | 2 +- .../impl/solver/mld_z_base_solver_apply.f90 | 2 +- .../impl/solver/mld_z_diag_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_z_id_solver_apply.f90 | 2 +- mlprec/impl/solver/mld_z_ilu_solver_apply.f90 | 2 +- mlprec/mld_c_as_smoother.f90 | 2 +- mlprec/mld_c_base_smoother_mod.f90 | 2 +- mlprec/mld_c_base_solver_mod.f90 | 2 +- mlprec/mld_c_diag_solver.f90 | 2 +- mlprec/mld_c_id_solver.f90 | 2 +- mlprec/mld_c_ilu_solver.f90 | 2 +- mlprec/mld_c_jac_smoother.f90 | 2 +- mlprec/mld_c_prec_type.f90 | 4 +- mlprec/mld_c_slu_solver.F90 | 2 +- mlprec/mld_c_sludist_solver.F90 | 57 +++--- mlprec/mld_c_umf_solver.F90 | 2 +- mlprec/mld_d_as_smoother.f90 | 2 +- mlprec/mld_d_base_smoother_mod.f90 | 2 +- mlprec/mld_d_base_solver_mod.f90 | 2 +- mlprec/mld_d_diag_solver.f90 | 2 +- mlprec/mld_d_id_solver.f90 | 2 +- mlprec/mld_d_ilu_solver.f90 | 2 +- mlprec/mld_d_jac_smoother.f90 | 2 +- mlprec/mld_d_prec_type.f90 | 4 +- mlprec/mld_d_slu_solver.F90 | 8 +- mlprec/mld_d_sludist_solver.F90 | 19 +- mlprec/mld_d_umf_solver.F90 | 2 +- mlprec/mld_s_as_smoother.f90 | 2 +- mlprec/mld_s_base_smoother_mod.f90 | 2 +- mlprec/mld_s_base_solver_mod.f90 | 2 +- mlprec/mld_s_diag_solver.f90 | 2 +- mlprec/mld_s_id_solver.f90 | 2 +- mlprec/mld_s_ilu_solver.f90 | 2 +- mlprec/mld_s_jac_smoother.f90 | 2 +- mlprec/mld_s_prec_type.f90 | 4 +- mlprec/mld_s_slu_solver.F90 | 2 +- mlprec/mld_s_sludist_solver.F90 | 57 +++--- mlprec/mld_s_umf_solver.F90 | 2 +- mlprec/mld_z_as_smoother.f90 | 2 +- mlprec/mld_z_base_smoother_mod.f90 | 2 +- mlprec/mld_z_base_solver_mod.f90 | 2 +- mlprec/mld_z_diag_solver.f90 | 2 +- mlprec/mld_z_id_solver.f90 | 2 +- mlprec/mld_z_ilu_solver.f90 | 2 +- mlprec/mld_z_jac_smoother.f90 | 2 +- mlprec/mld_z_prec_type.f90 | 4 +- mlprec/mld_z_slu_solver.F90 | 2 +- mlprec/mld_z_sludist_solver.F90 | 57 +++--- mlprec/mld_z_umf_solver.F90 | 2 +- 85 files changed, 371 insertions(+), 533 deletions(-) diff --git a/README b/README index 827496d7..fcc0e1e2 100644 --- a/README +++ b/README @@ -12,7 +12,7 @@ Version 2.0. addition (and in front of) libpsb_prec.a, and no longer in place of it. 3. As for the basic usage, this is practically identical to the - previous version(s). + previous version(s). The Fortran 2003 support means that it is far easier to develop and integrate new solvers and smoothers; you need to take one of the @@ -20,8 +20,8 @@ Version 2.0. changing/replacing the model contents, and then pass the new object to the PREC%SET() method which will copy into the internals, as per the PROTOTYPE design pattern. It's easier done than said! - Note in this beta version SuperLU_Dist support is broken, we'll fix - it asap. + Note in this beta version SuperLU and SuperLU_Dist are not + working completely, we'll fix them asap. In version 1.1: diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 36e77d71..eeb9d199 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -316,7 +316,7 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: p + type(mld_cprec_type), intent(inout) :: p complex(psb_spk_),intent(in) :: alpha,beta complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) @@ -410,7 +410,7 @@ contains ! Arguments integer(psb_ipk_) :: level - type(mld_cprec_type), intent(in) :: p + type(mld_cprec_type), intent(inout) :: p type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) character, intent(in) :: trans complex(psb_spk_),target :: work(:) diff --git a/mlprec/impl/mld_cprecaply.f90 b/mlprec/impl/mld_cprecaply.f90 index c09fed03..eb2eda33 100644 --- a/mlprec/impl/mld_cprecaply.f90 +++ b/mlprec/impl/mld_cprecaply.f90 @@ -80,7 +80,7 @@ subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec + type(mld_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -215,7 +215,7 @@ subroutine mld_cprecaply1(prec,x,desc_data,info,trans) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec + type(mld_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/impl/mld_cslud_interface.c b/mlprec/impl/mld_cslud_interface.c index 40feae95..0268d67b 100644 --- a/mlprec/impl/mld_cslud_interface.c +++ b/mlprec/impl/mld_cslud_interface.c @@ -36,9 +36,9 @@ * POSSIBILITY OF SUCH DAMAGE. * * - * File: mld_zslud_interface.c + * File: mld_cslud_interface.c * - * Functions: mld_zsludist_fact_, mld_zsludist_solve_, mld_zsludist_free_. + * Functions: mld_csludist_fact, mld_csludist_solve, mld_csludist_free. * * This file is an interface to the SuperLU_dist routines for sparse factorization and * solve. It was obtained by modifying the c_fortran_zgssv.c file from the SuperLU_dist @@ -87,23 +87,17 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * */ -/* No single complex version in SuperLU_Dist */ -#ifdef Have_SLUDist_ -#undef Have_SLUDist_ +/* as of v 3.3 SLUDist does not have a single precision interface */ +#ifdef Have_SLUDist_ +#undef Have_SLUDist_ #endif + #ifdef Have_SLUDist_ #include #include "superlu_zdefs.h" #define HANDLE_SIZE 8 -/* kind of integer to hold a pointer. Use int. - This might need to be changed on 64-bit systems. */ -#ifdef Ptr64Bits -typedef long long fptr; -#else -typedef int fptr; /* 32-bit by default */ -#endif typedef struct { SuperMatrix *A; @@ -120,58 +114,22 @@ typedef struct { #endif -#ifdef LowerUnderscore -#define mld_csludist_fact_ mld_csludist_fact_ -#define mld_csludist_solve_ mld_csludist_solve_ -#define mld_csludist_free_ mld_csludist_free_ -#endif -#ifdef LowerDoubleUnderscore -#define mld_csludist_fact_ mld_csludist_fact__ -#define mld_csludist_solve_ mld_csludist_solve__ -#define mld_csludist_free_ mld_csludist_free__ -#endif -#ifdef LowerCase -#define mld_csludist_fact_ mld_csludist_fact -#define mld_csludist_solve_ mld_csludist_solve -#define mld_csludist_free_ mld_csludist_free -#endif -#ifdef UpperUnderscore -#define mld_csludist_fact_ MLD_CSLUDIST_FACT_ -#define mld_csludist_solve_ MLD_CSLUDIST_SOLVE_ -#define mld_csludist_free_ MLD_CSLUDIST_FREE_ -#endif -#ifdef UpperDoubleUnderscore -#define mld_csludist_fact_ MLD_CSLUDIST_FACT__ -#define mld_csludist_solve_ MLD_CSLUDIST_SOLVE__ -#define mld_csludist_free_ MLD_CSLUDIST_FREE__ -#endif -#ifdef UpperCase -#define mld_csludist_fact_ MLD_CSLUDIST_FACT -#define mld_csludist_solve_ MLD_CSLUDIST_SOLVE -#define mld_csludist_free_ MLD_CSLUDIST_FREE -#endif - - - - -void -mld_csludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, +int mld_csludist_fact(int n, int nl, int nnzl, int ffstr, #ifdef Have_SLUDist_ - complex *values, int *rowptr, int *colind, - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ + complex *values, int *rowptr, int *colind, + void **f_factors, #else - void *values, int *rowptr, int *colind, - void *f_factors, + void *values, int *rowptr, int *colind, + void **f_factors, #endif - int *nprow, int *npcol, int *info) - + int nprow, int npcol) + { /* * This routine can be called from Fortran. * performs LU decomposition. * - * f_factors (input/output) fptr* + * f_factors (input/output) void** * On output contains the pointer pointing to * the structure of the factored matrices. * @@ -185,7 +143,7 @@ mld_csludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; float drop_tol = 0.0,berr[1]; mem_usage_t mem_usage; @@ -198,42 +156,35 @@ mld_csludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, trans = NOTRANS; grid = (gridinfo_t *) SUPERLU_MALLOC(sizeof(gridinfo_t)); - superlu_gridinit(MPI_COMM_WORLD, *nprow, *npcol, grid); + superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, grid); /* Initialize the statistics variables. */ PStatInit(&stat); - fst_row = (*ffstr) -1; - /* Adjust to 0-based indexing */ - icol = (int *) malloc((*nnzl)*sizeof(int)); - irpt = (int *) malloc(((*nl)+1)*sizeof(int)); - ival = (complex *) malloc((*nnzl)*sizeof(doublecomplex)); - for (i = 0; i < *nnzl; ++i) ival[i] = values[i]; - for (i = 0; i < *nnzl; ++i) icol[i] = colind[i] -1; - for (i = 0; i <= *nl; ++i) irpt[i] = rowptr[i] -1; + fst_row = (ffstr); A = (SuperMatrix *) malloc(sizeof(SuperMatrix)); - zCreate_CompRowLoc_Matrix_dist(A, *n, *n, *nnzl, *nl, fst_row, - ival, icol, irpt, + zCreate_CompRowLoc_Matrix_dist(A, n, n, nnzl, nl, fst_row, + values, colind, rowptr, SLU_NR_loc, SLU_Z, SLU_GE); /* Initialize ScalePermstruct and LUstruct. */ ScalePermstruct = (ScalePermstruct_t *) SUPERLU_MALLOC(sizeof(ScalePermstruct_t)); LUstruct = (LUstruct_t *) SUPERLU_MALLOC(sizeof(LUstruct_t)); - ScalePermstructInit(*n,*n, ScalePermstruct); - LUstructInit(*n,*n, LUstruct); + ScalePermstructInit(n,n, ScalePermstruct); + LUstructInit(n,n, LUstruct); /* Set the default input options. */ set_default_options_dist(&options); options.IterRefine=NO; options.PrintStat=NO; - pzgssvx(&options, A, ScalePermstruct, b, *nl, 0, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pzgssvx(&options, A, ScalePermstruct, b, nl, 0, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); - if ( *info == 0 ) { + if ( info == 0 ) { ; } else { - printf("pzgssvx() error returns INFO= %d\n", *info); - if ( *info <= *n ) { /* factorization completes */ + printf("pzgssvx() error returns INFO= %d\n", info); + if ( info <= n ) { /* factorization completes */ ; } } @@ -252,28 +203,24 @@ mld_csludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, /* fprintf(stderr,"slud factor: A %p %p\n",A,LUfactors->A); */ /* fprintf(stderr,"slud factor: grid %p %p\n",grid,LUfactors->grid); */ /* fprintf(stderr,"slud factor: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - *f_factors = (fptr) LUfactors; - + *f_factors = (void *) LUfactors; PStatFree(&stat); + return(info); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } -void -mld_csludist_solve_(int *itrans, int *n, int *nrhs, +int mld_csludist_solve(int itrans, int n, int nrhs, #ifdef Have_SLUDist_ - doublecomplex *b, int *ldb, - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ + complex *b, #else - void *b, int *ldb, - void *f_factors, + void *b, #endif - int *info) - + int ldb, void *f_factors) + { /* * This routine can be called from Fortran. @@ -286,16 +233,16 @@ mld_csludist_solve_(int *itrans, int *n, int *nrhs, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; - double drop_tol = 0.0; - double *berr; + float drop_tol = 0.0; + float *berr; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; @@ -307,18 +254,18 @@ mld_csludist_solve_(int *itrans, int *n, int *nrhs, /* fprintf(stderr,"slud solve: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - if (*itrans == 0) { + if (itrans == 0) { trans = NOTRANS; - } else if (*itrans ==1) { + } else if (itrans ==1) { trans = TRANS; - } else if (*itrans ==2) { + } else if (itrans ==2) { trans = CONJ; } else { trans = NOTRANS; } /* fprintf(stderr,"Entry to sludist_solve\n"); */ - berr = (double *) malloc((*nrhs) *sizeof(double)); + berr = (float *) malloc((nrhs) *sizeof(float)); /* Initialize the statistics variables. */ PStatInit(&stat); @@ -329,33 +276,25 @@ mld_csludist_solve_(int *itrans, int *n, int *nrhs, options.Fact = FACTORED; options.PrintStat = NO; - pzgssvx(&options, A, ScalePermstruct, b, *ldb, *nrhs, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pzgssvx(&options, A, ScalePermstruct, b, ldb, nrhs, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); -/* fprintf(stderr,"Double check: after solve %d %lf\n",*info,berr[0]); */ +/* fprintf(stderr,"Float check: after solve %d %lf\n",*info,berr[0]); */ if (options.SolveInitialized) { zSolveFinalize(&options,&SOLVEstruct); } PStatFree(&stat); free(berr); + return(info); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } -void -mld_csludist_free_( -#ifdef Have_SLUDist_ - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ -#else - void *f_factors, -#endif - int *info) - +int mld_csludist_free(void *f_factors) { /* * This routine can be called from Fortran. @@ -371,20 +310,26 @@ mld_csludist_free_( gridinfo_t *grid; int i, panel_size, permc_spec, relax; trans_t trans; - double drop_tol = 0.0; - double *berr; + float drop_tol = 0.0; + float *berr; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + + if (f_factors == NULL) + return(0); + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; - Destroy_CompRowLoc_Matrix_dist(A); + // Memory leak: with SuperLU_Dist 3.3 + // we either have a leak or a segfault here. + // To be investigated further. + //Destroy_CompRowLoc_Matrix_dist(A); ScalePermstructFree(ScalePermstruct); LUstructFree(LUstruct); superlu_gridexit(grid); @@ -392,10 +337,11 @@ mld_csludist_free_( free(grid); free(LUstruct); free(LUfactors); + return(0); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index bac8406a..5333f556 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -316,7 +316,7 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: p + type(mld_dprec_type), intent(inout) :: p real(psb_dpk_),intent(in) :: alpha,beta real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) @@ -410,7 +410,7 @@ contains ! Arguments integer(psb_ipk_) :: level - type(mld_dprec_type), intent(in) :: p + type(mld_dprec_type), intent(inout) :: p type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) character, intent(in) :: trans real(psb_dpk_),target :: work(:) diff --git a/mlprec/impl/mld_dprecaply.f90 b/mlprec/impl/mld_dprecaply.f90 index 5a3d0059..4b679425 100644 --- a/mlprec/impl/mld_dprecaply.f90 +++ b/mlprec/impl/mld_dprecaply.f90 @@ -80,7 +80,7 @@ subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec + type(mld_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -215,7 +215,7 @@ subroutine mld_dprecaply1(prec,x,desc_data,info,trans) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec + type(mld_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/impl/mld_dslud_interface.c b/mlprec/impl/mld_dslud_interface.c index 074acbdc..e5749e0e 100644 --- a/mlprec/impl/mld_dslud_interface.c +++ b/mlprec/impl/mld_dslud_interface.c @@ -36,7 +36,7 @@ * POSSIBILITY OF SUCH DAMAGE. * * - * File: mld_slud_interface.c + * File: mld_dslud_interface.c * * Functions: mld_dsludist_fact, mld_dsludist_solve, mld_dsludist_free. * @@ -142,7 +142,6 @@ int mld_dsludist_fact(int n, int nl, int nnzl, int ffstr, double *ival; trans = NOTRANS; -/* fprintf(stderr,"Entry to sludist_fact\n"); */ grid = (gridinfo_t *) SUPERLU_MALLOC(sizeof(gridinfo_t)); superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, grid); /* Initialize the statistics variables. */ @@ -202,7 +201,7 @@ int mld_dsludist_fact(int n, int nl, int nnzl, int ffstr, int mld_dsludist_solve(int itrans, int n, int nrhs, - double *b, int ldb, void *f_factors) + double *b, int ldb, void *f_factors) { /* @@ -231,7 +230,6 @@ int mld_dsludist_solve(int itrans, int n, int nrhs, grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; - fprintf(stderr,"slud solve: ldb %d n %d \n",ldb,n); /* fprintf(stderr,"slud solve: LUFactors %p \n",LUfactors); */ /* fprintf(stderr,"slud solve: A %p %p\n",A,LUfactors->A); */ /* fprintf(stderr,"slud solve: grid %p %p\n",grid,LUfactors->grid); */ @@ -279,8 +277,6 @@ int mld_dsludist_solve(int itrans, int n, int nrhs, int mld_dsludist_free(void *f_factors) - - { /* * This routine can be called from Fortran. @@ -312,7 +308,10 @@ int mld_dsludist_free(void *f_factors) grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; - Destroy_CompRowLoc_Matrix_dist(A); + // Memory leak: with SuperLU_Dist 3.3 + // we either have a leak or a segfault here. + // To be investigated further. + //Destroy_CompRowLoc_Matrix_dist(A); ScalePermstructFree(ScalePermstruct); LUstructFree(LUstruct); superlu_gridexit(grid); diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 1a84c9d6..b5ef086e 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -316,7 +316,7 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: p + type(mld_sprec_type), intent(inout) :: p real(psb_spk_),intent(in) :: alpha,beta real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) @@ -410,7 +410,7 @@ contains ! Arguments integer(psb_ipk_) :: level - type(mld_sprec_type), intent(in) :: p + type(mld_sprec_type), intent(inout) :: p type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) character, intent(in) :: trans real(psb_spk_),target :: work(:) diff --git a/mlprec/impl/mld_sprecaply.f90 b/mlprec/impl/mld_sprecaply.f90 index 8068e4be..fa079c35 100644 --- a/mlprec/impl/mld_sprecaply.f90 +++ b/mlprec/impl/mld_sprecaply.f90 @@ -80,7 +80,7 @@ subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec + type(mld_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -215,7 +215,7 @@ subroutine mld_sprecaply1(prec,x,desc_data,info,trans) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec + type(mld_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/impl/mld_sslud_interface.c b/mlprec/impl/mld_sslud_interface.c index 0b98301b..ff1a1491 100644 --- a/mlprec/impl/mld_sslud_interface.c +++ b/mlprec/impl/mld_sslud_interface.c @@ -38,7 +38,7 @@ * * File: mld_slud_interface.c * - * Functions: mld_ssludist_fact_, mld_ssludist_solve_, mld_ssludist_free_. + * Functions: mld_ssludist_fact, mld_ssludist_solve, mld_ssludist_free. * * This file is an interface to the SuperLU_dist routines for sparse factorization and * solve. It was obtained by modifying the c_fortran_dgssv.c file from the SuperLU_dist @@ -87,23 +87,16 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. * */ -/* as of v 2.1 SLUDist does not have a single precision interface */ +/* as of v 3.3 SLUDist does not have a single precision interface */ #ifdef Have_SLUDist_ #undef Have_SLUDist_ #endif -#ifdef Have_SLUDist_ +#ifdef Have_SLUDist_ #include #include "superlu_sdefs.h" #define HANDLE_SIZE 8 -/* kind of integer to hold a pointer. Use int. - This might need to be changed on 64-bit systems. */ -#ifdef Ptr64Bits -typedef long long fptr; -#else -typedef int fptr; /* 32-bit by default */ -#endif typedef struct { SuperMatrix *A; @@ -120,57 +113,15 @@ typedef struct { #endif -#ifdef LowerUnderscore -#define mld_ssludist_fact_ mld_ssludist_fact_ -#define mld_ssludist_solve_ mld_ssludist_solve_ -#define mld_ssludist_free_ mld_ssludist_free_ -#endif -#ifdef LowerDoubleUnderscore -#define mld_ssludist_fact_ mld_ssludist_fact__ -#define mld_ssludist_solve_ mld_ssludist_solve__ -#define mld_ssludist_free_ mld_ssludist_free__ -#endif -#ifdef LowerCase -#define mld_ssludist_fact_ mld_ssludist_fact -#define mld_ssludist_solve_ mld_ssludist_solve -#define mld_ssludist_free_ mld_ssludist_free -#endif -#ifdef UpperUnderscore -#define mld_ssludist_fact_ MLD_SSLUDIST_FACT_ -#define mld_ssludist_solve_ MLD_SSLUDIST_SOLVE_ -#define mld_ssludist_free_ MLD_SSLUDIST_FREE_ -#endif -#ifdef UpperDoubleUnderscore -#define mld_ssludist_fact_ MLD_SSLUDIST_FACT__ -#define mld_ssludist_solve_ MLD_SSLUDIST_SOLVE__ -#define mld_ssludist_free_ MLD_SSLUDIST_FREE__ -#endif -#ifdef UpperCase -#define mld_ssludist_fact_ MLD_SSLUDIST_FACT -#define mld_ssludist_solve_ MLD_SSLUDIST_SOLVE -#define mld_ssludist_free_ MLD_SSLUDIST_FREE -#endif - - - - -void -mld_ssludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, - double *values, int *rowptr, int *colind, -#ifdef Have_SLUDist_ - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ -#else - void *f_factors, -#endif - int *nprow, int *npcol, int *info) - +int mld_ssludist_fact(int n, int nl, int nnzl, int ffstr, + float *values, int *rowptr, int *colind, + void **f_factors, int nprow, int npcol) { /* * This routine can be called from Fortran. * performs LU decomposition. * - * f_factors (input/output) fptr* + * f_factors (input/output) void** * On output contains the pointer pointing to * the structure of the factored matrices. * @@ -184,56 +135,48 @@ mld_ssludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; - double drop_tol = 0.0,b[1],berr[1]; + float drop_tol = 0.0, b[1], berr[1]; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; factors_t *LUfactors; int fst_row; int *icol,*irpt; - double *ival; + float *ival; trans = NOTRANS; -/* fprintf(stderr,"Entry to sludist_fact\n"); */ grid = (gridinfo_t *) SUPERLU_MALLOC(sizeof(gridinfo_t)); - superlu_gridinit(MPI_COMM_WORLD, *nprow, *npcol, grid); + superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, grid); /* Initialize the statistics variables. */ PStatInit(&stat); - fst_row = (*ffstr) -1; - /* Adjust to 0-based indexing */ - icol = (int *) malloc((*nnzl)*sizeof(int)); - irpt = (int *) malloc(((*nl)+1)*sizeof(int)); - ival = (double *) malloc((*nnzl)*sizeof(double)); - for (i = 0; i < *nnzl; ++i) ival[i] = values[i]; - for (i = 0; i < *nnzl; ++i) icol[i] = colind[i] -1; - for (i = 0; i <= *nl; ++i) irpt[i] = rowptr[i] -1; + fst_row = (ffstr); A = (SuperMatrix *) malloc(sizeof(SuperMatrix)); - dCreate_CompRowLoc_Matrix_dist(A, *n, *n, *nnzl, *nl, fst_row, - ival, icol, irpt, + dCreate_CompRowLoc_Matrix_dist(A, n, n, nnzl, nl, fst_row, + values, colind, rowptr, SLU_NR_loc, SLU_D, SLU_GE); /* Initialize ScalePermstruct and LUstruct. */ ScalePermstruct = (ScalePermstruct_t *) SUPERLU_MALLOC(sizeof(ScalePermstruct_t)); LUstruct = (LUstruct_t *) SUPERLU_MALLOC(sizeof(LUstruct_t)); - ScalePermstructInit(*n,*n, ScalePermstruct); - LUstructInit(*n,*n, LUstruct); + ScalePermstructInit(n,n, ScalePermstruct); + LUstructInit(n,n, LUstruct); /* Set the default input options. */ set_default_options_dist(&options); options.IterRefine=NO; options.PrintStat=NO; - pdgssvx(&options, A, ScalePermstruct, b, *nl, 0, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pdgssvx(&options, A, ScalePermstruct, b, nl, 0, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); - if ( *info == 0 ) { + if ( info == 0 ) { ; } else { - printf("pdgssvx() error returns INFO= %d\n", *info); - if ( *info <= *n ) { /* factorization completes */ + printf("pdgssvx() error returns INFO= %d\n", info); + if ( info <= n ) { /* factorization completes */ ; } } @@ -252,26 +195,18 @@ mld_ssludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, /* fprintf(stderr,"slud factor: A %p %p\n",A,LUfactors->A); */ /* fprintf(stderr,"slud factor: grid %p %p\n",grid,LUfactors->grid); */ /* fprintf(stderr,"slud factor: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - *f_factors = (fptr) LUfactors; - + *f_factors = (void *) LUfactors; PStatFree(&stat); + return(info); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } -void -mld_ssludist_solve_(int *itrans, int *n, int *nrhs, - double *b, int *ldb, -#ifdef Have_SLUDist_ - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ -#else - void *f_factors, -#endif - int *info) +int mld_ssludist_solve(int itrans, int n, int nrhs, + float *b, int ldb, void *f_factors) { /* @@ -285,39 +220,40 @@ mld_ssludist_solve_(int *itrans, int *n, int *nrhs, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; - double drop_tol = 0.0; - double *berr; + float drop_tol = 0.0; + float *berr; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; + fprintf(stderr,"slud solve: ldb %d n %d \n",ldb,n); /* fprintf(stderr,"slud solve: LUFactors %p \n",LUfactors); */ /* fprintf(stderr,"slud solve: A %p %p\n",A,LUfactors->A); */ /* fprintf(stderr,"slud solve: grid %p %p\n",grid,LUfactors->grid); */ /* fprintf(stderr,"slud solve: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - if (*itrans == 0) { + if (itrans == 0) { trans = NOTRANS; - } else if (*itrans ==1) { + } else if (itrans ==1) { trans = TRANS; - } else if (*itrans ==2) { + } else if (itrans ==2) { trans = CONJ; } else { trans = NOTRANS; } /* fprintf(stderr,"Entry to sludist_solve\n"); */ - berr = (double *) malloc((*nrhs) *sizeof(double)); + berr = (float *) malloc((nrhs) *sizeof(float)); /* Initialize the statistics variables. */ PStatInit(&stat); @@ -328,32 +264,26 @@ mld_ssludist_solve_(int *itrans, int *n, int *nrhs, options.Fact = FACTORED; options.PrintStat = NO; - pdgssvx(&options, A, ScalePermstruct, b, *ldb, *nrhs, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pdgssvx(&options, A, ScalePermstruct, b, ldb, nrhs, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); -/* fprintf(stderr,"Double check: after solve %d %lf\n",*info,berr[0]); */ +/* fprintf(stderr,"Float check: after solve %d %lf\n",*info,berr[0]); */ if (options.SolveInitialized) { dSolveFinalize(&options,&SOLVEstruct); } PStatFree(&stat); free(berr); + return(info); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } -void -mld_ssludist_free_( -#ifdef Have_SLUDist_ - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ -#else - void *f_factors, -#endif - int *info) +int mld_ssludist_free(void *f_factors) + { /* @@ -370,20 +300,26 @@ mld_ssludist_free_( gridinfo_t *grid; int i, panel_size, permc_spec, relax; trans_t trans; - double drop_tol = 0.0; - double *berr; + float drop_tol = 0.0; + float *berr; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + + if (f_factors == NULL) + return(0); + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; - Destroy_CompRowLoc_Matrix_dist(A); + // Memory leak: with SuperLU_Dist 3.3 + // we either have a leak or a segfault here. + // To be investigated further. + //Destroy_CompRowLoc_Matrix_dist(A); ScalePermstructFree(ScalePermstruct); LUstructFree(LUstruct); superlu_gridexit(grid); @@ -391,11 +327,11 @@ mld_ssludist_free_( free(grid); free(LUstruct); free(LUfactors); + return(0); #else - fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + fprintf(stderr," SLUDist does not have single precision, sorry.\n"); + return(-1); #endif } - diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index cd43a381..a02962f6 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -316,7 +316,7 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: p + type(mld_zprec_type), intent(inout) :: p complex(psb_dpk_),intent(in) :: alpha,beta complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) @@ -410,7 +410,7 @@ contains ! Arguments integer(psb_ipk_) :: level - type(mld_zprec_type), intent(in) :: p + type(mld_zprec_type), intent(inout) :: p type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) character, intent(in) :: trans complex(psb_dpk_),target :: work(:) diff --git a/mlprec/impl/mld_zprecaply.f90 b/mlprec/impl/mld_zprecaply.f90 index b5219d00..5067d546 100644 --- a/mlprec/impl/mld_zprecaply.f90 +++ b/mlprec/impl/mld_zprecaply.f90 @@ -80,7 +80,7 @@ subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec + type(mld_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -215,7 +215,7 @@ subroutine mld_zprecaply1(prec,x,desc_data,info,trans) ! Arguments type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec + type(mld_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/impl/mld_zslud_interface.c b/mlprec/impl/mld_zslud_interface.c index 03c095d2..664f6400 100644 --- a/mlprec/impl/mld_zslud_interface.c +++ b/mlprec/impl/mld_zslud_interface.c @@ -38,7 +38,7 @@ * * File: mld_zslud_interface.c * - * Functions: mld_zsludist_fact_, mld_zsludist_solve_, mld_zsludist_free_. + * Functions: mld_zsludist_fact, mld_zsludist_solve, mld_zsludist_free. * * This file is an interface to the SuperLU_dist routines for sparse factorization and * solve. It was obtained by modifying the c_fortran_zgssv.c file from the SuperLU_dist @@ -92,13 +92,6 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "superlu_zdefs.h" #define HANDLE_SIZE 8 -/* kind of integer to hold a pointer. Use int. - This might need to be changed on 64-bit systems. */ -#ifdef Ptr64Bits -typedef long long fptr; -#else -typedef int fptr; /* 32-bit by default */ -#endif typedef struct { SuperMatrix *A; @@ -115,58 +108,22 @@ typedef struct { #endif -#ifdef LowerUnderscore -#define mld_zsludist_fact_ mld_zsludist_fact_ -#define mld_zsludist_solve_ mld_zsludist_solve_ -#define mld_zsludist_free_ mld_zsludist_free_ -#endif -#ifdef LowerDoubleUnderscore -#define mld_zsludist_fact_ mld_zsludist_fact__ -#define mld_zsludist_solve_ mld_zsludist_solve__ -#define mld_zsludist_free_ mld_zsludist_free__ -#endif -#ifdef LowerCase -#define mld_zsludist_fact_ mld_zsludist_fact -#define mld_zsludist_solve_ mld_zsludist_solve -#define mld_zsludist_free_ mld_zsludist_free -#endif -#ifdef UpperUnderscore -#define mld_zsludist_fact_ MLD_ZSLUDIST_FACT_ -#define mld_zsludist_solve_ MLD_ZSLUDIST_SOLVE_ -#define mld_zsludist_free_ MLD_ZSLUDIST_FREE_ -#endif -#ifdef UpperDoubleUnderscore -#define mld_zsludist_fact_ MLD_ZSLUDIST_FACT__ -#define mld_zsludist_solve_ MLD_ZSLUDIST_SOLVE__ -#define mld_zsludist_free_ MLD_ZSLUDIST_FREE__ -#endif -#ifdef UpperCase -#define mld_zsludist_fact_ MLD_ZSLUDIST_FACT -#define mld_zsludist_solve_ MLD_ZSLUDIST_SOLVE -#define mld_zsludist_free_ MLD_ZSLUDIST_FREE -#endif - - - - -void -mld_zsludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, +int mld_zsludist_fact(int n, int nl, int nnzl, int ffstr, #ifdef Have_SLUDist_ - doublecomplex *values, int *rowptr, int *colind, - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ + doublecomplex *values, int *rowptr, int *colind, + void **f_factors, #else - void *values, int *rowptr, int *colind, - void *f_factors, + void *values, int *rowptr, int *colind, + void **f_factors, #endif - int *nprow, int *npcol, int *info) - + int nprow, int npcol) + { /* * This routine can be called from Fortran. * performs LU decomposition. * - * f_factors (input/output) fptr* + * f_factors (input/output) void** * On output contains the pointer pointing to * the structure of the factored matrices. * @@ -180,7 +137,7 @@ mld_zsludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; double drop_tol = 0.0,berr[1]; mem_usage_t mem_usage; @@ -193,42 +150,35 @@ mld_zsludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, trans = NOTRANS; grid = (gridinfo_t *) SUPERLU_MALLOC(sizeof(gridinfo_t)); - superlu_gridinit(MPI_COMM_WORLD, *nprow, *npcol, grid); + superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, grid); /* Initialize the statistics variables. */ PStatInit(&stat); - fst_row = (*ffstr) -1; - /* Adjust to 0-based indexing */ - icol = (int *) malloc((*nnzl)*sizeof(int)); - irpt = (int *) malloc(((*nl)+1)*sizeof(int)); - ival = (doublecomplex *) malloc((*nnzl)*sizeof(doublecomplex)); - for (i = 0; i < *nnzl; ++i) ival[i] = values[i]; - for (i = 0; i < *nnzl; ++i) icol[i] = colind[i] -1; - for (i = 0; i <= *nl; ++i) irpt[i] = rowptr[i] -1; + fst_row = (ffstr); A = (SuperMatrix *) malloc(sizeof(SuperMatrix)); - zCreate_CompRowLoc_Matrix_dist(A, *n, *n, *nnzl, *nl, fst_row, - ival, icol, irpt, + zCreate_CompRowLoc_Matrix_dist(A, n, n, nnzl, nl, fst_row, + values, colind, rowptr, SLU_NR_loc, SLU_Z, SLU_GE); /* Initialize ScalePermstruct and LUstruct. */ ScalePermstruct = (ScalePermstruct_t *) SUPERLU_MALLOC(sizeof(ScalePermstruct_t)); LUstruct = (LUstruct_t *) SUPERLU_MALLOC(sizeof(LUstruct_t)); - ScalePermstructInit(*n,*n, ScalePermstruct); - LUstructInit(*n,*n, LUstruct); + ScalePermstructInit(n,n, ScalePermstruct); + LUstructInit(n,n, LUstruct); /* Set the default input options. */ set_default_options_dist(&options); options.IterRefine=NO; options.PrintStat=NO; - pzgssvx(&options, A, ScalePermstruct, b, *nl, 0, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pzgssvx(&options, A, ScalePermstruct, b, nl, 0, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); - if ( *info == 0 ) { + if ( info == 0 ) { ; } else { - printf("pzgssvx() error returns INFO= %d\n", *info); - if ( *info <= *n ) { /* factorization completes */ + printf("pzgssvx() error returns INFO= %d\n", info); + if ( info <= n ) { /* factorization completes */ ; } } @@ -247,28 +197,24 @@ mld_zsludist_fact_(int *n, int *nl, int *nnzl, int *ffstr, /* fprintf(stderr,"slud factor: A %p %p\n",A,LUfactors->A); */ /* fprintf(stderr,"slud factor: grid %p %p\n",grid,LUfactors->grid); */ /* fprintf(stderr,"slud factor: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - *f_factors = (fptr) LUfactors; - + *f_factors = (void *) LUfactors; PStatFree(&stat); + return(info); #else fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + return(-1); #endif } -void -mld_zsludist_solve_(int *itrans, int *n, int *nrhs, +int mld_zsludist_solve(int itrans, int n, int nrhs, #ifdef Have_SLUDist_ - doublecomplex *b, int *ldb, - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ + doublecomplex *b, #else - void *b, int *ldb, - void *f_factors, + void *b, #endif - int *info) - + int ldb, void *f_factors) + { /* * This routine can be called from Fortran. @@ -281,7 +227,7 @@ mld_zsludist_solve_(int *itrans, int *n, int *nrhs, LUstruct_t *LUstruct; SOLVEstruct_t SOLVEstruct; gridinfo_t *grid; - int i, panel_size, permc_spec, relax; + int i, panel_size, permc_spec, relax, info; trans_t trans; double drop_tol = 0.0; double *berr; @@ -290,7 +236,7 @@ mld_zsludist_solve_(int *itrans, int *n, int *nrhs, SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; @@ -302,18 +248,18 @@ mld_zsludist_solve_(int *itrans, int *n, int *nrhs, /* fprintf(stderr,"slud solve: LUstruct %p %p\n",LUstruct,LUfactors->LUstruct); */ - if (*itrans == 0) { + if (itrans == 0) { trans = NOTRANS; - } else if (*itrans ==1) { + } else if (itrans ==1) { trans = TRANS; - } else if (*itrans ==2) { + } else if (itrans ==2) { trans = CONJ; } else { trans = NOTRANS; } /* fprintf(stderr,"Entry to sludist_solve\n"); */ - berr = (double *) malloc((*nrhs) *sizeof(double)); + berr = (double *) malloc((nrhs) *sizeof(double)); /* Initialize the statistics variables. */ PStatInit(&stat); @@ -324,8 +270,8 @@ mld_zsludist_solve_(int *itrans, int *n, int *nrhs, options.Fact = FACTORED; options.PrintStat = NO; - pzgssvx(&options, A, ScalePermstruct, b, *ldb, *nrhs, - grid, LUstruct, &SOLVEstruct, berr, &stat, info); + pzgssvx(&options, A, ScalePermstruct, b, ldb, nrhs, + grid, LUstruct, &SOLVEstruct, berr, &stat, &info); /* fprintf(stderr,"Double check: after solve %d %lf\n",*info,berr[0]); */ if (options.SolveInitialized) { @@ -333,24 +279,16 @@ mld_zsludist_solve_(int *itrans, int *n, int *nrhs, } PStatFree(&stat); free(berr); + return(info); #else fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + return(-1); #endif } -void -mld_zsludist_free_( -#ifdef Have_SLUDist_ - fptr *f_factors, /* a handle containing the address - pointing to the factored matrices */ -#else - void *f_factors, -#endif - int *info) - +int mld_zsludist_free(void *f_factors) { /* * This routine can be called from Fortran. @@ -373,13 +311,19 @@ mld_zsludist_free_( SuperLUStat_t stat; factors_t *LUfactors; - LUfactors = (factors_t *) *f_factors ; + + if (f_factors == NULL) + return(0); + LUfactors = (factors_t *) f_factors ; A = LUfactors->A ; LUstruct = LUfactors->LUstruct ; grid = LUfactors->grid ; ScalePermstruct = LUfactors->ScalePermstruct; - Destroy_CompRowLoc_Matrix_dist(A); + // Memory leak: with SuperLU_Dist 3.3 + // we either have a leak or a segfault here. + // To be investigated further. + //Destroy_CompRowLoc_Matrix_dist(A); ScalePermstructFree(ScalePermstruct); LUstructFree(LUstruct); superlu_gridexit(grid); @@ -387,10 +331,11 @@ mld_zsludist_free_( free(grid); free(LUstruct); free(LUfactors); + return(0); #else fprintf(stderr," SLUDist Not Configured, fix make.inc and recompile\n"); - *info=-1; + return(-1); #endif } diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 index 5c63236d..b1124cbf 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_as_smoother_type), intent(in) :: sm + class(mld_c_as_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_c_base_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_base_smoother_apply.f90 index 7ad50777..fabfa841 100644 --- a/mlprec/impl/smoother/mld_c_base_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_base_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo use mld_c_base_smoother_mod, mld_protect_name => mld_c_base_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_base_smoother_type), intent(in) :: sm + class(mld_c_base_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 index 2ab38aa3..c17a3c4d 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_jac_smoother_type), intent(in) :: sm + class(mld_c_jac_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 index 4d2379e7..2f1830a9 100644 --- a/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_as_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_as_smoother_type), intent(in) :: sm + class(mld_d_as_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_d_base_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_base_smoother_apply.f90 index f1654f25..d821e11a 100644 --- a/mlprec/impl/smoother/mld_d_base_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_base_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo use mld_d_base_smoother_mod, mld_protect_name => mld_d_base_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_base_smoother_type), intent(in) :: sm + class(mld_d_base_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 index 943614a6..1ed15e12 100644 --- a/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_d_jac_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_jac_smoother_type), intent(in) :: sm + class(mld_d_jac_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 index efc977a1..67987669 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_as_smoother_type), intent(in) :: sm + class(mld_s_as_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_s_base_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_base_smoother_apply.f90 index 23021001..a43e0c6f 100644 --- a/mlprec/impl/smoother/mld_s_base_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_base_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo use mld_s_base_smoother_mod, mld_protect_name => mld_s_base_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_base_smoother_type), intent(in) :: sm + class(mld_s_base_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 index 75ff99e6..b9017c55 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_jac_smoother_type), intent(in) :: sm + class(mld_s_jac_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 index 4eb4a8cb..4231d08a 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_as_smoother_type), intent(in) :: sm + class(mld_z_as_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_z_base_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_base_smoother_apply.f90 index 4fe675e4..d08e4e96 100644 --- a/mlprec/impl/smoother/mld_z_base_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_base_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wo use mld_z_base_smoother_mod, mld_protect_name => mld_z_base_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_base_smoother_type), intent(in) :: sm + class(mld_z_base_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 index 8d5f976c..295eedbf 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,wor use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_jac_smoother_type), intent(in) :: sm + class(mld_z_jac_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_c_base_solver_apply.f90 b/mlprec/impl/solver/mld_c_base_solver_apply.f90 index e32066ba..7ce49077 100644 --- a/mlprec/impl/solver/mld_c_base_solver_apply.f90 +++ b/mlprec/impl/solver/mld_c_base_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_c_base_solver_mod, mld_protect_name => mld_c_base_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_base_solver_type), intent(in) :: sv + class(mld_c_base_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_c_diag_solver_apply.f90 b/mlprec/impl/solver/mld_c_diag_solver_apply.f90 index 291ca984..a93b5392 100644 --- a/mlprec/impl/solver/mld_c_diag_solver_apply.f90 +++ b/mlprec/impl/solver/mld_c_diag_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_c_diag_solver, mld_protect_name => mld_c_diag_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_diag_solver_type), intent(in) :: sv + class(mld_c_diag_solver_type), intent(inout) :: sv complex(psb_spk_), intent(inout) :: x(:) complex(psb_spk_), intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_c_id_solver_apply.f90 b/mlprec/impl/solver/mld_c_id_solver_apply.f90 index 57429983..ce4c4cab 100644 --- a/mlprec/impl/solver/mld_c_id_solver_apply.f90 +++ b/mlprec/impl/solver/mld_c_id_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_c_id_solver, mld_protect_name => mld_c_id_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_id_solver_type), intent(in) :: sv + class(mld_c_id_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_c_ilu_solver_apply.f90 b/mlprec/impl/solver/mld_c_ilu_solver_apply.f90 index e2c1aa91..ea5b86de 100644 --- a/mlprec/impl/solver/mld_c_ilu_solver_apply.f90 +++ b/mlprec/impl/solver/mld_c_ilu_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_c_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_c_ilu_solver, mld_protect_name => mld_c_ilu_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_ilu_solver_type), intent(in) :: sv + class(mld_c_ilu_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_d_base_solver_apply.f90 b/mlprec/impl/solver/mld_d_base_solver_apply.f90 index 5fb7158b..7551be37 100644 --- a/mlprec/impl/solver/mld_d_base_solver_apply.f90 +++ b/mlprec/impl/solver/mld_d_base_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_d_base_solver_mod, mld_protect_name => mld_d_base_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_base_solver_type), intent(in) :: sv + class(mld_d_base_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_d_diag_solver_apply.f90 b/mlprec/impl/solver/mld_d_diag_solver_apply.f90 index 1889cbe1..32703cd3 100644 --- a/mlprec/impl/solver/mld_d_diag_solver_apply.f90 +++ b/mlprec/impl/solver/mld_d_diag_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_d_diag_solver, mld_protect_name => mld_d_diag_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_diag_solver_type), intent(in) :: sv + class(mld_d_diag_solver_type), intent(inout) :: sv real(psb_dpk_), intent(inout) :: x(:) real(psb_dpk_), intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_d_id_solver_apply.f90 b/mlprec/impl/solver/mld_d_id_solver_apply.f90 index 6fce6b55..2e2bb5f0 100644 --- a/mlprec/impl/solver/mld_d_id_solver_apply.f90 +++ b/mlprec/impl/solver/mld_d_id_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_d_id_solver, mld_protect_name => mld_d_id_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_id_solver_type), intent(in) :: sv + class(mld_d_id_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_d_ilu_solver_apply.f90 b/mlprec/impl/solver/mld_d_ilu_solver_apply.f90 index 6d089e64..dc3216bb 100644 --- a/mlprec/impl/solver/mld_d_ilu_solver_apply.f90 +++ b/mlprec/impl/solver/mld_d_ilu_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_d_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_d_ilu_solver, mld_protect_name => mld_d_ilu_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_ilu_solver_type), intent(in) :: sv + class(mld_d_ilu_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_s_base_solver_apply.f90 b/mlprec/impl/solver/mld_s_base_solver_apply.f90 index 131696f7..eac44b9b 100644 --- a/mlprec/impl/solver/mld_s_base_solver_apply.f90 +++ b/mlprec/impl/solver/mld_s_base_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_s_base_solver_mod, mld_protect_name => mld_s_base_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_base_solver_type), intent(in) :: sv + class(mld_s_base_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_s_diag_solver_apply.f90 b/mlprec/impl/solver/mld_s_diag_solver_apply.f90 index 2ca24cb8..7e11839b 100644 --- a/mlprec/impl/solver/mld_s_diag_solver_apply.f90 +++ b/mlprec/impl/solver/mld_s_diag_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_s_diag_solver, mld_protect_name => mld_s_diag_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_diag_solver_type), intent(in) :: sv + class(mld_s_diag_solver_type), intent(inout) :: sv real(psb_spk_), intent(inout) :: x(:) real(psb_spk_), intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_s_id_solver_apply.f90 b/mlprec/impl/solver/mld_s_id_solver_apply.f90 index e9c0a6e5..c2bef19d 100644 --- a/mlprec/impl/solver/mld_s_id_solver_apply.f90 +++ b/mlprec/impl/solver/mld_s_id_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_s_id_solver, mld_protect_name => mld_s_id_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_id_solver_type), intent(in) :: sv + class(mld_s_id_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_s_ilu_solver_apply.f90 b/mlprec/impl/solver/mld_s_ilu_solver_apply.f90 index 01b6afb6..b636ee32 100644 --- a/mlprec/impl/solver/mld_s_ilu_solver_apply.f90 +++ b/mlprec/impl/solver/mld_s_ilu_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_s_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_s_ilu_solver, mld_protect_name => mld_s_ilu_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_ilu_solver_type), intent(in) :: sv + class(mld_s_ilu_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_z_base_solver_apply.f90 b/mlprec/impl/solver/mld_z_base_solver_apply.f90 index 571a6718..59ad6751 100644 --- a/mlprec/impl/solver/mld_z_base_solver_apply.f90 +++ b/mlprec/impl/solver/mld_z_base_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_base_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_z_base_solver_mod, mld_protect_name => mld_z_base_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_base_solver_type), intent(in) :: sv + class(mld_z_base_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_z_diag_solver_apply.f90 b/mlprec/impl/solver/mld_z_diag_solver_apply.f90 index 57c471eb..c6c21a26 100644 --- a/mlprec/impl/solver/mld_z_diag_solver_apply.f90 +++ b/mlprec/impl/solver/mld_z_diag_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_diag_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_z_diag_solver, mld_protect_name => mld_z_diag_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_diag_solver_type), intent(in) :: sv + class(mld_z_diag_solver_type), intent(inout) :: sv complex(psb_dpk_), intent(inout) :: x(:) complex(psb_dpk_), intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_z_id_solver_apply.f90 b/mlprec/impl/solver/mld_z_id_solver_apply.f90 index 6f1001a3..a75497c3 100644 --- a/mlprec/impl/solver/mld_z_id_solver_apply.f90 +++ b/mlprec/impl/solver/mld_z_id_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_id_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_z_id_solver, mld_protect_name => mld_z_id_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_id_solver_type), intent(in) :: sv + class(mld_z_id_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/impl/solver/mld_z_ilu_solver_apply.f90 b/mlprec/impl/solver/mld_z_ilu_solver_apply.f90 index 285fd425..8ca1583f 100644 --- a/mlprec/impl/solver/mld_z_ilu_solver_apply.f90 +++ b/mlprec/impl/solver/mld_z_ilu_solver_apply.f90 @@ -42,7 +42,7 @@ subroutine mld_z_ilu_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) use mld_z_ilu_solver, mld_protect_name => mld_z_ilu_solver_apply implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_ilu_solver_type), intent(in) :: sv + class(mld_z_ilu_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_as_smoother.f90 b/mlprec/mld_c_as_smoother.f90 index 82d656a7..e1dcabc9 100644 --- a/mlprec/mld_c_as_smoother.f90 +++ b/mlprec/mld_c_as_smoother.f90 @@ -117,7 +117,7 @@ module mld_c_as_smoother & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, psb_desc_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_as_smoother_type), intent(in) :: sm + class(mld_c_as_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_base_smoother_mod.f90 b/mlprec/mld_c_base_smoother_mod.f90 index 464251c6..63ce0f88 100644 --- a/mlprec/mld_c_base_smoother_mod.f90 +++ b/mlprec/mld_c_base_smoother_mod.f90 @@ -127,7 +127,7 @@ module mld_c_base_smoother_mod & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & & mld_c_base_smoother_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_c_base_smoother_type), intent(in) :: sm + class(mld_c_base_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_base_solver_mod.f90 b/mlprec/mld_c_base_solver_mod.f90 index 76b2eb15..5762113c 100644 --- a/mlprec/mld_c_base_solver_mod.f90 +++ b/mlprec/mld_c_base_solver_mod.f90 @@ -118,7 +118,7 @@ module mld_c_base_solver_mod & mld_c_base_solver_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_base_solver_type), intent(in) :: sv + class(mld_c_base_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_diag_solver.f90 b/mlprec/mld_c_diag_solver.f90 index 99dd173b..f0e551b1 100644 --- a/mlprec/mld_c_diag_solver.f90 +++ b/mlprec/mld_c_diag_solver.f90 @@ -90,7 +90,7 @@ module mld_c_diag_solver & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & & mld_c_diag_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_c_diag_solver_type), intent(in) :: sv + class(mld_c_diag_solver_type), intent(inout) :: sv complex(psb_spk_), intent(inout) :: x(:) complex(psb_spk_), intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_id_solver.f90 b/mlprec/mld_c_id_solver.f90 index 0df6cb41..8b3e1fd0 100644 --- a/mlprec/mld_c_id_solver.f90 +++ b/mlprec/mld_c_id_solver.f90 @@ -85,7 +85,7 @@ module mld_c_id_solver & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & & mld_c_id_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_c_id_solver_type), intent(in) :: sv + class(mld_c_id_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_ilu_solver.f90 b/mlprec/mld_c_ilu_solver.f90 index 70cbf638..5018cdc8 100644 --- a/mlprec/mld_c_ilu_solver.f90 +++ b/mlprec/mld_c_ilu_solver.f90 @@ -115,7 +115,7 @@ module mld_c_ilu_solver & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_ilu_solver_type), intent(in) :: sv + class(mld_c_ilu_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_jac_smoother.f90 b/mlprec/mld_c_jac_smoother.f90 index 35ca1a5d..c6ea1489 100644 --- a/mlprec/mld_c_jac_smoother.f90 +++ b/mlprec/mld_c_jac_smoother.f90 @@ -93,7 +93,7 @@ module mld_c_jac_smoother import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_c_jac_smoother_type), intent(in) :: sm + class(mld_c_jac_smoother_type), intent(inout) :: sm complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index 0efb49ce..5f4bd1c9 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -157,7 +157,7 @@ module mld_c_prec_type subroutine mld_cprecaply(prec,x,y,desc_data,info,trans,work) import :: psb_cspmat_type, psb_desc_type, psb_spk_, mld_cprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec + type(mld_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -167,7 +167,7 @@ module mld_c_prec_type subroutine mld_cprecaply1(prec,x,desc_data,info,trans) import :: psb_cspmat_type, psb_desc_type, psb_spk_, mld_cprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_cprec_type), intent(in) :: prec + type(mld_cprec_type), intent(inout) :: prec complex(psb_spk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/mld_c_slu_solver.F90 b/mlprec/mld_c_slu_solver.F90 index 8e85421f..76774159 100644 --- a/mlprec/mld_c_slu_solver.F90 +++ b/mlprec/mld_c_slu_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_slu_solver_type), intent(in) :: sv + class(mld_c_slu_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_c_sludist_solver.F90 b/mlprec/mld_c_sludist_solver.F90 index 37f45be6..fac53f56 100644 --- a/mlprec/mld_c_sludist_solver.F90 +++ b/mlprec/mld_c_sludist_solver.F90 @@ -79,13 +79,12 @@ module mld_c_sludist_solver interface - function mld_csludist_fact(n,nnz,values,rowptr,colind,& - & lufactors)& + function mld_csludist_fact(n,nl,nnz,ifrst, & + & values,rowptr,colind,lufactors,npr,npc) & & bind(c,name='mld_csludist_fact') result(info) use iso_c_binding - integer(c_int), value :: n,nnz + integer(c_int), value :: n,nl,nnz,ifrst,npr,npc integer(c_int) :: info - !integer(c_long_long) :: ssize, nsize integer(c_int) :: rowptr(*),colind(*) complex(c_float_complex) :: values(*) type(c_ptr) :: lufactors @@ -93,12 +92,12 @@ module mld_c_sludist_solver end interface interface - function mld_csludist_solve(itrans,n,x, b, ldb, lufactors)& + function mld_csludist_solve(itrans,n,nrhs, b, ldb, lufactors)& & bind(c,name='mld_csludist_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_csludist_solve end interface @@ -118,7 +117,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_sludist_solver_type), intent(in) :: sv + class(mld_c_sludist_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta @@ -162,21 +161,24 @@ contains select case(trans_) case('N') - info = mld_csludist_solve(0,n_row,ww,x,n_row,sv%lufactors) + info = mld_csludist_solve(0,n_row,1,ww,n_row,sv%lufactors) case('T') - info = mld_csludist_solve(1,n_row,ww,x,n_row,sv%lufactors) + info = mld_csludist_solve(1,n_row,1,ww,n_row,sv%lufactors) case('C') - info = mld_csludist_solve(2,n_row,ww,x,n_row,sv%lufactors) + info = mld_csludist_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 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 @@ -253,7 +255,8 @@ contains ! Local variables type(psb_cspmat_type) :: atmp type(psb_c_csr_sparse_mat) :: acsr - integer :: n_row,n_col, nrow_a, nztota + integer :: n_row,n_col, nrow_a, nztota, nglob, nzt, npr, npc + integer :: ifrst, ibcheck integer :: ictxt,np,me,i, err_act, debug_unit, debug_level character(len=20) :: name='c_sludist_solver_bld', ch_err @@ -263,19 +266,18 @@ contains debug_level = psb_get_debug_level() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) + npr = np + npc = 1 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - write(0,*) 'SLUDIST INTERFACE IS CURRENTLY BROKEN. TO BE FIXED' - info=psb_err_internal_error_ - call psb_errpush(info,name) - goto 9999 - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() if (psb_toupper(upd) == 'F') then + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + nglob = desc_a%get_global_rows() + call a%cscnv(atmp,info,type='coo') call psb_rwextd(n_row,atmp,info,b=b) call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_) @@ -283,10 +285,15 @@ contains nrow_a = acsr%get_nrows() nztota = acsr%get_nzeros() ! Fix the entries to call C-base SuperLU + call psb_loc_to_glob(1,ifrst,desc_a,info) + call psb_loc_to_glob(nrow_a,ibcheck,desc_a,info) + call psb_loc_to_glob(acsr%ja(1:nztota),desc_a,info,iact='I') acsr%ja(:) = acsr%ja(:) - 1 acsr%irp(:) = acsr%irp(:) - 1 - info = mld_csludist_fact(nrow_a,nztota,acsr%val,& - & acsr%irp,acsr%ja,sv%lufactors) + ifrst = ifrst - 1 + info = mld_csludist_fact(nglob,nrow_a,nztota,ifrst,& + & acsr%val,acsr%irp,acsr%ja,sv%lufactors,& + & npr,npc) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -393,7 +400,7 @@ contains iout_ = 6 endif - write(iout_,*) ' SuperLU Sparse Factorization Solver. ' + write(iout_,*) ' SuperLU_Dist Sparse Factorization Solver. ' call psb_erractionrestore(err_act) return diff --git a/mlprec/mld_c_umf_solver.F90 b/mlprec/mld_c_umf_solver.F90 index 80ae1481..15a3a1b8 100644 --- a/mlprec/mld_c_umf_solver.F90 +++ b/mlprec/mld_c_umf_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_c_umf_solver_type), intent(in) :: sv + class(mld_c_umf_solver_type), intent(inout) :: sv complex(psb_spk_),intent(inout) :: x(:) complex(psb_spk_),intent(inout) :: y(:) complex(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_as_smoother.f90 b/mlprec/mld_d_as_smoother.f90 index b54b7cad..02e1d4f2 100644 --- a/mlprec/mld_d_as_smoother.f90 +++ b/mlprec/mld_d_as_smoother.f90 @@ -117,7 +117,7 @@ module mld_d_as_smoother & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, psb_desc_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_as_smoother_type), intent(in) :: sm + class(mld_d_as_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_base_smoother_mod.f90 b/mlprec/mld_d_base_smoother_mod.f90 index ae4668f7..396c25cf 100644 --- a/mlprec/mld_d_base_smoother_mod.f90 +++ b/mlprec/mld_d_base_smoother_mod.f90 @@ -127,7 +127,7 @@ module mld_d_base_smoother_mod & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & & mld_d_base_smoother_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_d_base_smoother_type), intent(in) :: sm + class(mld_d_base_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_base_solver_mod.f90 b/mlprec/mld_d_base_solver_mod.f90 index 49d01d26..3731d143 100644 --- a/mlprec/mld_d_base_solver_mod.f90 +++ b/mlprec/mld_d_base_solver_mod.f90 @@ -118,7 +118,7 @@ module mld_d_base_solver_mod & mld_d_base_solver_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_base_solver_type), intent(in) :: sv + class(mld_d_base_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_diag_solver.f90 b/mlprec/mld_d_diag_solver.f90 index cf8e9686..089ea2c4 100644 --- a/mlprec/mld_d_diag_solver.f90 +++ b/mlprec/mld_d_diag_solver.f90 @@ -90,7 +90,7 @@ module mld_d_diag_solver & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & & mld_d_diag_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_d_diag_solver_type), intent(in) :: sv + class(mld_d_diag_solver_type), intent(inout) :: sv real(psb_dpk_), intent(inout) :: x(:) real(psb_dpk_), intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_id_solver.f90 b/mlprec/mld_d_id_solver.f90 index 8b80a8d0..e93ec10f 100644 --- a/mlprec/mld_d_id_solver.f90 +++ b/mlprec/mld_d_id_solver.f90 @@ -85,7 +85,7 @@ module mld_d_id_solver & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & & mld_d_id_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_d_id_solver_type), intent(in) :: sv + class(mld_d_id_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_ilu_solver.f90 b/mlprec/mld_d_ilu_solver.f90 index 318c907e..4e45fed5 100644 --- a/mlprec/mld_d_ilu_solver.f90 +++ b/mlprec/mld_d_ilu_solver.f90 @@ -115,7 +115,7 @@ module mld_d_ilu_solver & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_ilu_solver_type), intent(in) :: sv + class(mld_d_ilu_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_jac_smoother.f90 b/mlprec/mld_d_jac_smoother.f90 index dd15da23..fb487df9 100644 --- a/mlprec/mld_d_jac_smoother.f90 +++ b/mlprec/mld_d_jac_smoother.f90 @@ -93,7 +93,7 @@ module mld_d_jac_smoother import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_d_jac_smoother_type), intent(in) :: sm + class(mld_d_jac_smoother_type), intent(inout) :: sm real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index c960cc5b..0c7047b7 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -157,7 +157,7 @@ module mld_d_prec_type subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) import :: psb_dspmat_type, psb_desc_type, psb_dpk_, mld_dprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec + type(mld_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -167,7 +167,7 @@ module mld_d_prec_type subroutine mld_dprecaply1(prec,x,desc_data,info,trans) import :: psb_dspmat_type, psb_desc_type, psb_dpk_, mld_dprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_dprec_type), intent(in) :: prec + type(mld_dprec_type), intent(inout) :: prec real(psb_dpk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/mld_d_slu_solver.F90 b/mlprec/mld_d_slu_solver.F90 index 1cf536ad..b5e036f6 100644 --- a/mlprec/mld_d_slu_solver.F90 +++ b/mlprec/mld_d_slu_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_slu_solver_type), intent(in) :: sv + class(mld_d_slu_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta @@ -301,9 +301,9 @@ contains call atmp%free() else ! ? - info=psb_err_internal_error_ - call psb_errpush(info,name) - goto 9999 + info=psb_err_internal_error_ + call psb_errpush(info,name) + goto 9999 end if diff --git a/mlprec/mld_d_sludist_solver.F90 b/mlprec/mld_d_sludist_solver.F90 index db97e483..fbe501d1 100644 --- a/mlprec/mld_d_sludist_solver.F90 +++ b/mlprec/mld_d_sludist_solver.F90 @@ -86,18 +86,18 @@ module mld_d_sludist_solver integer(c_int), value :: n,nl,nnz,ifrst,npr,npc integer(c_int) :: info integer(c_int) :: rowptr(*),colind(*) - real(c_double) :: values(*) + real(c_double) :: values(*) type(c_ptr) :: lufactors end function mld_dsludist_fact end interface interface - function mld_dsludist_solve(itrans,n,nrhs,b,ldb,lufactors)& + function mld_dsludist_solve(itrans,n,nrhs, b, ldb, lufactors)& & bind(c,name='mld_dsludist_solve') result(info) use iso_c_binding integer(c_int) :: info integer(c_int), value :: itrans,n,nrhs,ldb - real(c_double) :: b(ldb,*) + real(c_double) :: b(ldb,*) type(c_ptr), value :: lufactors end function mld_dsludist_solve end interface @@ -117,7 +117,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_sludist_solver_type), intent(in) :: sv + class(mld_d_sludist_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta @@ -158,7 +158,7 @@ contains goto 9999 end if endif - ww(1:n_row)=x(1:n_row) + select case(trans_) case('N') info = mld_dsludist_solve(0,n_row,1,ww,n_row,sv%lufactors) @@ -218,10 +218,6 @@ contains info = psb_success_ -!!$ write(0,*) 'SLUDIST INTERFACE IS CURRENTLY BROKEN. TO BE FIXED' -!!$ info=psb_err_internal_error_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 call x%v%sync() call y%v%sync() call sv%apply(alpha,x%v%v,beta,y%v%v,desc_data,trans,work,info) @@ -275,10 +271,6 @@ contains if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' -!!$ write(0,*) 'SLUDIST INTERFACE IS CURRENTLY BROKEN. TO BE FIXED' -!!$ info=psb_err_internal_error_ -!!$ call psb_errpush(info,name) -!!$ goto 9999 if (psb_toupper(upd) == 'F') then @@ -298,7 +290,6 @@ contains call psb_loc_to_glob(acsr%ja(1:nztota),desc_a,info,iact='I') acsr%ja(:) = acsr%ja(:) - 1 acsr%irp(:) = acsr%irp(:) - 1 - write(0,*) 'ACSR ',maxval(acsr%ja),minval(acsr%ja),nrow_a,nztota ifrst = ifrst - 1 info = mld_dsludist_fact(nglob,nrow_a,nztota,ifrst,& & acsr%val,acsr%irp,acsr%ja,sv%lufactors,& diff --git a/mlprec/mld_d_umf_solver.F90 b/mlprec/mld_d_umf_solver.F90 index 9b1a1e18..16c12044 100644 --- a/mlprec/mld_d_umf_solver.F90 +++ b/mlprec/mld_d_umf_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_d_umf_solver_type), intent(in) :: sv + class(mld_d_umf_solver_type), intent(inout) :: sv real(psb_dpk_),intent(inout) :: x(:) real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_as_smoother.f90 b/mlprec/mld_s_as_smoother.f90 index fb0a64e3..64dc566a 100644 --- a/mlprec/mld_s_as_smoother.f90 +++ b/mlprec/mld_s_as_smoother.f90 @@ -117,7 +117,7 @@ module mld_s_as_smoother & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, psb_desc_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_as_smoother_type), intent(in) :: sm + class(mld_s_as_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_base_smoother_mod.f90 b/mlprec/mld_s_base_smoother_mod.f90 index 9ad581d5..270f7dea 100644 --- a/mlprec/mld_s_base_smoother_mod.f90 +++ b/mlprec/mld_s_base_smoother_mod.f90 @@ -127,7 +127,7 @@ module mld_s_base_smoother_mod & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & & mld_s_base_smoother_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_s_base_smoother_type), intent(in) :: sm + class(mld_s_base_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_base_solver_mod.f90 b/mlprec/mld_s_base_solver_mod.f90 index f8dc2f50..a2226e53 100644 --- a/mlprec/mld_s_base_solver_mod.f90 +++ b/mlprec/mld_s_base_solver_mod.f90 @@ -118,7 +118,7 @@ module mld_s_base_solver_mod & mld_s_base_solver_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_base_solver_type), intent(in) :: sv + class(mld_s_base_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_diag_solver.f90 b/mlprec/mld_s_diag_solver.f90 index b0f94be8..8f7d239b 100644 --- a/mlprec/mld_s_diag_solver.f90 +++ b/mlprec/mld_s_diag_solver.f90 @@ -90,7 +90,7 @@ module mld_s_diag_solver & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & & mld_s_diag_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_s_diag_solver_type), intent(in) :: sv + class(mld_s_diag_solver_type), intent(inout) :: sv real(psb_spk_), intent(inout) :: x(:) real(psb_spk_), intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_id_solver.f90 b/mlprec/mld_s_id_solver.f90 index 59fe657d..72bed12c 100644 --- a/mlprec/mld_s_id_solver.f90 +++ b/mlprec/mld_s_id_solver.f90 @@ -85,7 +85,7 @@ module mld_s_id_solver & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & & mld_s_id_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_s_id_solver_type), intent(in) :: sv + class(mld_s_id_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_ilu_solver.f90 b/mlprec/mld_s_ilu_solver.f90 index 4d1a09a8..68297ee7 100644 --- a/mlprec/mld_s_ilu_solver.f90 +++ b/mlprec/mld_s_ilu_solver.f90 @@ -115,7 +115,7 @@ module mld_s_ilu_solver & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_ilu_solver_type), intent(in) :: sv + class(mld_s_ilu_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_jac_smoother.f90 b/mlprec/mld_s_jac_smoother.f90 index 58400bbf..43b73a5e 100644 --- a/mlprec/mld_s_jac_smoother.f90 +++ b/mlprec/mld_s_jac_smoother.f90 @@ -93,7 +93,7 @@ module mld_s_jac_smoother import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_s_jac_smoother_type), intent(in) :: sm + class(mld_s_jac_smoother_type), intent(inout) :: sm real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index 7d4818d4..28df6ab1 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -157,7 +157,7 @@ module mld_s_prec_type subroutine mld_sprecaply(prec,x,y,desc_data,info,trans,work) import :: psb_sspmat_type, psb_desc_type, psb_spk_, mld_sprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec + type(mld_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -167,7 +167,7 @@ module mld_s_prec_type subroutine mld_sprecaply1(prec,x,desc_data,info,trans) import :: psb_sspmat_type, psb_desc_type, psb_spk_, mld_sprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_sprec_type), intent(in) :: prec + type(mld_sprec_type), intent(inout) :: prec real(psb_spk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/mld_s_slu_solver.F90 b/mlprec/mld_s_slu_solver.F90 index d6202c22..e32ca169 100644 --- a/mlprec/mld_s_slu_solver.F90 +++ b/mlprec/mld_s_slu_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_slu_solver_type), intent(in) :: sv + class(mld_s_slu_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_s_sludist_solver.F90 b/mlprec/mld_s_sludist_solver.F90 index 83fb40d2..5808c716 100644 --- a/mlprec/mld_s_sludist_solver.F90 +++ b/mlprec/mld_s_sludist_solver.F90 @@ -79,13 +79,12 @@ module mld_s_sludist_solver interface - function mld_ssludist_fact(n,nnz,values,rowptr,colind,& - & lufactors)& + function mld_ssludist_fact(n,nl,nnz,ifrst, & + & values,rowptr,colind,lufactors,npr,npc) & & bind(c,name='mld_ssludist_fact') result(info) use iso_c_binding - integer(c_int), value :: n,nnz + integer(c_int), value :: n,nl,nnz,ifrst,npr,npc integer(c_int) :: info - !integer(c_long_long) :: ssize, nsize integer(c_int) :: rowptr(*),colind(*) real(c_float) :: values(*) type(c_ptr) :: lufactors @@ -93,12 +92,12 @@ module mld_s_sludist_solver end interface interface - function mld_ssludist_solve(itrans,n,x, b, ldb, lufactors)& + function mld_ssludist_solve(itrans,n,nrhs, b, ldb, lufactors)& & bind(c,name='mld_ssludist_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_ssludist_solve end interface @@ -118,7 +117,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_sludist_solver_type), intent(in) :: sv + class(mld_s_sludist_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta @@ -162,21 +161,24 @@ contains select case(trans_) case('N') - info = mld_ssludist_solve(0,n_row,ww,x,n_row,sv%lufactors) + info = mld_ssludist_solve(0,n_row,1,ww,n_row,sv%lufactors) case('T') - info = mld_ssludist_solve(1,n_row,ww,x,n_row,sv%lufactors) + info = mld_ssludist_solve(1,n_row,1,ww,n_row,sv%lufactors) case('C') - info = mld_ssludist_solve(2,n_row,ww,x,n_row,sv%lufactors) + info = mld_ssludist_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 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 @@ -253,7 +255,8 @@ contains ! Local variables type(psb_sspmat_type) :: atmp type(psb_s_csr_sparse_mat) :: acsr - integer :: n_row,n_col, nrow_a, nztota + integer :: n_row,n_col, nrow_a, nztota, nglob, nzt, npr, npc + integer :: ifrst, ibcheck integer :: ictxt,np,me,i, err_act, debug_unit, debug_level character(len=20) :: name='s_sludist_solver_bld', ch_err @@ -263,19 +266,18 @@ contains debug_level = psb_get_debug_level() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) + npr = np + npc = 1 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - write(0,*) 'SLUDIST INTERFACE IS CURRENTLY BROKEN. TO BE FIXED' - info=psb_err_internal_error_ - call psb_errpush(info,name) - goto 9999 - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() if (psb_toupper(upd) == 'F') then + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + nglob = desc_a%get_global_rows() + call a%cscnv(atmp,info,type='coo') call psb_rwextd(n_row,atmp,info,b=b) call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_) @@ -283,10 +285,15 @@ contains nrow_a = acsr%get_nrows() nztota = acsr%get_nzeros() ! Fix the entries to call C-base SuperLU + call psb_loc_to_glob(1,ifrst,desc_a,info) + call psb_loc_to_glob(nrow_a,ibcheck,desc_a,info) + call psb_loc_to_glob(acsr%ja(1:nztota),desc_a,info,iact='I') acsr%ja(:) = acsr%ja(:) - 1 acsr%irp(:) = acsr%irp(:) - 1 - info = mld_ssludist_fact(nrow_a,nztota,acsr%val,& - & acsr%irp,acsr%ja,sv%lufactors) + ifrst = ifrst - 1 + info = mld_ssludist_fact(nglob,nrow_a,nztota,ifrst,& + & acsr%val,acsr%irp,acsr%ja,sv%lufactors,& + & npr,npc) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -393,7 +400,7 @@ contains iout_ = 6 endif - write(iout_,*) ' SuperLU Sparse Factorization Solver. ' + write(iout_,*) ' SuperLU_Dist Sparse Factorization Solver. ' call psb_erractionrestore(err_act) return diff --git a/mlprec/mld_s_umf_solver.F90 b/mlprec/mld_s_umf_solver.F90 index 6298523b..1cf3a41c 100644 --- a/mlprec/mld_s_umf_solver.F90 +++ b/mlprec/mld_s_umf_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_s_umf_solver_type), intent(in) :: sv + class(mld_s_umf_solver_type), intent(inout) :: sv real(psb_spk_),intent(inout) :: x(:) real(psb_spk_),intent(inout) :: y(:) real(psb_spk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_as_smoother.f90 b/mlprec/mld_z_as_smoother.f90 index 6fdd18f7..3cdceb85 100644 --- a/mlprec/mld_z_as_smoother.f90 +++ b/mlprec/mld_z_as_smoother.f90 @@ -117,7 +117,7 @@ module mld_z_as_smoother & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, psb_desc_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_as_smoother_type), intent(in) :: sm + class(mld_z_as_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_base_smoother_mod.f90 b/mlprec/mld_z_base_smoother_mod.f90 index e7b91656..d03a4836 100644 --- a/mlprec/mld_z_base_smoother_mod.f90 +++ b/mlprec/mld_z_base_smoother_mod.f90 @@ -127,7 +127,7 @@ module mld_z_base_smoother_mod & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & & mld_z_base_smoother_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_z_base_smoother_type), intent(in) :: sm + class(mld_z_base_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_base_solver_mod.f90 b/mlprec/mld_z_base_solver_mod.f90 index 55b9cb73..e4841d0a 100644 --- a/mlprec/mld_z_base_solver_mod.f90 +++ b/mlprec/mld_z_base_solver_mod.f90 @@ -118,7 +118,7 @@ module mld_z_base_solver_mod & mld_z_base_solver_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_base_solver_type), intent(in) :: sv + class(mld_z_base_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_diag_solver.f90 b/mlprec/mld_z_diag_solver.f90 index 52b6af96..982b2545 100644 --- a/mlprec/mld_z_diag_solver.f90 +++ b/mlprec/mld_z_diag_solver.f90 @@ -90,7 +90,7 @@ module mld_z_diag_solver & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & & mld_z_diag_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_z_diag_solver_type), intent(in) :: sv + class(mld_z_diag_solver_type), intent(inout) :: sv complex(psb_dpk_), intent(inout) :: x(:) complex(psb_dpk_), intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_id_solver.f90 b/mlprec/mld_z_id_solver.f90 index 6b6b6c2c..fc0d45c5 100644 --- a/mlprec/mld_z_id_solver.f90 +++ b/mlprec/mld_z_id_solver.f90 @@ -85,7 +85,7 @@ module mld_z_id_solver & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & & mld_z_id_solver_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_z_id_solver_type), intent(in) :: sv + class(mld_z_id_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_ilu_solver.f90 b/mlprec/mld_z_ilu_solver.f90 index 83e3ea17..9b87f37d 100644 --- a/mlprec/mld_z_ilu_solver.f90 +++ b/mlprec/mld_z_ilu_solver.f90 @@ -115,7 +115,7 @@ module mld_z_ilu_solver & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_ implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_ilu_solver_type), intent(in) :: sv + class(mld_z_ilu_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_jac_smoother.f90 b/mlprec/mld_z_jac_smoother.f90 index b61a8ec2..61246a3d 100644 --- a/mlprec/mld_z_jac_smoother.f90 +++ b/mlprec/mld_z_jac_smoother.f90 @@ -93,7 +93,7 @@ module mld_z_jac_smoother import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_ type(psb_desc_type), intent(in) :: desc_data - class(mld_z_jac_smoother_type), intent(in) :: sm + class(mld_z_jac_smoother_type), intent(inout) :: sm complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index 1763f447..f50f8d3f 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -157,7 +157,7 @@ module mld_z_prec_type subroutine mld_zprecaply(prec,x,y,desc_data,info,trans,work) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, mld_zprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec + type(mld_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) integer(psb_ipk_), intent(out) :: info @@ -167,7 +167,7 @@ module mld_z_prec_type subroutine mld_zprecaply1(prec,x,desc_data,info,trans) import :: psb_zspmat_type, psb_desc_type, psb_dpk_, mld_zprec_type, psb_ipk_ type(psb_desc_type),intent(in) :: desc_data - type(mld_zprec_type), intent(in) :: prec + type(mld_zprec_type), intent(inout) :: prec complex(psb_dpk_),intent(inout) :: x(:) integer(psb_ipk_), intent(out) :: info character(len=1), optional :: trans diff --git a/mlprec/mld_z_slu_solver.F90 b/mlprec/mld_z_slu_solver.F90 index d8c4d99f..c1e0e111 100644 --- a/mlprec/mld_z_slu_solver.F90 +++ b/mlprec/mld_z_slu_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_slu_solver_type), intent(in) :: sv + class(mld_z_slu_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta diff --git a/mlprec/mld_z_sludist_solver.F90 b/mlprec/mld_z_sludist_solver.F90 index d15f272b..e5c2c0c5 100644 --- a/mlprec/mld_z_sludist_solver.F90 +++ b/mlprec/mld_z_sludist_solver.F90 @@ -79,13 +79,12 @@ module mld_z_sludist_solver interface - function mld_zsludist_fact(n,nnz,values,rowptr,colind,& - & lufactors)& + function mld_zsludist_fact(n,nl,nnz,ifrst, & + & values,rowptr,colind,lufactors,npr,npc) & & bind(c,name='mld_zsludist_fact') result(info) use iso_c_binding - integer(c_int), value :: n,nnz + integer(c_int), value :: n,nl,nnz,ifrst,npr,npc integer(c_int) :: info - !integer(c_long_long) :: ssize, nsize integer(c_int) :: rowptr(*),colind(*) complex(c_double_complex) :: values(*) type(c_ptr) :: lufactors @@ -93,12 +92,12 @@ module mld_z_sludist_solver end interface interface - function mld_zsludist_solve(itrans,n,x, b, ldb, lufactors)& + function mld_zsludist_solve(itrans,n,nrhs, b, ldb, lufactors)& & bind(c,name='mld_zsludist_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_zsludist_solve end interface @@ -118,7 +117,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_sludist_solver_type), intent(in) :: sv + class(mld_z_sludist_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta @@ -162,21 +161,24 @@ contains select case(trans_) case('N') - info = mld_zsludist_solve(0,n_row,ww,x,n_row,sv%lufactors) + info = mld_zsludist_solve(0,n_row,1,ww,n_row,sv%lufactors) case('T') - info = mld_zsludist_solve(1,n_row,ww,x,n_row,sv%lufactors) + info = mld_zsludist_solve(1,n_row,1,ww,n_row,sv%lufactors) case('C') - info = mld_zsludist_solve(2,n_row,ww,x,n_row,sv%lufactors) + info = mld_zsludist_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 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 @@ -253,7 +255,8 @@ contains ! Local variables type(psb_zspmat_type) :: atmp type(psb_z_csr_sparse_mat) :: acsr - integer :: n_row,n_col, nrow_a, nztota + integer :: n_row,n_col, nrow_a, nztota, nglob, nzt, npr, npc + integer :: ifrst, ibcheck integer :: ictxt,np,me,i, err_act, debug_unit, debug_level character(len=20) :: name='z_sludist_solver_bld', ch_err @@ -263,19 +266,18 @@ contains debug_level = psb_get_debug_level() ictxt = desc_a%get_context() call psb_info(ictxt, me, np) + npr = np + npc = 1 if (debug_level >= psb_debug_outer_) & & write(debug_unit,*) me,' ',trim(name),' start' - write(0,*) 'SLUDIST INTERFACE IS CURRENTLY BROKEN. TO BE FIXED' - info=psb_err_internal_error_ - call psb_errpush(info,name) - goto 9999 - - n_row = desc_a%get_local_rows() - n_col = desc_a%get_local_cols() if (psb_toupper(upd) == 'F') then + n_row = desc_a%get_local_rows() + n_col = desc_a%get_local_cols() + nglob = desc_a%get_global_rows() + call a%cscnv(atmp,info,type='coo') call psb_rwextd(n_row,atmp,info,b=b) call atmp%cscnv(info,type='csr',dupl=psb_dupl_add_) @@ -283,10 +285,15 @@ contains nrow_a = acsr%get_nrows() nztota = acsr%get_nzeros() ! Fix the entries to call C-base SuperLU + call psb_loc_to_glob(1,ifrst,desc_a,info) + call psb_loc_to_glob(nrow_a,ibcheck,desc_a,info) + call psb_loc_to_glob(acsr%ja(1:nztota),desc_a,info,iact='I') acsr%ja(:) = acsr%ja(:) - 1 acsr%irp(:) = acsr%irp(:) - 1 - info = mld_zsludist_fact(nrow_a,nztota,acsr%val,& - & acsr%irp,acsr%ja,sv%lufactors) + ifrst = ifrst - 1 + info = mld_zsludist_fact(nglob,nrow_a,nztota,ifrst,& + & acsr%val,acsr%irp,acsr%ja,sv%lufactors,& + & npr,npc) if (info /= psb_success_) then info=psb_err_from_subroutine_ @@ -393,7 +400,7 @@ contains iout_ = 6 endif - write(iout_,*) ' SuperLU Sparse Factorization Solver. ' + write(iout_,*) ' SuperLU_Dist Sparse Factorization Solver. ' call psb_erractionrestore(err_act) return diff --git a/mlprec/mld_z_umf_solver.F90 b/mlprec/mld_z_umf_solver.F90 index e017d882..d9352cf4 100644 --- a/mlprec/mld_z_umf_solver.F90 +++ b/mlprec/mld_z_umf_solver.F90 @@ -120,7 +120,7 @@ contains use psb_base_mod implicit none type(psb_desc_type), intent(in) :: desc_data - class(mld_z_umf_solver_type), intent(in) :: sv + class(mld_z_umf_solver_type), intent(inout) :: sv complex(psb_dpk_),intent(inout) :: x(:) complex(psb_dpk_),intent(inout) :: y(:) complex(psb_dpk_),intent(in) :: alpha,beta