From da1c129ce592853b2c2c1ca0880ccb20580f7a51 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Sun, 28 Feb 2016 15:18:29 +0000 Subject: [PATCH] mld2p4-2: mlprec/Makefile mlprec/impl/mld_ccprecset.F90 mlprec/impl/mld_cprecset.F90 mlprec/impl/mld_dcprecset.F90 mlprec/impl/mld_dprecset.F90 mlprec/impl/mld_scprecset.F90 mlprec/impl/mld_sprecset.F90 mlprec/impl/mld_zcprecset.F90 mlprec/impl/mld_zprecset.F90 mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 mlprec/impl/solver/Makefile mlprec/impl/solver/mld_c_gs_solver_apply.f90 mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 mlprec/impl/solver/mld_c_gs_solver_bld.f90 mlprec/impl/solver/mld_c_gs_solver_clone.f90 mlprec/impl/solver/mld_c_gs_solver_cnv.f90 mlprec/impl/solver/mld_c_gs_solver_dmp.f90 mlprec/impl/solver/mld_d_gs_solver_bld.f90 mlprec/impl/solver/mld_s_gs_solver_apply.f90 mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 mlprec/impl/solver/mld_s_gs_solver_bld.f90 mlprec/impl/solver/mld_s_gs_solver_clone.f90 mlprec/impl/solver/mld_s_gs_solver_cnv.f90 mlprec/impl/solver/mld_s_gs_solver_dmp.f90 mlprec/impl/solver/mld_z_gs_solver_apply.f90 mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 mlprec/impl/solver/mld_z_gs_solver_bld.f90 mlprec/impl/solver/mld_z_gs_solver_clone.f90 mlprec/impl/solver/mld_z_gs_solver_cnv.f90 mlprec/impl/solver/mld_z_gs_solver_dmp.f90 mlprec/mld_c_base_solver_mod.f90 mlprec/mld_c_gs_solver.f90 mlprec/mld_c_prec_mod.f90 mlprec/mld_c_prec_type.f90 mlprec/mld_d_base_solver_mod.f90 mlprec/mld_s_base_solver_mod.f90 mlprec/mld_s_gs_solver.f90 mlprec/mld_s_prec_mod.f90 mlprec/mld_s_prec_type.f90 mlprec/mld_z_base_solver_mod.f90 mlprec/mld_z_gs_solver.f90 mlprec/mld_z_prec_mod.f90 mlprec/mld_z_prec_type.f90 Templated GS solver. --- mlprec/Makefile | 12 +- mlprec/impl/mld_ccprecset.F90 | 170 +++--- mlprec/impl/mld_cprecset.F90 | 170 +++--- mlprec/impl/mld_dcprecset.F90 | 42 +- mlprec/impl/mld_dprecset.F90 | 27 +- mlprec/impl/mld_scprecset.F90 | 170 +++--- mlprec/impl/mld_sprecset.F90 | 170 +++--- mlprec/impl/mld_zcprecset.F90 | 240 ++++---- mlprec/impl/mld_zprecset.F90 | 240 ++++---- .../smoother/mld_c_as_smoother_apply_vect.f90 | 2 +- .../mld_c_jac_smoother_apply_vect.f90 | 7 +- .../smoother/mld_s_as_smoother_apply_vect.f90 | 2 +- .../mld_s_jac_smoother_apply_vect.f90 | 7 +- .../smoother/mld_z_as_smoother_apply_vect.f90 | 2 +- .../mld_z_jac_smoother_apply_vect.f90 | 7 +- mlprec/impl/solver/Makefile | 18 + mlprec/impl/solver/mld_c_gs_solver_apply.f90 | 196 +++++++ .../solver/mld_c_gs_solver_apply_vect.f90 | 200 +++++++ mlprec/impl/solver/mld_c_gs_solver_bld.f90 | 117 ++++ mlprec/impl/solver/mld_c_gs_solver_clone.f90 | 88 +++ mlprec/impl/solver/mld_c_gs_solver_cnv.f90 | 75 +++ mlprec/impl/solver/mld_c_gs_solver_dmp.f90 | 87 +++ mlprec/impl/solver/mld_d_gs_solver_bld.f90 | 1 - mlprec/impl/solver/mld_s_gs_solver_apply.f90 | 196 +++++++ .../solver/mld_s_gs_solver_apply_vect.f90 | 200 +++++++ mlprec/impl/solver/mld_s_gs_solver_bld.f90 | 117 ++++ mlprec/impl/solver/mld_s_gs_solver_clone.f90 | 88 +++ mlprec/impl/solver/mld_s_gs_solver_cnv.f90 | 75 +++ mlprec/impl/solver/mld_s_gs_solver_dmp.f90 | 87 +++ mlprec/impl/solver/mld_z_gs_solver_apply.f90 | 196 +++++++ .../solver/mld_z_gs_solver_apply_vect.f90 | 200 +++++++ mlprec/impl/solver/mld_z_gs_solver_bld.f90 | 117 ++++ mlprec/impl/solver/mld_z_gs_solver_clone.f90 | 88 +++ mlprec/impl/solver/mld_z_gs_solver_cnv.f90 | 75 +++ mlprec/impl/solver/mld_z_gs_solver_dmp.f90 | 87 +++ mlprec/mld_c_base_solver_mod.f90 | 14 +- mlprec/mld_c_gs_solver.f90 | 518 ++++++++++++++++++ mlprec/mld_c_prec_mod.f90 | 1 + mlprec/mld_c_prec_type.f90 | 2 +- mlprec/mld_d_base_solver_mod.f90 | 4 +- mlprec/mld_s_base_solver_mod.f90 | 14 +- mlprec/mld_s_gs_solver.f90 | 518 ++++++++++++++++++ mlprec/mld_s_prec_mod.f90 | 1 + mlprec/mld_s_prec_type.f90 | 2 +- mlprec/mld_z_base_solver_mod.f90 | 14 +- mlprec/mld_z_gs_solver.f90 | 518 ++++++++++++++++++ mlprec/mld_z_prec_mod.f90 | 1 + mlprec/mld_z_prec_type.f90 | 2 +- 48 files changed, 4684 insertions(+), 501 deletions(-) create mode 100644 mlprec/impl/solver/mld_c_gs_solver_apply.f90 create mode 100644 mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 create mode 100644 mlprec/impl/solver/mld_c_gs_solver_bld.f90 create mode 100644 mlprec/impl/solver/mld_c_gs_solver_clone.f90 create mode 100644 mlprec/impl/solver/mld_c_gs_solver_cnv.f90 create mode 100644 mlprec/impl/solver/mld_c_gs_solver_dmp.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_apply.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_bld.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_clone.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_cnv.f90 create mode 100644 mlprec/impl/solver/mld_s_gs_solver_dmp.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_apply.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_bld.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_clone.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_cnv.f90 create mode 100644 mlprec/impl/solver/mld_z_gs_solver_dmp.f90 create mode 100644 mlprec/mld_c_gs_solver.f90 create mode 100644 mlprec/mld_s_gs_solver.f90 create mode 100644 mlprec/mld_z_gs_solver.f90 diff --git a/mlprec/Makefile b/mlprec/Makefile index c32d6d4f..5fd1e7e1 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -16,17 +16,17 @@ DMODOBJS=mld_d_prec_type.o mld_d_ilu_fact_mod.o \ SMODOBJS=mld_s_prec_type.o mld_s_ilu_fact_mod.o \ mld_s_inner_mod.o mld_s_ilu_solver.o mld_s_diag_solver.o mld_s_jac_smoother.o mld_s_as_smoother.o \ mld_s_slu_solver.o mld_s_sludist_solver.o mld_s_id_solver.o\ - mld_s_base_solver_mod.o mld_s_base_smoother_mod.o mld_s_onelev_mod.o + mld_s_base_solver_mod.o mld_s_base_smoother_mod.o mld_s_onelev_mod.o mld_s_gs_solver.o ZMODOBJS=mld_z_prec_type.o mld_z_ilu_fact_mod.o \ mld_z_inner_mod.o mld_z_ilu_solver.o mld_z_diag_solver.o mld_z_jac_smoother.o mld_z_as_smoother.o \ mld_z_umf_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o mld_z_id_solver.o\ - mld_z_base_solver_mod.o mld_z_base_smoother_mod.o mld_z_onelev_mod.o + mld_z_base_solver_mod.o mld_z_base_smoother_mod.o mld_z_onelev_mod.o mld_z_gs_solver.o CMODOBJS=mld_c_prec_type.o mld_c_ilu_fact_mod.o \ mld_c_inner_mod.o mld_c_ilu_solver.o mld_c_diag_solver.o mld_c_jac_smoother.o mld_c_as_smoother.o \ mld_c_slu_solver.o mld_c_sludist_solver.o mld_c_id_solver.o\ - mld_c_base_solver_mod.o mld_c_base_smoother_mod.o mld_c_onelev_mod.o + mld_c_base_solver_mod.o mld_c_base_smoother_mod.o mld_c_onelev_mod.o mld_c_gs_solver.o @@ -104,7 +104,7 @@ mld_dprecinit.o mld_dprecset.o: mld_d_diag_solver.o mld_d_ilu_solver.o \ mld_d_umf_solver.o mld_d_as_smoother.o mld_d_jac_smoother.o \ mld_d_id_solver.o mld_d_slu_solver.o mld_d_sludist_solver.o -mld_s_id_solver.o mld_s_sludist_solver.o mld_s_slu_solver.o \ +mld_s_gs_solver.o mld_s_id_solver.o mld_s_sludist_solver.o mld_s_slu_solver.o \ mld_s_diag_solver.o mld_s_ilu_solver.o: mld_s_base_solver_mod.o mld_s_prec_type.o mld_s_ilu_fact_mod.o: mld_base_prec_type.o mld_s_base_solver_mod.o mld_s_ilu_solver.o mld_s_iluk_fact.o: mld_s_ilu_fact_mod.o @@ -114,7 +114,7 @@ mld_sprecinit.o mld_sprecset.o: mld_s_diag_solver.o mld_s_ilu_solver.o \ mld_s_as_smoother.o mld_s_jac_smoother.o \ mld_s_id_solver.o mld_s_slu_solver.o mld_s_sludist_solver.o -mld_z_id_solver.o mld_z_sludist_solver.o mld_z_slu_solver.o \ +mld_z_gs_solver.o mld_z_id_solver.o mld_z_sludist_solver.o mld_z_slu_solver.o \ mld_z_umf_solver.o mld_z_diag_solver.o mld_z_ilu_solver.o: mld_z_base_solver_mod.o mld_z_prec_type.o mld_z_ilu_fact_mod.o: mld_base_prec_type.o mld_z_base_solver_mod.o mld_z_ilu_solver.o mld_z_iluk_fact.o: mld_z_ilu_fact_mod.o @@ -124,7 +124,7 @@ mld_zprecinit.o mld_zprecset.o: mld_z_diag_solver.o mld_z_ilu_solver.o \ mld_z_umf_solver.o mld_z_as_smoother.o mld_z_jac_smoother.o \ mld_z_id_solver.o mld_z_slu_solver.o mld_z_sludist_solver.o -mld_c_id_solver.o mld_c_sludist_solver.o mld_c_slu_solver.o \ +mld_c_gs_solver.o mld_c_id_solver.o mld_c_sludist_solver.o mld_c_slu_solver.o \ mld_c_diag_solver.o mld_c_ilu_solver.o: mld_c_base_solver_mod.o mld_c_prec_type.o mld_c_ilu_fact_mod.o: mld_base_prec_type.o mld_c_base_solver_mod.o mld_c_ilu_solver.o mld_c_iluk_fact.o: mld_c_ilu_fact_mod.o diff --git a/mlprec/impl/mld_ccprecset.F90 b/mlprec/impl/mld_ccprecset.F90 index 2ab44055..bffa0160 100644 --- a/mlprec/impl/mld_ccprecset.F90 +++ b/mlprec/impl/mld_ccprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) use mld_c_diag_solver use mld_c_ilu_solver use mld_c_id_solver + use mld_c_gs_solver #if defined(HAVE_SLU_) use mld_c_slu_solver #endif @@ -445,83 +446,126 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_diag_solver_type) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_diag_solver_type) ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_diag_solver_type :: level%sm%sv, stat=info) - endif + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + else + endif + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set('SUB_SOLVE',val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set('SUB_SOLVE',val,info) - + #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/mld_cprecset.F90 b/mlprec/impl/mld_cprecset.F90 index 8efd77f8..bd297d58 100644 --- a/mlprec/impl/mld_cprecset.F90 +++ b/mlprec/impl/mld_cprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_cprecseti(p,what,val,info,ilev) use mld_c_diag_solver use mld_c_ilu_solver use mld_c_id_solver + use mld_c_gs_solver #if defined(HAVE_SLU_) use mld_c_slu_solver #endif @@ -446,82 +447,125 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - + case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_diag_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_diag_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_diag_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + end if + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set(mld_sub_solve_,val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set(mld_sub_solve_,val,info) #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_c_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_c_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_c_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_c_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_c_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_c_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/mld_dcprecset.F90 b/mlprec/impl/mld_dcprecset.F90 index 648a7a6f..c3c3d66f 100644 --- a/mlprec/impl/mld_dcprecset.F90 +++ b/mlprec/impl/mld_dcprecset.F90 @@ -459,9 +459,9 @@ contains if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_id_solver_type) - ! do nothing - class default + class is (mld_d_id_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_id_solver_type ::& @@ -478,14 +478,14 @@ contains write(0,*) 'Calling set_solver without a smoother?' info = -5 end if - + case (mld_diag_scale_) if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_diag_solver_type) - ! do nothing - class default + class is (mld_d_diag_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_diag_solver_type ::& @@ -528,15 +528,13 @@ contains info = -5 end if - - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_ilu_solver_type) - ! do nothing - class default + class is (mld_d_ilu_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_ilu_solver_type ::& @@ -554,15 +552,15 @@ contains write(0,*) 'Calling set_solver without a smoother?' info = -5 end if - + #ifdef HAVE_SLU_ case (mld_slu_) if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_slu_solver_type) - ! do nothing - class default + class is (mld_d_slu_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_slu_solver_type ::& @@ -585,9 +583,9 @@ contains if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_umf_solver_type) - ! do nothing - class default + class is (mld_d_umf_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_umf_solver_type ::& @@ -610,9 +608,9 @@ contains if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_sludist_solver_type) - ! do nothing - class default + class is (mld_d_sludist_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_sludist_solver_type ::& diff --git a/mlprec/impl/mld_dprecset.F90 b/mlprec/impl/mld_dprecset.F90 index 41cc1700..cc26afd0 100644 --- a/mlprec/impl/mld_dprecset.F90 +++ b/mlprec/impl/mld_dprecset.F90 @@ -529,7 +529,6 @@ contains info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) if (allocated(level%sm)) then if (allocated(level%sm%sv)) then @@ -559,9 +558,9 @@ contains if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_umf_solver_type) - ! do nothing - class default + class is (mld_d_umf_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_umf_solver_type ::& @@ -584,9 +583,9 @@ contains if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_sludist_solver_type) - ! do nothing - class default + class is (mld_d_sludist_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_sludist_solver_type ::& @@ -606,13 +605,12 @@ contains #endif #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm)) then if (allocated(level%sm%sv)) then select type (sv => level%sm%sv) - class is (mld_d_slu_solver_type) - ! do nothing - class default + class is (mld_d_slu_solver_type) + ! do nothing + class default call level%sm%sv%free(info) if (info == 0) deallocate(level%sm%sv) if (info == 0) allocate(mld_d_slu_solver_type ::& @@ -621,12 +619,14 @@ contains else allocate(mld_d_slu_solver_type :: level%sm%sv, stat=info) endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if else write(0,*) 'Calling set_solver without a smoother?' info = -5 end if - - #endif case default ! @@ -634,7 +634,6 @@ contains ! end select - end subroutine onelev_set_solver diff --git a/mlprec/impl/mld_scprecset.F90 b/mlprec/impl/mld_scprecset.F90 index dfff440e..6e7cc39f 100644 --- a/mlprec/impl/mld_scprecset.F90 +++ b/mlprec/impl/mld_scprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_scprecseti(p,what,val,info,ilev) use mld_s_diag_solver use mld_s_ilu_solver use mld_s_id_solver + use mld_s_gs_solver #if defined(HAVE_SLU_) use mld_s_slu_solver #endif @@ -445,83 +446,126 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_diag_solver_type) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_diag_solver_type) ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_diag_solver_type :: level%sm%sv, stat=info) - endif + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + else + endif + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set('SUB_SOLVE',val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set('SUB_SOLVE',val,info) - + #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/mld_sprecset.F90 b/mlprec/impl/mld_sprecset.F90 index 50f22fb6..00df9c25 100644 --- a/mlprec/impl/mld_sprecset.F90 +++ b/mlprec/impl/mld_sprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_sprecseti(p,what,val,info,ilev) use mld_s_diag_solver use mld_s_ilu_solver use mld_s_id_solver + use mld_s_gs_solver #if defined(HAVE_SLU_) use mld_s_slu_solver #endif @@ -446,82 +447,125 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - + case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_diag_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_diag_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_diag_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + end if + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set(mld_sub_solve_,val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set(mld_sub_solve_,val,info) #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_s_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_s_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_s_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_s_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_s_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_s_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/mld_zcprecset.F90 b/mlprec/impl/mld_zcprecset.F90 index 1fc07169..b0e14b0c 100644 --- a/mlprec/impl/mld_zcprecset.F90 +++ b/mlprec/impl/mld_zcprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) use mld_z_diag_solver use mld_z_ilu_solver use mld_z_id_solver + use mld_z_gs_solver #if defined(HAVE_UMF_) use mld_z_umf_solver #endif @@ -455,123 +456,176 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_diag_solver_type) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_diag_solver_type) ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_diag_solver_type :: level%sm%sv, stat=info) - endif + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + else + endif + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set('SUB_SOLVE',val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set('SUB_SOLVE',val,info) - + #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif #ifdef HAVE_UMF_ case (mld_umf_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_umf_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_umf_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_umf_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_umf_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_umf_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_umf_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif #ifdef HAVE_SLUDIST_ case (mld_sludist_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_sludist_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_sludist_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_sludist_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_sludist_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_sludist_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_sludist_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/mld_zprecset.F90 b/mlprec/impl/mld_zprecset.F90 index c561b31b..4ba61c4f 100644 --- a/mlprec/impl/mld_zprecset.F90 +++ b/mlprec/impl/mld_zprecset.F90 @@ -85,6 +85,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev) use mld_z_diag_solver use mld_z_ilu_solver use mld_z_id_solver + use mld_z_gs_solver #if defined(HAVE_UMF_) use mld_z_umf_solver #endif @@ -456,122 +457,175 @@ contains ! select case (val) case (mld_f_none_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_id_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_id_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_id_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_id_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_id_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_id_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - + case (mld_diag_scale_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_diag_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_diag_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_diag_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_diag_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_diag_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_diag_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 + end if + + case (mld_gs_) + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_gs_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_gs_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_gs_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm%sv)) then + call level%sm%sv%default() + end if + + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - case (mld_ilu_n_,mld_milu_n_,mld_ilu_t_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_ilu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_ilu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_ilu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_ilu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_ilu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_ilu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + call level%sm%sv%set(mld_sub_solve_,val,info) + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if - call level%sm%sv%set(mld_sub_solve_,val,info) #ifdef HAVE_UMF_ case (mld_umf_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_umf_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_umf_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_umf_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_umf_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_umf_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_umf_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif #ifdef HAVE_SLUDIST_ case (mld_sludist_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_sludist_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_sludist_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_sludist_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_sludist_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_sludist_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_sludist_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif #ifdef HAVE_SLU_ case (mld_slu_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_slu_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_slu_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_slu_solver_type :: level%sm%sv, stat=info) - endif if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_slu_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_slu_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_slu_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if + else + write(0,*) 'Calling set_solver without a smoother?' + info = -5 end if #endif case default diff --git a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 index 3ae21122..88b06239 100644 --- a/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90 @@ -124,7 +124,7 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& endif - if ((sm%novr == 0).and.(sweeps == 1)) then + if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then ! ! Shortcut: in this case it's just the same ! as Block Jacobi. diff --git a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 index 8096238b..9bd519e6 100644 --- a/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90 @@ -107,7 +107,9 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then +!!$ write(0,*) 'Jacobi smoother with ',sweeps + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -117,8 +119,7 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else if (sweeps > 1) then - + else if (sweeps >= 1) then ! ! ! Apply multiple sweeps of a block-Jacobi solver diff --git a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 index 36d542d3..31883887 100644 --- a/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90 @@ -124,7 +124,7 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& endif - if ((sm%novr == 0).and.(sweeps == 1)) then + if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then ! ! Shortcut: in this case it's just the same ! as Block Jacobi. diff --git a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 index 74e057b2..adc4f331 100644 --- a/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90 @@ -107,7 +107,9 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then +!!$ write(0,*) 'Jacobi smoother with ',sweeps + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -117,8 +119,7 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else if (sweeps > 1) then - + else if (sweeps >= 1) then ! ! ! Apply multiple sweeps of a block-Jacobi solver diff --git a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 index c36ea5ac..6e72e59b 100644 --- a/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90 @@ -124,7 +124,7 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& endif - if ((sm%novr == 0).and.(sweeps == 1)) then + if ((sm%novr == 0).and.(sweeps == 1).and.(.not.sm%sv%is_iterative())) then ! ! Shortcut: in this case it's just the same ! as Block Jacobi. diff --git a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 index f509550e..28f75dc8 100644 --- a/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 +++ b/mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90 @@ -107,7 +107,9 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& end if endif - if ((sweeps == 1).or.(sm%nnz_nd_tot==0)) then +!!$ write(0,*) 'Jacobi smoother with ',sweeps + + if ((.not.sm%sv%is_iterative()).and.((sweeps == 1).or.(sm%nnz_nd_tot==0))) then call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -117,8 +119,7 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& goto 9999 endif - else if (sweeps > 1) then - + else if (sweeps >= 1) then ! ! ! Apply multiple sweeps of a block-Jacobi solver diff --git a/mlprec/impl/solver/Makefile b/mlprec/impl/solver/Makefile index 7a728329..9c863c33 100644 --- a/mlprec/impl/solver/Makefile +++ b/mlprec/impl/solver/Makefile @@ -28,6 +28,12 @@ mld_c_diag_solver_apply_vect.o \ mld_c_diag_solver_bld.o \ mld_c_diag_solver_clone.o \ mld_c_diag_solver_cnv.o \ +mld_c_gs_solver_bld.o \ +mld_c_gs_solver_clone.o \ +mld_c_gs_solver_cnv.o \ +mld_c_gs_solver_dmp.o \ +mld_c_gs_solver_apply.o \ +mld_c_gs_solver_apply_vect.o \ mld_c_id_solver_apply.o \ mld_c_id_solver_apply_vect.o \ mld_c_id_solver_clone.o \ @@ -92,6 +98,12 @@ mld_s_diag_solver_apply_vect.o \ mld_s_diag_solver_bld.o \ mld_s_diag_solver_clone.o \ mld_s_diag_solver_cnv.o \ +mld_s_gs_solver_bld.o \ +mld_s_gs_solver_clone.o \ +mld_s_gs_solver_cnv.o \ +mld_s_gs_solver_dmp.o \ +mld_s_gs_solver_apply.o \ +mld_s_gs_solver_apply_vect.o \ mld_s_id_solver_apply.o \ mld_s_id_solver_apply_vect.o \ mld_s_id_solver_clone.o \ @@ -121,6 +133,12 @@ mld_z_diag_solver_apply_vect.o \ mld_z_diag_solver_bld.o \ mld_z_diag_solver_clone.o \ mld_z_diag_solver_cnv.o \ +mld_z_gs_solver_bld.o \ +mld_z_gs_solver_clone.o \ +mld_z_gs_solver_cnv.o \ +mld_z_gs_solver_dmp.o \ +mld_z_gs_solver_apply.o \ +mld_z_gs_solver_apply_vect.o \ mld_z_id_solver_apply.o \ mld_z_id_solver_apply_vect.o \ mld_z_id_solver_clone.o \ diff --git a/mlprec/impl/solver/mld_c_gs_solver_apply.f90 b/mlprec/impl/solver/mld_c_gs_solver_apply.f90 new file mode 100644 index 00000000..acd2366c --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_apply.f90 @@ -0,0 +1,196 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_gs_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + complex(psb_spk_), allocatable :: temp(:),wv(:),xit(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/4*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info) + call psb_geasb(xit,desc_data,info) + + select case(trans_) + case('N') + if (sv%eps <=szero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(cone,y,czero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(cone,x,czero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-cone,sv%u,xit,cone,wv,desc_data,info,doswap=.false.) + call psb_spsm(cone,sv%l,wv,czero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(cone,sv%u,x,czero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(cone,sv%u,x,czero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(cone,sv%dv,wv,czero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_gs_solver_apply diff --git a/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 new file mode 100644 index 00000000..6a8bfca2 --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_apply_vect.f90 @@ -0,0 +1,200 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_c_gs_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + type(psb_c_vect_type) :: wv, xit + complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + complex(psb_spk_), allocatable :: temp(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='c_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/itwo,n_row,izero,izero,izero/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/ithree,n_row,izero,izero,izero/)) + goto 9999 + end if + + + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) + + select case(trans_) + case('N') + if (sv%eps <=szero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(cone,y,czero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(cone,x,czero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-cone,sv%u,xit,cone,wv,desc_data,info,doswap=.false.) + call psb_spsm(cone,sv%l,wv,czero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(cone,sv%u,x,czero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(cone,sv%u,x,czero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(cone,sv%dv,wv,czero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + call wv%free(info) + call xit%free(info) + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_c_gs_solver_apply_vect diff --git a/mlprec/impl/solver/mld_c_gs_solver_bld.f90 b/mlprec/impl/solver/mld_c_gs_solver_bld.f90 new file mode 100644 index 00000000..62d9122b --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_bld.f90 @@ -0,0 +1,117 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_bld + + Implicit None + + ! Arguments + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='c_gs_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + + if (psb_toupper(upd) == 'F') then + nrow_a = a%get_nrows() + nztota = a%get_nzeros() +!!$ if (present(b)) then +!!$ nztota = nztota + b%get_nzeros() +!!$ end if + if (sv%eps <= dzero) then + ! + ! This cuts out the off-diagonal part, because it's supposed to + ! be handled by the outer Jacobi smoother. + ! + call a%tril(sv%l,info) + call a%triu(sv%u,info,diag=1,jmax=nrow_a) + + else + + info = psb_err_missing_override_method_ + call psb_errpush(info,name) + goto 9999 + end if + + + + call sv%l%set_asb() + call sv%l%trim() + call sv%u%set_asb() + call sv%u%trim() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_c_gs_solver_bld diff --git a/mlprec/impl/solver/mld_c_gs_solver_clone.f90 b/mlprec/impl/solver/mld_c_gs_solver_clone.f90 new file mode 100644 index 00000000..29af4419 --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_clone.f90 @@ -0,0 +1,88 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_clone(sv,svout,info) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_clone + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + class(mld_c_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + ! Local variables + integer(psb_ipk_) :: err_act + + + info=psb_success_ + call psb_erractionsave(err_act) + if (allocated(svout)) then + call svout%free(info) + if (info == psb_success_) deallocate(svout, stat=info) + end if + if (info == psb_success_) & + & allocate(mld_c_gs_solver_type :: svout, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + goto 9999 + end if + + select type(svo => svout) + type is (mld_c_gs_solver_type) + svo%sweeps = sv%sweeps + svo%eps = sv%eps + if (info == psb_success_) & + & call sv%l%clone(svo%l,info) + if (info == psb_success_) & + & call sv%u%clone(svo%u,info) + + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_c_gs_solver_clone diff --git a/mlprec/impl/solver/mld_c_gs_solver_cnv.f90 b/mlprec/impl/solver/mld_c_gs_solver_cnv.f90 new file mode 100644 index 00000000..4640b3e7 --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_cnv.f90 @@ -0,0 +1,75 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_cnv(sv,info,amold,vmold,imold) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_cnv + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: err_act, debug_unit, debug_level + character(len=20) :: name='c_gs_solver_cnv', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_c_gs_solver_cnv diff --git a/mlprec/impl/solver/mld_c_gs_solver_dmp.f90 b/mlprec/impl/solver/mld_c_gs_solver_dmp.f90 new file mode 100644 index 00000000..dd1b061a --- /dev/null +++ b/mlprec/impl/solver/mld_c_gs_solver_dmp.f90 @@ -0,0 +1,87 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_c_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + + use psb_base_mod + use mld_c_gs_solver, mld_protect_name => mld_c_gs_solver_dmp + implicit none + class(mld_c_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt,level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + integer(psb_ipk_) :: i, j, il1, iln, lname, lev + integer(psb_ipk_) :: icontxt,iam, np + character(len=80) :: prefix_ + character(len=120) :: fname ! len should be at least 20 more than + logical :: solver_ + ! len of prefix_ + + info = 0 + + + call psb_info(ictxt,iam,np) + + if (present(solver)) then + solver_ = solver + else + solver_ = .false. + end if + + if (solver_) then + if (present(prefix)) then + prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) + else + prefix_ = "dump_slv_d" + end if + lname = len_trim(prefix_) + fname = trim(prefix_) + write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam + lname = lname + 5 + + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_lower.mtx' + if (sv%l%is_asb()) & + & call sv%l%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_upper.mtx' + if (sv%u%is_asb()) & + & call sv%u%print(fname,head=head) + + end if + +end subroutine mld_c_gs_solver_dmp diff --git a/mlprec/impl/solver/mld_d_gs_solver_bld.f90 b/mlprec/impl/solver/mld_d_gs_solver_bld.f90 index 0b79b5df..b0b0f762 100644 --- a/mlprec/impl/solver/mld_d_gs_solver_bld.f90 +++ b/mlprec/impl/solver/mld_d_gs_solver_bld.f90 @@ -55,7 +55,6 @@ subroutine mld_d_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) class(psb_i_base_vect_type), intent(in), optional :: imold ! Local variables integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota -!!$ real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level character(len=20) :: name='d_gs_solver_bld', ch_err diff --git a/mlprec/impl/solver/mld_s_gs_solver_apply.f90 b/mlprec/impl/solver/mld_s_gs_solver_apply.f90 new file mode 100644 index 00000000..fa452a4d --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_apply.f90 @@ -0,0 +1,196 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_gs_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + real(psb_spk_), allocatable :: temp(:),wv(:),xit(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/4*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info) + call psb_geasb(xit,desc_data,info) + + select case(trans_) + case('N') + if (sv%eps <=szero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(sone,y,szero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(sone,x,szero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-sone,sv%u,xit,sone,wv,desc_data,info,doswap=.false.) + call psb_spsm(sone,sv%l,wv,szero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(sone,sv%u,x,szero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(sone,sv%u,x,szero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(sone,sv%dv,wv,szero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_gs_solver_apply diff --git a/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 new file mode 100644 index 00000000..71698d13 --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_apply_vect.f90 @@ -0,0 +1,200 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_s_gs_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + type(psb_s_vect_type) :: wv, xit + real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:) + real(psb_spk_), allocatable :: temp(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='s_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/itwo,n_row,izero,izero,izero/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/ithree,n_row,izero,izero,izero/)) + goto 9999 + end if + + + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) + + select case(trans_) + case('N') + if (sv%eps <=szero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(sone,y,szero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(sone,x,szero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-sone,sv%u,xit,sone,wv,desc_data,info,doswap=.false.) + call psb_spsm(sone,sv%l,wv,szero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(sone,sv%u,x,szero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(sone,sv%u,x,szero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(sone,sv%dv,wv,szero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + call wv%free(info) + call xit%free(info) + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_s_gs_solver_apply_vect diff --git a/mlprec/impl/solver/mld_s_gs_solver_bld.f90 b/mlprec/impl/solver/mld_s_gs_solver_bld.f90 new file mode 100644 index 00000000..91979c61 --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_bld.f90 @@ -0,0 +1,117 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_bld + + Implicit None + + ! Arguments + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_s_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='s_gs_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + + if (psb_toupper(upd) == 'F') then + nrow_a = a%get_nrows() + nztota = a%get_nzeros() +!!$ if (present(b)) then +!!$ nztota = nztota + b%get_nzeros() +!!$ end if + if (sv%eps <= dzero) then + ! + ! This cuts out the off-diagonal part, because it's supposed to + ! be handled by the outer Jacobi smoother. + ! + call a%tril(sv%l,info) + call a%triu(sv%u,info,diag=1,jmax=nrow_a) + + else + + info = psb_err_missing_override_method_ + call psb_errpush(info,name) + goto 9999 + end if + + + + call sv%l%set_asb() + call sv%l%trim() + call sv%u%set_asb() + call sv%u%trim() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_s_gs_solver_bld diff --git a/mlprec/impl/solver/mld_s_gs_solver_clone.f90 b/mlprec/impl/solver/mld_s_gs_solver_clone.f90 new file mode 100644 index 00000000..feeaeb33 --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_clone.f90 @@ -0,0 +1,88 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_clone(sv,svout,info) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_clone + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + class(mld_s_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + ! Local variables + integer(psb_ipk_) :: err_act + + + info=psb_success_ + call psb_erractionsave(err_act) + if (allocated(svout)) then + call svout%free(info) + if (info == psb_success_) deallocate(svout, stat=info) + end if + if (info == psb_success_) & + & allocate(mld_s_gs_solver_type :: svout, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + goto 9999 + end if + + select type(svo => svout) + type is (mld_s_gs_solver_type) + svo%sweeps = sv%sweeps + svo%eps = sv%eps + if (info == psb_success_) & + & call sv%l%clone(svo%l,info) + if (info == psb_success_) & + & call sv%u%clone(svo%u,info) + + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_s_gs_solver_clone diff --git a/mlprec/impl/solver/mld_s_gs_solver_cnv.f90 b/mlprec/impl/solver/mld_s_gs_solver_cnv.f90 new file mode 100644 index 00000000..152d4684 --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_cnv.f90 @@ -0,0 +1,75 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_cnv(sv,info,amold,vmold,imold) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_cnv + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: err_act, debug_unit, debug_level + character(len=20) :: name='s_gs_solver_cnv', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_s_gs_solver_cnv diff --git a/mlprec/impl/solver/mld_s_gs_solver_dmp.f90 b/mlprec/impl/solver/mld_s_gs_solver_dmp.f90 new file mode 100644 index 00000000..cd97cd25 --- /dev/null +++ b/mlprec/impl/solver/mld_s_gs_solver_dmp.f90 @@ -0,0 +1,87 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_s_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + + use psb_base_mod + use mld_s_gs_solver, mld_protect_name => mld_s_gs_solver_dmp + implicit none + class(mld_s_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt,level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + integer(psb_ipk_) :: i, j, il1, iln, lname, lev + integer(psb_ipk_) :: icontxt,iam, np + character(len=80) :: prefix_ + character(len=120) :: fname ! len should be at least 20 more than + logical :: solver_ + ! len of prefix_ + + info = 0 + + + call psb_info(ictxt,iam,np) + + if (present(solver)) then + solver_ = solver + else + solver_ = .false. + end if + + if (solver_) then + if (present(prefix)) then + prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) + else + prefix_ = "dump_slv_d" + end if + lname = len_trim(prefix_) + fname = trim(prefix_) + write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam + lname = lname + 5 + + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_lower.mtx' + if (sv%l%is_asb()) & + & call sv%l%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_upper.mtx' + if (sv%u%is_asb()) & + & call sv%u%print(fname,head=head) + + end if + +end subroutine mld_s_gs_solver_dmp diff --git a/mlprec/impl/solver/mld_z_gs_solver_apply.f90 b/mlprec/impl/solver/mld_z_gs_solver_apply.f90 new file mode 100644 index 00000000..e11bc09d --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_apply.f90 @@ -0,0 +1,196 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_apply + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_gs_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + complex(psb_dpk_), allocatable :: temp(:),wv(:),xit(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/4*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info) + call psb_geasb(xit,desc_data,info) + + select case(trans_) + case('N') + if (sv%eps <=dzero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(zone,y,zzero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(zone,x,zzero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-zone,sv%u,xit,zone,wv,desc_data,info,doswap=.false.) + call psb_spsm(zone,sv%l,wv,zzero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(zone,sv%u,x,zzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(zone,sv%u,x,zzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(zone,sv%dv,wv,zzero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_gs_solver_apply diff --git a/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 b/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 new file mode 100644 index 00000000..5140cdd7 --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_apply_vect.f90 @@ -0,0 +1,200 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_apply_vect + implicit none + type(psb_desc_type), intent(in) :: desc_data + class(mld_z_gs_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + + integer(psb_ipk_) :: n_row,n_col, itx + type(psb_z_vect_type) :: wv, xit + complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:) + complex(psb_dpk_), allocatable :: temp(:) + integer(psb_ipk_) :: ictxt,np,me,i, err_act + character :: trans_ + character(len=20) :: name='z_gs_solver_apply' + + call psb_erractionsave(err_act) + ictxt = desc_data%get_ctxt() + call psb_info(ictxt,me,np) + info = psb_success_ + + trans_ = psb_toupper(trans) + select case(trans_) + case('N') +!!$ case('T') +!!$ case('C') + case default + call psb_errpush(psb_err_iarg_invalid_i_,name) + goto 9999 + end select + + n_row = desc_data%get_local_rows() + n_col = desc_data%get_local_cols() + + + if (x%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/itwo,n_row,izero,izero,izero/)) + goto 9999 + end if + if (y%get_nrows() < n_row) then + info = 36 + call psb_errpush(info,name,& + & i_err=(/ithree,n_row,izero,izero,izero/)) + goto 9999 + end if + + + + if (n_col <= size(work)) then + ww => work(1:n_col) + if ((4*n_col+n_col) <= size(work)) then + aux => work(n_col+1:) + else + allocate(aux(4*n_col),stat=info) + endif + else + allocate(ww(n_col),aux(4*n_col),stat=info) + endif + + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,& + & i_err=(/5*n_col,izero,izero,izero,izero/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + + call psb_geasb(wv,desc_data,info,mold=x%v,scratch=.true.) + call psb_geasb(xit,desc_data,info,mold=x%v,scratch=.true.) + + select case(trans_) + case('N') + if (sv%eps <=dzero) then + ! + ! Fixed number of iterations + ! + ! + ! WARNING: this is not completely satisfactory. We are assuming here Y + ! as the initial guess, but this is only working if we are called from the + ! current JAC smoother loop. A good solution would be to have a separate + ! input argument as the initial guess + ! +!!$ write(0,*) 'GS Iteration with ',sv%sweeps + call psb_geaxpby(zone,y,zzero,xit,desc_data,info) + do itx=1,sv%sweeps + call psb_geaxpby(zone,x,zzero,wv,desc_data,info) + ! Update with U. The off-diagonal block is taken care + ! from the Jacobi smoother, hence this is purely local. + call psb_spmm(-zone,sv%u,xit,zone,wv,desc_data,info,doswap=.false.) + call psb_spsm(zone,sv%l,wv,zzero,xit,desc_data,info) +!!$ temp = xit%get_vect() +!!$ write(0,*) me,'GS Iteration ',itx,':',temp(1:n_row) + end do + + call psb_geaxpby(alpha,xit,beta,y,desc_data,info) + + else + ! + ! Iterations to convergence, not implemented right now. + ! + info = psb_err_internal_error_ + call psb_errpush(info,name,a_err='EPS>0 not implemented in GS subsolve') + goto 9999 + + end if +!!$ case('T') +!!$ call psb_spsm(zone,sv%u,x,zzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='L',diag=sv%dv,choice=psb_none_,work=aux) +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ case('C') +!!$ +!!$ call psb_spsm(zone,sv%u,x,zzero,wv,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) +!!$ +!!$ call wv1%mlt(zone,sv%dv,wv,zzero,info,conjgx=trans_) +!!$ +!!$ if (info == psb_success_) call psb_spsm(alpha,sv%l,wv1,beta,y,desc_data,info,& +!!$ & trans=trans_,scale='U',choice=psb_none_,work=aux) + + case default + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid TRANS in GS subsolve') + goto 9999 + end select + + + if (info /= psb_success_) then + + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in subsolve') + goto 9999 + endif + call wv%free(info) + call xit%free(info) + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then + else + deallocate(aux) + endif + else + deallocate(ww,aux) + endif + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return + +end subroutine mld_z_gs_solver_apply_vect diff --git a/mlprec/impl/solver/mld_z_gs_solver_bld.f90 b/mlprec/impl/solver/mld_z_gs_solver_bld.f90 new file mode 100644 index 00000000..e9dee83a --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_bld.f90 @@ -0,0 +1,117 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_bld + + Implicit None + + ! Arguments + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota + integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level + character(len=20) :: name='z_gs_solver_bld', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ictxt = desc_a%get_context() + call psb_info(ictxt, me, np) + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' start' + + + n_row = desc_a%get_local_rows() + + if (psb_toupper(upd) == 'F') then + nrow_a = a%get_nrows() + nztota = a%get_nzeros() +!!$ if (present(b)) then +!!$ nztota = nztota + b%get_nzeros() +!!$ end if + if (sv%eps <= dzero) then + ! + ! This cuts out the off-diagonal part, because it's supposed to + ! be handled by the outer Jacobi smoother. + ! + call a%tril(sv%l,info) + call a%triu(sv%u,info,diag=1,jmax=nrow_a) + + else + + info = psb_err_missing_override_method_ + call psb_errpush(info,name) + goto 9999 + end if + + + + call sv%l%set_asb() + call sv%l%trim() + call sv%u%set_asb() + call sv%u%trim() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) me,' ',trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_z_gs_solver_bld diff --git a/mlprec/impl/solver/mld_z_gs_solver_clone.f90 b/mlprec/impl/solver/mld_z_gs_solver_clone.f90 new file mode 100644 index 00000000..01f19fd8 --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_clone.f90 @@ -0,0 +1,88 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_clone(sv,svout,info) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_clone + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + class(mld_z_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + ! Local variables + integer(psb_ipk_) :: err_act + + + info=psb_success_ + call psb_erractionsave(err_act) + if (allocated(svout)) then + call svout%free(info) + if (info == psb_success_) deallocate(svout, stat=info) + end if + if (info == psb_success_) & + & allocate(mld_z_gs_solver_type :: svout, stat=info) + if (info /= 0) then + info = psb_err_alloc_dealloc_ + goto 9999 + end if + + select type(svo => svout) + type is (mld_z_gs_solver_type) + svo%sweeps = sv%sweeps + svo%eps = sv%eps + if (info == psb_success_) & + & call sv%l%clone(svo%l,info) + if (info == psb_success_) & + & call sv%u%clone(svo%u,info) + + class default + info = psb_err_internal_error_ + end select + + if (info /= 0) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_z_gs_solver_clone diff --git a/mlprec/impl/solver/mld_z_gs_solver_cnv.f90 b/mlprec/impl/solver/mld_z_gs_solver_cnv.f90 new file mode 100644 index 00000000..34c6ce39 --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_cnv.f90 @@ -0,0 +1,75 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_cnv(sv,info,amold,vmold,imold) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_cnv + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + ! Local variables + integer(psb_ipk_) :: err_act, debug_unit, debug_level + character(len=20) :: name='z_gs_solver_cnv', ch_err + + info=psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + + if (present(amold)) then + call sv%l%cscnv(info,mold=amold) + call sv%u%cscnv(info,mold=amold) + end if + + if (debug_level >= psb_debug_outer_) & + & write(debug_unit,*) trim(name),' end' + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + + return +end subroutine mld_z_gs_solver_cnv diff --git a/mlprec/impl/solver/mld_z_gs_solver_dmp.f90 b/mlprec/impl/solver/mld_z_gs_solver_dmp.f90 new file mode 100644 index 00000000..282801a8 --- /dev/null +++ b/mlprec/impl/solver/mld_z_gs_solver_dmp.f90 @@ -0,0 +1,87 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +subroutine mld_z_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + + use psb_base_mod + use mld_z_gs_solver, mld_protect_name => mld_z_gs_solver_dmp + implicit none + class(mld_z_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt,level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + integer(psb_ipk_) :: i, j, il1, iln, lname, lev + integer(psb_ipk_) :: icontxt,iam, np + character(len=80) :: prefix_ + character(len=120) :: fname ! len should be at least 20 more than + logical :: solver_ + ! len of prefix_ + + info = 0 + + + call psb_info(ictxt,iam,np) + + if (present(solver)) then + solver_ = solver + else + solver_ = .false. + end if + + if (solver_) then + if (present(prefix)) then + prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) + else + prefix_ = "dump_slv_d" + end if + lname = len_trim(prefix_) + fname = trim(prefix_) + write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam + lname = lname + 5 + + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_lower.mtx' + if (sv%l%is_asb()) & + & call sv%l%print(fname,head=head) + write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_upper.mtx' + if (sv%u%is_asb()) & + & call sv%u%print(fname,head=head) + + end if + +end subroutine mld_z_gs_solver_dmp diff --git a/mlprec/mld_c_base_solver_mod.f90 b/mlprec/mld_c_base_solver_mod.f90 index 8777716e..653b1530 100644 --- a/mlprec/mld_c_base_solver_mod.f90 +++ b/mlprec/mld_c_base_solver_mod.f90 @@ -107,10 +107,12 @@ module mld_c_base_solver_mod procedure, pass(sv) :: get_nzeros => c_base_solver_get_nzeros procedure, nopass :: stringval => mld_stringval procedure, nopass :: get_fmt => c_base_solver_get_fmt + procedure, nopass :: is_iterative => c_base_solver_is_iterative end type mld_c_base_solver_type private :: c_base_solver_sizeof, c_base_solver_default,& - & c_base_solver_get_nzeros, c_base_solver_get_fmt + & c_base_solver_get_nzeros, c_base_solver_get_fmt, & + & c_base_solver_is_iterative interface @@ -386,5 +388,15 @@ contains val = "Base solver" end function c_base_solver_get_fmt + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function c_base_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .false. + end function c_base_solver_is_iterative end module mld_c_base_solver_mod diff --git a/mlprec/mld_c_gs_solver.f90 b/mlprec/mld_c_gs_solver.f90 new file mode 100644 index 00000000..7b074d78 --- /dev/null +++ b/mlprec/mld_c_gs_solver.f90 @@ -0,0 +1,518 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! +! +! +! +! + +module mld_c_gs_solver + + use mld_c_base_solver_mod + + type, extends(mld_c_base_solver_type) :: mld_c_gs_solver_type + type(psb_cspmat_type) :: l, u + integer(psb_ipk_) :: sweeps + real(psb_spk_) :: eps + contains + procedure, pass(sv) :: dump => mld_c_gs_solver_dmp + procedure, pass(sv) :: ccheck => d_gs_solver_check + procedure, pass(sv) :: clone => mld_c_gs_solver_clone + procedure, pass(sv) :: build => mld_c_gs_solver_bld + procedure, pass(sv) :: cnv => mld_c_gs_solver_cnv + procedure, pass(sv) :: apply_v => mld_c_gs_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_c_gs_solver_apply + procedure, pass(sv) :: free => d_gs_solver_free + procedure, pass(sv) :: seti => d_gs_solver_seti + procedure, pass(sv) :: setc => d_gs_solver_setc + procedure, pass(sv) :: setr => d_gs_solver_setr + procedure, pass(sv) :: cseti => d_gs_solver_cseti + procedure, pass(sv) :: csetc => d_gs_solver_csetc + procedure, pass(sv) :: csetr => d_gs_solver_csetr + procedure, pass(sv) :: descr => d_gs_solver_descr + procedure, pass(sv) :: default => d_gs_solver_default + procedure, pass(sv) :: sizeof => d_gs_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros + procedure, nopass :: get_fmt => d_gs_solver_get_fmt + procedure, nopass :: is_iterative => d_gs_solver_is_iterative + end type mld_c_gs_solver_type + + + private :: d_gs_solver_bld, d_gs_solver_apply, & + & d_gs_solver_free, d_gs_solver_seti, & + & d_gs_solver_setc, d_gs_solver_setr,& + & d_gs_solver_descr, d_gs_solver_sizeof, & + & d_gs_solver_default, d_gs_solver_dmp, & + & d_gs_solver_apply_vect, d_gs_solver_get_nzeros, & + & d_gs_solver_get_fmt, d_gs_solver_check,& + & d_gs_solver_is_iterative + + + interface + subroutine mld_c_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, & + & 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_gs_solver_type), intent(inout) :: sv + type(psb_c_vect_type),intent(inout) :: x + type(psb_c_vect_type),intent(inout) :: y + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_gs_solver_apply_vect + end interface + + interface + subroutine mld_c_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, & + & 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_gs_solver_type), intent(inout) :: sv + complex(psb_spk_),intent(inout) :: x(:) + complex(psb_spk_),intent(inout) :: y(:) + complex(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_gs_solver_apply + end interface + + interface + subroutine mld_c_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_cspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_c_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_cspmat_type), intent(in), target, optional :: b + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_c_gs_solver_bld + end interface + + interface + subroutine mld_c_gs_solver_cnv(sv,info,amold,vmold,imold) + import :: mld_c_gs_solver_type, psb_spk_, & + & psb_c_base_sparse_mat, psb_c_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_c_base_sparse_mat), intent(in), optional :: amold + class(psb_c_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_c_gs_solver_cnv + end interface + + interface + subroutine mld_c_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + import :: psb_desc_type, mld_c_gs_solver_type, psb_c_vect_type, psb_spk_, & + & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type, & + & psb_ipk_ + implicit none + class(mld_c_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + end subroutine mld_c_gs_solver_dmp + end interface + + interface + subroutine mld_c_gs_solver_clone(sv,svout,info) + import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & + & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & + & mld_c_base_solver_type, mld_c_gs_solver_type, psb_ipk_ + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + class(mld_c_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + end subroutine mld_c_gs_solver_clone + end interface + +contains + + subroutine d_gs_solver_default(sv) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + + sv%sweeps = ione + sv%eps = dzero + + return + end subroutine d_gs_solver_default + + subroutine d_gs_solver_check(sv,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + call mld_check_def(sv%sweeps,& + & 'GS Sweeps',ione,is_int_positive) + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine d_gs_solver_check + + + subroutine d_gs_solver_seti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_seti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(what) + case(mld_solver_sweeps_) + sv%sweeps = val + case default + call sv%mld_c_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_seti + + subroutine d_gs_solver_setc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_setc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setc + + subroutine d_gs_solver_setr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_setr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(what) + case(mld_solver_eps_) + sv%eps = val + case default + call sv%mld_c_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setr + + subroutine d_gs_solver_cseti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('SOLVER_SWEEPS') + sv%sweeps = val + case default + call sv%mld_c_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_cseti + + subroutine d_gs_solver_csetc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetc + + subroutine d_gs_solver_csetr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('SOLVER_EPS') + sv%eps = val + case default + call sv%mld_c_base_solver_type%set(what,val,info) + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetr + + subroutine d_gs_solver_free(sv,info) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%l%free() + call sv%u%free() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_free + + subroutine d_gs_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(mld_c_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='mld_c_gs_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + if (sv%eps<=dzero) then + write(iout_,*) ' Gauss-Seidel iterative solver with ',& + & sv%sweeps,' sweeps' + else + write(iout_,*) ' Gauss-Seidel iterative solver with tolerance',& + & sv%eps,' and maxit', sv%sweeps + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_descr + + function d_gs_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(mld_c_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = 0 + val = val + sv%l%get_nzeros() + val = val + sv%u%get_nzeros() + + return + end function d_gs_solver_get_nzeros + + function d_gs_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(mld_c_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = psb_sizeof_int + val = val + sv%l%sizeof() + val = val + sv%u%sizeof() + + return + end function d_gs_solver_sizeof + + function d_gs_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Gauss-Seidel solver" + end function d_gs_solver_get_fmt + + + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function d_gs_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function d_gs_solver_is_iterative + +end module mld_c_gs_solver diff --git a/mlprec/mld_c_prec_mod.f90 b/mlprec/mld_c_prec_mod.f90 index 2d2839c0..b5e7e0bb 100644 --- a/mlprec/mld_c_prec_mod.f90 +++ b/mlprec/mld_c_prec_mod.f90 @@ -51,6 +51,7 @@ module mld_c_prec_mod use mld_c_id_solver use mld_c_diag_solver use mld_c_ilu_solver + use mld_c_gs_solver interface mld_precinit subroutine mld_cprecinit(p,ptype,info,nlev) diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index 18bec836..6c61d95e 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -471,7 +471,7 @@ contains call p%precv(1)%sm%descr(info,iout=iout_) if (nlev == 1) then if (p%precv(1)%parms%sweeps > 1) then - write(iout_,*) ' Number of sweeps : ',& + write(iout_,*) ' Number of smoother sweeps : ',& & p%precv(1)%parms%sweeps end if write(iout_,*) diff --git a/mlprec/mld_d_base_solver_mod.f90 b/mlprec/mld_d_base_solver_mod.f90 index be2c1b0c..89f6cd3c 100644 --- a/mlprec/mld_d_base_solver_mod.f90 +++ b/mlprec/mld_d_base_solver_mod.f90 @@ -111,8 +111,8 @@ module mld_d_base_solver_mod end type mld_d_base_solver_type private :: d_base_solver_sizeof, d_base_solver_default,& - & d_base_solver_get_nzeros, d_base_solver_get_fmt, & - & d_base_solver_is_iterative + & d_base_solver_get_nzeros, d_base_solver_get_fmt, & + & d_base_solver_is_iterative interface diff --git a/mlprec/mld_s_base_solver_mod.f90 b/mlprec/mld_s_base_solver_mod.f90 index e892a3f6..ef281740 100644 --- a/mlprec/mld_s_base_solver_mod.f90 +++ b/mlprec/mld_s_base_solver_mod.f90 @@ -107,10 +107,12 @@ module mld_s_base_solver_mod procedure, pass(sv) :: get_nzeros => s_base_solver_get_nzeros procedure, nopass :: stringval => mld_stringval procedure, nopass :: get_fmt => s_base_solver_get_fmt + procedure, nopass :: is_iterative => s_base_solver_is_iterative end type mld_s_base_solver_type private :: s_base_solver_sizeof, s_base_solver_default,& - & s_base_solver_get_nzeros, s_base_solver_get_fmt + & s_base_solver_get_nzeros, s_base_solver_get_fmt, & + & s_base_solver_is_iterative interface @@ -386,5 +388,15 @@ contains val = "Base solver" end function s_base_solver_get_fmt + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function s_base_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .false. + end function s_base_solver_is_iterative end module mld_s_base_solver_mod diff --git a/mlprec/mld_s_gs_solver.f90 b/mlprec/mld_s_gs_solver.f90 new file mode 100644 index 00000000..7e040560 --- /dev/null +++ b/mlprec/mld_s_gs_solver.f90 @@ -0,0 +1,518 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! +! +! +! +! + +module mld_s_gs_solver + + use mld_s_base_solver_mod + + type, extends(mld_s_base_solver_type) :: mld_s_gs_solver_type + type(psb_sspmat_type) :: l, u + integer(psb_ipk_) :: sweeps + real(psb_spk_) :: eps + contains + procedure, pass(sv) :: dump => mld_s_gs_solver_dmp + procedure, pass(sv) :: ccheck => d_gs_solver_check + procedure, pass(sv) :: clone => mld_s_gs_solver_clone + procedure, pass(sv) :: build => mld_s_gs_solver_bld + procedure, pass(sv) :: cnv => mld_s_gs_solver_cnv + procedure, pass(sv) :: apply_v => mld_s_gs_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_s_gs_solver_apply + procedure, pass(sv) :: free => d_gs_solver_free + procedure, pass(sv) :: seti => d_gs_solver_seti + procedure, pass(sv) :: setc => d_gs_solver_setc + procedure, pass(sv) :: setr => d_gs_solver_setr + procedure, pass(sv) :: cseti => d_gs_solver_cseti + procedure, pass(sv) :: csetc => d_gs_solver_csetc + procedure, pass(sv) :: csetr => d_gs_solver_csetr + procedure, pass(sv) :: descr => d_gs_solver_descr + procedure, pass(sv) :: default => d_gs_solver_default + procedure, pass(sv) :: sizeof => d_gs_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros + procedure, nopass :: get_fmt => d_gs_solver_get_fmt + procedure, nopass :: is_iterative => d_gs_solver_is_iterative + end type mld_s_gs_solver_type + + + private :: d_gs_solver_bld, d_gs_solver_apply, & + & d_gs_solver_free, d_gs_solver_seti, & + & d_gs_solver_setc, d_gs_solver_setr,& + & d_gs_solver_descr, d_gs_solver_sizeof, & + & d_gs_solver_default, d_gs_solver_dmp, & + & d_gs_solver_apply_vect, d_gs_solver_get_nzeros, & + & d_gs_solver_get_fmt, d_gs_solver_check,& + & d_gs_solver_is_iterative + + + interface + subroutine mld_s_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_s_gs_solver_type, psb_s_vect_type, psb_spk_, & + & 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_gs_solver_type), intent(inout) :: sv + type(psb_s_vect_type),intent(inout) :: x + type(psb_s_vect_type),intent(inout) :: y + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_gs_solver_apply_vect + end interface + + interface + subroutine mld_s_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_s_gs_solver_type, psb_s_vect_type, psb_spk_, & + & 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_gs_solver_type), intent(inout) :: sv + real(psb_spk_),intent(inout) :: x(:) + real(psb_spk_),intent(inout) :: y(:) + real(psb_spk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + real(psb_spk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_gs_solver_apply + end interface + + interface + subroutine mld_s_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + import :: psb_desc_type, mld_s_gs_solver_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_sspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_s_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_sspmat_type), intent(in), target, optional :: b + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_s_gs_solver_bld + end interface + + interface + subroutine mld_s_gs_solver_cnv(sv,info,amold,vmold,imold) + import :: mld_s_gs_solver_type, psb_spk_, & + & psb_s_base_sparse_mat, psb_s_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_s_base_sparse_mat), intent(in), optional :: amold + class(psb_s_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_s_gs_solver_cnv + end interface + + interface + subroutine mld_s_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + import :: psb_desc_type, mld_s_gs_solver_type, psb_s_vect_type, psb_spk_, & + & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type, & + & psb_ipk_ + implicit none + class(mld_s_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + end subroutine mld_s_gs_solver_dmp + end interface + + interface + subroutine mld_s_gs_solver_clone(sv,svout,info) + import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & + & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & + & mld_s_base_solver_type, mld_s_gs_solver_type, psb_ipk_ + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + class(mld_s_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + end subroutine mld_s_gs_solver_clone + end interface + +contains + + subroutine d_gs_solver_default(sv) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + + sv%sweeps = ione + sv%eps = dzero + + return + end subroutine d_gs_solver_default + + subroutine d_gs_solver_check(sv,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + call mld_check_def(sv%sweeps,& + & 'GS Sweeps',ione,is_int_positive) + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine d_gs_solver_check + + + subroutine d_gs_solver_seti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_seti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(what) + case(mld_solver_sweeps_) + sv%sweeps = val + case default + call sv%mld_s_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_seti + + subroutine d_gs_solver_setc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_setc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setc + + subroutine d_gs_solver_setr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_setr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(what) + case(mld_solver_eps_) + sv%eps = val + case default + call sv%mld_s_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setr + + subroutine d_gs_solver_cseti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('SOLVER_SWEEPS') + sv%sweeps = val + case default + call sv%mld_s_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_cseti + + subroutine d_gs_solver_csetc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetc + + subroutine d_gs_solver_csetr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('SOLVER_EPS') + sv%eps = val + case default + call sv%mld_s_base_solver_type%set(what,val,info) + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetr + + subroutine d_gs_solver_free(sv,info) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%l%free() + call sv%u%free() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_free + + subroutine d_gs_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(mld_s_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='mld_s_gs_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + if (sv%eps<=dzero) then + write(iout_,*) ' Gauss-Seidel iterative solver with ',& + & sv%sweeps,' sweeps' + else + write(iout_,*) ' Gauss-Seidel iterative solver with tolerance',& + & sv%eps,' and maxit', sv%sweeps + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_descr + + function d_gs_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(mld_s_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = 0 + val = val + sv%l%get_nzeros() + val = val + sv%u%get_nzeros() + + return + end function d_gs_solver_get_nzeros + + function d_gs_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(mld_s_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = psb_sizeof_int + val = val + sv%l%sizeof() + val = val + sv%u%sizeof() + + return + end function d_gs_solver_sizeof + + function d_gs_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Gauss-Seidel solver" + end function d_gs_solver_get_fmt + + + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function d_gs_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function d_gs_solver_is_iterative + +end module mld_s_gs_solver diff --git a/mlprec/mld_s_prec_mod.f90 b/mlprec/mld_s_prec_mod.f90 index 98e8831a..c8c0e4b2 100644 --- a/mlprec/mld_s_prec_mod.f90 +++ b/mlprec/mld_s_prec_mod.f90 @@ -51,6 +51,7 @@ module mld_s_prec_mod use mld_s_id_solver use mld_s_diag_solver use mld_s_ilu_solver + use mld_s_gs_solver interface mld_precinit subroutine mld_sprecinit(p,ptype,info,nlev) diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index 299d7a04..3b1d30f7 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -471,7 +471,7 @@ contains call p%precv(1)%sm%descr(info,iout=iout_) if (nlev == 1) then if (p%precv(1)%parms%sweeps > 1) then - write(iout_,*) ' Number of sweeps : ',& + write(iout_,*) ' Number of smoother sweeps : ',& & p%precv(1)%parms%sweeps end if write(iout_,*) diff --git a/mlprec/mld_z_base_solver_mod.f90 b/mlprec/mld_z_base_solver_mod.f90 index 9871f521..9e2d1223 100644 --- a/mlprec/mld_z_base_solver_mod.f90 +++ b/mlprec/mld_z_base_solver_mod.f90 @@ -107,10 +107,12 @@ module mld_z_base_solver_mod procedure, pass(sv) :: get_nzeros => z_base_solver_get_nzeros procedure, nopass :: stringval => mld_stringval procedure, nopass :: get_fmt => z_base_solver_get_fmt + procedure, nopass :: is_iterative => z_base_solver_is_iterative end type mld_z_base_solver_type private :: z_base_solver_sizeof, z_base_solver_default,& - & z_base_solver_get_nzeros, z_base_solver_get_fmt + & z_base_solver_get_nzeros, z_base_solver_get_fmt, & + & z_base_solver_is_iterative interface @@ -386,5 +388,15 @@ contains val = "Base solver" end function z_base_solver_get_fmt + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function z_base_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .false. + end function z_base_solver_is_iterative end module mld_z_base_solver_mod diff --git a/mlprec/mld_z_gs_solver.f90 b/mlprec/mld_z_gs_solver.f90 new file mode 100644 index 00000000..1bbac45e --- /dev/null +++ b/mlprec/mld_z_gs_solver.f90 @@ -0,0 +1,518 @@ +!!$ +!!$ +!!$ MLD2P4 version 2.0 +!!$ MultiLevel Domain Decomposition Parallel Preconditioners Package +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) +!!$ +!!$ (C) Copyright 2008, 2010, 2012, 2015 +!!$ +!!$ Salvatore Filippone University of Rome Tor Vergata +!!$ Alfredo Buttari CNRS-IRIT, Toulouse +!!$ Pasqua D'Ambra ICAR-CNR, Naples +!!$ Daniela di Serafino Second University of Naples +!!$ +!!$ Redistribution and use in source and binary forms, with or without +!!$ modification, are permitted provided that the following conditions +!!$ are met: +!!$ 1. Redistributions of source code must retain the above copyright +!!$ notice, this list of conditions and the following disclaimer. +!!$ 2. Redistributions in binary form must reproduce the above copyright +!!$ notice, this list of conditions, and the following disclaimer in the +!!$ documentation and/or other materials provided with the distribution. +!!$ 3. The name of the MLD2P4 group or the names of its contributors may +!!$ not be used to endorse or promote products derived from this +!!$ software without specific written permission. +!!$ +!!$ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +!!$ ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +!!$ TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +!!$ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE MLD2P4 GROUP OR ITS CONTRIBUTORS +!!$ BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +!!$ CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +!!$ SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +!!$ INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +!!$ CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +!!$ ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +!!$ POSSIBILITY OF SUCH DAMAGE. +!!$ +!!$ +! +! +! +! +! +! + +module mld_z_gs_solver + + use mld_z_base_solver_mod + + type, extends(mld_z_base_solver_type) :: mld_z_gs_solver_type + type(psb_zspmat_type) :: l, u + integer(psb_ipk_) :: sweeps + real(psb_dpk_) :: eps + contains + procedure, pass(sv) :: dump => mld_z_gs_solver_dmp + procedure, pass(sv) :: ccheck => d_gs_solver_check + procedure, pass(sv) :: clone => mld_z_gs_solver_clone + procedure, pass(sv) :: build => mld_z_gs_solver_bld + procedure, pass(sv) :: cnv => mld_z_gs_solver_cnv + procedure, pass(sv) :: apply_v => mld_z_gs_solver_apply_vect + procedure, pass(sv) :: apply_a => mld_z_gs_solver_apply + procedure, pass(sv) :: free => d_gs_solver_free + procedure, pass(sv) :: seti => d_gs_solver_seti + procedure, pass(sv) :: setc => d_gs_solver_setc + procedure, pass(sv) :: setr => d_gs_solver_setr + procedure, pass(sv) :: cseti => d_gs_solver_cseti + procedure, pass(sv) :: csetc => d_gs_solver_csetc + procedure, pass(sv) :: csetr => d_gs_solver_csetr + procedure, pass(sv) :: descr => d_gs_solver_descr + procedure, pass(sv) :: default => d_gs_solver_default + procedure, pass(sv) :: sizeof => d_gs_solver_sizeof + procedure, pass(sv) :: get_nzeros => d_gs_solver_get_nzeros + procedure, nopass :: get_fmt => d_gs_solver_get_fmt + procedure, nopass :: is_iterative => d_gs_solver_is_iterative + end type mld_z_gs_solver_type + + + private :: d_gs_solver_bld, d_gs_solver_apply, & + & d_gs_solver_free, d_gs_solver_seti, & + & d_gs_solver_setc, d_gs_solver_setr,& + & d_gs_solver_descr, d_gs_solver_sizeof, & + & d_gs_solver_default, d_gs_solver_dmp, & + & d_gs_solver_apply_vect, d_gs_solver_get_nzeros, & + & d_gs_solver_get_fmt, d_gs_solver_check,& + & d_gs_solver_is_iterative + + + interface + subroutine mld_z_gs_solver_apply_vect(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_z_gs_solver_type, psb_z_vect_type, psb_dpk_, & + & 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_gs_solver_type), intent(inout) :: sv + type(psb_z_vect_type),intent(inout) :: x + type(psb_z_vect_type),intent(inout) :: y + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_gs_solver_apply_vect + end interface + + interface + subroutine mld_z_gs_solver_apply(alpha,sv,x,beta,y,desc_data,trans,work,info) + import :: psb_desc_type, mld_z_gs_solver_type, psb_z_vect_type, psb_dpk_, & + & 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_gs_solver_type), intent(inout) :: sv + complex(psb_dpk_),intent(inout) :: x(:) + complex(psb_dpk_),intent(inout) :: y(:) + complex(psb_dpk_),intent(in) :: alpha,beta + character(len=1),intent(in) :: trans + complex(psb_dpk_),target, intent(inout) :: work(:) + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_gs_solver_apply + end interface + + interface + subroutine mld_z_gs_solver_bld(a,desc_a,sv,upd,info,b,amold,vmold,imold) + import :: psb_desc_type, mld_z_gs_solver_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + type(psb_zspmat_type), intent(in), target :: a + Type(psb_desc_type), Intent(in) :: desc_a + class(mld_z_gs_solver_type), intent(inout) :: sv + character, intent(in) :: upd + integer(psb_ipk_), intent(out) :: info + type(psb_zspmat_type), intent(in), target, optional :: b + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_z_gs_solver_bld + end interface + + interface + subroutine mld_z_gs_solver_cnv(sv,info,amold,vmold,imold) + import :: mld_z_gs_solver_type, psb_dpk_, & + & psb_z_base_sparse_mat, psb_z_base_vect_type,& + & psb_ipk_, psb_i_base_vect_type + implicit none + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + class(psb_z_base_sparse_mat), intent(in), optional :: amold + class(psb_z_base_vect_type), intent(in), optional :: vmold + class(psb_i_base_vect_type), intent(in), optional :: imold + end subroutine mld_z_gs_solver_cnv + end interface + + interface + subroutine mld_z_gs_solver_dmp(sv,ictxt,level,info,prefix,head,solver) + import :: psb_desc_type, mld_z_gs_solver_type, psb_z_vect_type, psb_dpk_, & + & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type, & + & psb_ipk_ + implicit none + class(mld_z_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(in) :: ictxt + integer(psb_ipk_), intent(in) :: level + integer(psb_ipk_), intent(out) :: info + character(len=*), intent(in), optional :: prefix, head + logical, optional, intent(in) :: solver + end subroutine mld_z_gs_solver_dmp + end interface + + interface + subroutine mld_z_gs_solver_clone(sv,svout,info) + import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & + & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & + & mld_z_base_solver_type, mld_z_gs_solver_type, psb_ipk_ + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + class(mld_z_base_solver_type), allocatable, intent(inout) :: svout + integer(psb_ipk_), intent(out) :: info + end subroutine mld_z_gs_solver_clone + end interface + +contains + + subroutine d_gs_solver_default(sv) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + + sv%sweeps = ione + sv%eps = dzero + + return + end subroutine d_gs_solver_default + + subroutine d_gs_solver_check(sv,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_check' + + call psb_erractionsave(err_act) + info = psb_success_ + + call mld_check_def(sv%sweeps,& + & 'GS Sweeps',ione,is_int_positive) + + if (info /= psb_success_) goto 9999 + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine d_gs_solver_check + + + subroutine d_gs_solver_seti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_seti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(what) + case(mld_solver_sweeps_) + sv%sweeps = val + case default + call sv%mld_z_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_seti + + subroutine d_gs_solver_setc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_setc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setc + + subroutine d_gs_solver_setr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_setr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(what) + case(mld_solver_eps_) + sv%eps = val + case default + call sv%mld_z_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_setr + + subroutine d_gs_solver_cseti(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_cseti' + + info = psb_success_ + call psb_erractionsave(err_act) + + select case(psb_toupper(what)) + case('SOLVER_SWEEPS') + sv%sweeps = val + case default + call sv%mld_z_base_solver_type%set(what,val,info) + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_cseti + + subroutine d_gs_solver_csetc(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + character(len=*), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, ival + character(len=20) :: name='d_gs_solver_csetc' + + info = psb_success_ + call psb_erractionsave(err_act) + + + ival = sv%stringval(val) + if (ival >= 0) then + call sv%set(what,ival,info) + end if + + if (info /= psb_success_) then + info = psb_err_from_subroutine_ + call psb_errpush(info, name) + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetc + + subroutine d_gs_solver_csetr(sv,what,val,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + character(len=*), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_csetr' + + call psb_erractionsave(err_act) + info = psb_success_ + + select case(psb_toupper(what)) + case('SOLVER_EPS') + sv%eps = val + case default + call sv%mld_z_base_solver_type%set(what,val,info) + end select + + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_csetr + + subroutine d_gs_solver_free(sv,info) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(inout) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act + character(len=20) :: name='d_gs_solver_free' + + call psb_erractionsave(err_act) + info = psb_success_ + + call sv%l%free() + call sv%u%free() + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_free + + subroutine d_gs_solver_descr(sv,info,iout,coarse) + + Implicit None + + ! Arguments + class(mld_z_gs_solver_type), intent(in) :: sv + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_), intent(in), optional :: iout + logical, intent(in), optional :: coarse + + ! Local variables + integer(psb_ipk_) :: err_act + character(len=20), parameter :: name='mld_z_gs_solver_descr' + integer(psb_ipk_) :: iout_ + + call psb_erractionsave(err_act) + info = psb_success_ + if (present(iout)) then + iout_ = iout + else + iout_ = 6 + endif + + if (sv%eps<=dzero) then + write(iout_,*) ' Gauss-Seidel iterative solver with ',& + & sv%sweeps,' sweeps' + else + write(iout_,*) ' Gauss-Seidel iterative solver with tolerance',& + & sv%eps,' and maxit', sv%sweeps + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + end subroutine d_gs_solver_descr + + function d_gs_solver_get_nzeros(sv) result(val) + + implicit none + ! Arguments + class(mld_z_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = 0 + val = val + sv%l%get_nzeros() + val = val + sv%u%get_nzeros() + + return + end function d_gs_solver_get_nzeros + + function d_gs_solver_sizeof(sv) result(val) + + implicit none + ! Arguments + class(mld_z_gs_solver_type), intent(in) :: sv + integer(psb_long_int_k_) :: val + integer(psb_ipk_) :: i + + val = psb_sizeof_int + val = val + sv%l%sizeof() + val = val + sv%u%sizeof() + + return + end function d_gs_solver_sizeof + + function d_gs_solver_get_fmt() result(val) + implicit none + character(len=32) :: val + + val = "Gauss-Seidel solver" + end function d_gs_solver_get_fmt + + + ! + ! If this is true, then the solver needs a starting + ! guess. Currently only handled in JAC smoother. + ! + function d_gs_solver_is_iterative() result(val) + implicit none + logical :: val + + val = .true. + end function d_gs_solver_is_iterative + +end module mld_z_gs_solver diff --git a/mlprec/mld_z_prec_mod.f90 b/mlprec/mld_z_prec_mod.f90 index 80be7a79..b1259522 100644 --- a/mlprec/mld_z_prec_mod.f90 +++ b/mlprec/mld_z_prec_mod.f90 @@ -51,6 +51,7 @@ module mld_z_prec_mod use mld_z_id_solver use mld_z_diag_solver use mld_z_ilu_solver + use mld_z_gs_solver interface mld_precinit subroutine mld_zprecinit(p,ptype,info,nlev) diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index 3e9c955b..0eda8dd0 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -471,7 +471,7 @@ contains call p%precv(1)%sm%descr(info,iout=iout_) if (nlev == 1) then if (p%precv(1)%parms%sweeps > 1) then - write(iout_,*) ' Number of sweeps : ',& + write(iout_,*) ' Number of smoother sweeps : ',& & p%precv(1)%parms%sweeps end if write(iout_,*)