diff --git a/base/modules/serial/psb_c_mat_mod.F90 b/base/modules/serial/psb_c_mat_mod.F90 index 76225758..5e889da2 100644 --- a/base/modules/serial/psb_c_mat_mod.F90 +++ b/base/modules/serial/psb_c_mat_mod.F90 @@ -128,6 +128,7 @@ module psb_c_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_c_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_c_free procedure, pass(a) :: trim => psb_c_trim procedure, pass(a) :: csput_a => psb_c_csput_a @@ -326,6 +327,7 @@ module psb_c_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_lc_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_lc_free procedure, pass(a) :: trim => psb_lc_trim procedure, pass(a) :: csput_a => psb_lc_csput_a @@ -604,12 +606,14 @@ module psb_c_mat_mod end interface interface - subroutine psb_c_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_cspmat_type + subroutine psb_c_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_cspmat_type, psb_c_base_sparse_mat class(psb_cspmat_type), intent(inout) :: a integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_c_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_c_csall end interface @@ -1384,12 +1388,14 @@ module psb_c_mat_mod end interface interface - subroutine psb_lc_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_lcspmat_type + subroutine psb_lc_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_lcspmat_type, psb_lc_base_sparse_mat class(psb_lcspmat_type), intent(inout) :: a integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_lc_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_lc_csall end interface diff --git a/base/modules/serial/psb_d_mat_mod.F90 b/base/modules/serial/psb_d_mat_mod.F90 index 878d099f..caf03994 100644 --- a/base/modules/serial/psb_d_mat_mod.F90 +++ b/base/modules/serial/psb_d_mat_mod.F90 @@ -128,6 +128,7 @@ module psb_d_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_d_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_d_free procedure, pass(a) :: trim => psb_d_trim procedure, pass(a) :: csput_a => psb_d_csput_a @@ -326,6 +327,7 @@ module psb_d_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_ld_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_ld_free procedure, pass(a) :: trim => psb_ld_trim procedure, pass(a) :: csput_a => psb_ld_csput_a @@ -604,12 +606,14 @@ module psb_d_mat_mod end interface interface - subroutine psb_d_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_dspmat_type + subroutine psb_d_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_dspmat_type, psb_d_base_sparse_mat class(psb_dspmat_type), intent(inout) :: a integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_d_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_d_csall end interface @@ -1384,12 +1388,14 @@ module psb_d_mat_mod end interface interface - subroutine psb_ld_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_ldspmat_type + subroutine psb_ld_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_ldspmat_type, psb_ld_base_sparse_mat class(psb_ldspmat_type), intent(inout) :: a integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_ld_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_ld_csall end interface diff --git a/base/modules/serial/psb_s_mat_mod.F90 b/base/modules/serial/psb_s_mat_mod.F90 index 3553e96b..8e3934b8 100644 --- a/base/modules/serial/psb_s_mat_mod.F90 +++ b/base/modules/serial/psb_s_mat_mod.F90 @@ -128,6 +128,7 @@ module psb_s_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_s_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_s_free procedure, pass(a) :: trim => psb_s_trim procedure, pass(a) :: csput_a => psb_s_csput_a @@ -326,6 +327,7 @@ module psb_s_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_ls_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_ls_free procedure, pass(a) :: trim => psb_ls_trim procedure, pass(a) :: csput_a => psb_ls_csput_a @@ -604,12 +606,14 @@ module psb_s_mat_mod end interface interface - subroutine psb_s_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_sspmat_type + subroutine psb_s_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_sspmat_type, psb_s_base_sparse_mat class(psb_sspmat_type), intent(inout) :: a integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_s_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_s_csall end interface @@ -1384,12 +1388,14 @@ module psb_s_mat_mod end interface interface - subroutine psb_ls_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_lsspmat_type + subroutine psb_ls_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_lsspmat_type, psb_ls_base_sparse_mat class(psb_lsspmat_type), intent(inout) :: a integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_ls_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_ls_csall end interface diff --git a/base/modules/serial/psb_z_mat_mod.F90 b/base/modules/serial/psb_z_mat_mod.F90 index 35586b3e..ed3338f9 100644 --- a/base/modules/serial/psb_z_mat_mod.F90 +++ b/base/modules/serial/psb_z_mat_mod.F90 @@ -128,6 +128,7 @@ module psb_z_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_z_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_z_free procedure, pass(a) :: trim => psb_z_trim procedure, pass(a) :: csput_a => psb_z_csput_a @@ -326,6 +327,7 @@ module psb_z_mat_mod ! Memory/data management procedure, pass(a) :: csall => psb_lz_csall + generic, public :: allocate => csall procedure, pass(a) :: free => psb_lz_free procedure, pass(a) :: trim => psb_lz_trim procedure, pass(a) :: csput_a => psb_lz_csput_a @@ -604,12 +606,14 @@ module psb_z_mat_mod end interface interface - subroutine psb_z_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_zspmat_type + subroutine psb_z_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_zspmat_type, psb_z_base_sparse_mat class(psb_zspmat_type), intent(inout) :: a integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_z_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_z_csall end interface @@ -1384,12 +1388,14 @@ module psb_z_mat_mod end interface interface - subroutine psb_lz_csall(nr,nc,a,info,nz) - import :: psb_ipk_, psb_lpk_, psb_lzspmat_type + subroutine psb_lz_csall(nr,nc,a,info,nz,type,mold) + import :: psb_ipk_, psb_lpk_, psb_lzspmat_type, psb_lz_base_sparse_mat class(psb_lzspmat_type), intent(inout) :: a integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_lz_base_sparse_mat), optional, intent(in) :: mold end subroutine psb_lz_csall end interface diff --git a/base/serial/impl/psb_c_mat_impl.F90 b/base/serial/impl/psb_c_mat_impl.F90 index 69c67d02..cc112015 100644 --- a/base/serial/impl/psb_c_mat_impl.F90 +++ b/base/serial/impl/psb_c_mat_impl.F90 @@ -582,7 +582,7 @@ end subroutine psb_c_get_neigh -subroutine psb_c_csall(nr,nc,a,info,nz) +subroutine psb_c_csall(nr,nc,a,info,nz,type,mold) use psb_c_mat_mod, psb_protect_name => psb_c_csall use psb_c_base_mat_mod use psb_error_mod @@ -591,6 +591,8 @@ subroutine psb_c_csall(nr,nc,a,info,nz) integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_c_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -601,7 +603,23 @@ subroutine psb_c_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_c_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_c_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_c_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_c_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_c_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_c_coo_sparse_mat :: a%a, stat=info) + end if + if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) @@ -3381,7 +3399,7 @@ end subroutine psb_lc_get_neigh -subroutine psb_lc_csall(nr,nc,a,info,nz) +subroutine psb_lc_csall(nr,nc,a,info,nz,type,mold) use psb_c_mat_mod, psb_protect_name => psb_lc_csall use psb_c_base_mat_mod use psb_error_mod @@ -3390,6 +3408,8 @@ subroutine psb_lc_csall(nr,nc,a,info,nz) integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_lc_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -3400,7 +3420,22 @@ subroutine psb_lc_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_lc_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_lc_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_lc_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_lc_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_lc_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_lc_coo_sparse_mat :: a%a, stat=info) + end if if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) diff --git a/base/serial/impl/psb_d_mat_impl.F90 b/base/serial/impl/psb_d_mat_impl.F90 index 86de5536..7f4ac0c1 100644 --- a/base/serial/impl/psb_d_mat_impl.F90 +++ b/base/serial/impl/psb_d_mat_impl.F90 @@ -582,7 +582,7 @@ end subroutine psb_d_get_neigh -subroutine psb_d_csall(nr,nc,a,info,nz) +subroutine psb_d_csall(nr,nc,a,info,nz,type,mold) use psb_d_mat_mod, psb_protect_name => psb_d_csall use psb_d_base_mat_mod use psb_error_mod @@ -591,6 +591,8 @@ subroutine psb_d_csall(nr,nc,a,info,nz) integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_d_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -601,7 +603,23 @@ subroutine psb_d_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_d_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_d_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_d_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_d_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_d_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_d_coo_sparse_mat :: a%a, stat=info) + end if + if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) @@ -3381,7 +3399,7 @@ end subroutine psb_ld_get_neigh -subroutine psb_ld_csall(nr,nc,a,info,nz) +subroutine psb_ld_csall(nr,nc,a,info,nz,type,mold) use psb_d_mat_mod, psb_protect_name => psb_ld_csall use psb_d_base_mat_mod use psb_error_mod @@ -3390,6 +3408,8 @@ subroutine psb_ld_csall(nr,nc,a,info,nz) integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_ld_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -3400,7 +3420,22 @@ subroutine psb_ld_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_ld_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_ld_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_ld_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_ld_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_ld_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_ld_coo_sparse_mat :: a%a, stat=info) + end if if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) diff --git a/base/serial/impl/psb_s_mat_impl.F90 b/base/serial/impl/psb_s_mat_impl.F90 index 867f9fa4..806a08e3 100644 --- a/base/serial/impl/psb_s_mat_impl.F90 +++ b/base/serial/impl/psb_s_mat_impl.F90 @@ -582,7 +582,7 @@ end subroutine psb_s_get_neigh -subroutine psb_s_csall(nr,nc,a,info,nz) +subroutine psb_s_csall(nr,nc,a,info,nz,type,mold) use psb_s_mat_mod, psb_protect_name => psb_s_csall use psb_s_base_mat_mod use psb_error_mod @@ -591,6 +591,8 @@ subroutine psb_s_csall(nr,nc,a,info,nz) integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_s_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -601,7 +603,23 @@ subroutine psb_s_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_s_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_s_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_s_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_s_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_s_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_s_coo_sparse_mat :: a%a, stat=info) + end if + if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) @@ -3381,7 +3399,7 @@ end subroutine psb_ls_get_neigh -subroutine psb_ls_csall(nr,nc,a,info,nz) +subroutine psb_ls_csall(nr,nc,a,info,nz,type,mold) use psb_s_mat_mod, psb_protect_name => psb_ls_csall use psb_s_base_mat_mod use psb_error_mod @@ -3390,6 +3408,8 @@ subroutine psb_ls_csall(nr,nc,a,info,nz) integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_ls_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -3400,7 +3420,22 @@ subroutine psb_ls_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_ls_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_ls_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_ls_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_ls_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_ls_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_ls_coo_sparse_mat :: a%a, stat=info) + end if if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) diff --git a/base/serial/impl/psb_z_mat_impl.F90 b/base/serial/impl/psb_z_mat_impl.F90 index 07616c05..422a664d 100644 --- a/base/serial/impl/psb_z_mat_impl.F90 +++ b/base/serial/impl/psb_z_mat_impl.F90 @@ -582,7 +582,7 @@ end subroutine psb_z_get_neigh -subroutine psb_z_csall(nr,nc,a,info,nz) +subroutine psb_z_csall(nr,nc,a,info,nz,type,mold) use psb_z_mat_mod, psb_protect_name => psb_z_csall use psb_z_base_mat_mod use psb_error_mod @@ -591,6 +591,8 @@ subroutine psb_z_csall(nr,nc,a,info,nz) integer(psb_ipk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_z_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -601,7 +603,23 @@ subroutine psb_z_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_z_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_z_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_z_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_z_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_z_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_z_coo_sparse_mat :: a%a, stat=info) + end if + if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) @@ -3381,7 +3399,7 @@ end subroutine psb_lz_get_neigh -subroutine psb_lz_csall(nr,nc,a,info,nz) +subroutine psb_lz_csall(nr,nc,a,info,nz,type,mold) use psb_z_mat_mod, psb_protect_name => psb_lz_csall use psb_z_base_mat_mod use psb_error_mod @@ -3390,6 +3408,8 @@ subroutine psb_lz_csall(nr,nc,a,info,nz) integer(psb_lpk_), intent(in) :: nr,nc integer(psb_ipk_), intent(out) :: info integer(psb_lpk_), intent(in), optional :: nz + character(len=*), intent(in), optional :: type + class(psb_lz_base_sparse_mat), optional, intent(in) :: mold integer(psb_ipk_) :: err_act character(len=20) :: name='csall' @@ -3400,7 +3420,22 @@ subroutine psb_lz_csall(nr,nc,a,info,nz) call a%free() info = psb_success_ - allocate(psb_lz_coo_sparse_mat :: a%a, stat=info) + if (present(mold)) then + allocate(a%a, stat=info, mold=mold) + else if (present(type)) then + select case (type) + case('CSR') + allocate(psb_lz_csr_sparse_mat :: a%a, stat=info) + case('COO') + allocate(psb_lz_coo_sparse_mat :: a%a, stat=info) + case('CSC') + allocate(psb_lz_csc_sparse_mat :: a%a, stat=info) + case default + allocate(psb_lz_coo_sparse_mat :: a%a, stat=info) + end select + else + allocate(psb_lz_coo_sparse_mat :: a%a, stat=info) + end if if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_errpush(info, name) diff --git a/prec/impl/psb_c_bjacprec_impl.f90 b/prec/impl/psb_c_bjacprec_impl.f90 index de453684..3cef79dc 100644 --- a/prec/impl/psb_c_bjacprec_impl.f90 +++ b/prec/impl/psb_c_bjacprec_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,13 +27,13 @@ ! 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 psb_c_bjac_dump(prec,info,prefix,head) use psb_base_mod use psb_c_bjacprec, psb_protect_name => psb_c_bjac_dump - implicit none + implicit none class(psb_c_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head @@ -42,13 +42,13 @@ subroutine psb_c_bjac_dump(prec,info,prefix,head) character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - ! len of prefix_ + ! len of prefix_ info = 0 ictxt = prec%get_ctxt() call psb_info(ictxt,iam,np) - if (present(prefix)) then + if (present(prefix)) then prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) else prefix_ = "dump_fact_c" @@ -73,7 +73,7 @@ end subroutine psb_c_bjac_dump subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_c_bjacprec, psb_protect_name => psb_c_bjac_apply_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_c_bjac_prec_type), intent(inout) :: prec complex(psb_spk_),intent(in) :: alpha,beta @@ -116,12 +116,12 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then + if (x%get_nrows() < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (y%get_nrows() < n_row) then + if (y%get_nrows() < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -138,9 +138,9 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then aux => work(n_col+1:) else allocate(aux(4*n_col),stat=info) @@ -150,19 +150,19 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) allocate(ww(n_col),aux(4*n_col),stat=info) endif - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if do_alloc_wrk = .not.prec%is_allocated_wrk() if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) - + select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) + select case(trans_) case('N') call psb_spsm(cone,prec%av(psb_l_pr_),x,czero,wv,desc_data,info,& @@ -170,31 +170,31 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + case('T') call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + case('C') - + call psb_spsm(cone,prec%av(psb_u_pr_),x,czero,wv,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + call wv1%mlt(cone,prec%dv,wv,czero,info,conjgx=trans_) - + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if - + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -202,12 +202,12 @@ subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end select end associate - + call psb_halo(y,desc_data,info,data=psb_comm_mov_) if (do_alloc_wrk) call prec%free_wrk(info) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -229,7 +229,7 @@ end subroutine psb_c_bjac_apply_vect subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_c_bjacprec, psb_protect_name => psb_c_bjac_apply - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_c_bjac_prec_type), intent(inout) :: prec complex(psb_spk_),intent(in) :: alpha,beta @@ -270,12 +270,12 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (size(x) < n_row) then + if (size(x) < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (size(y) < n_row) then + if (size(y) < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -292,29 +292,29 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + 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 + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif else allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) select case(trans_) case('N') @@ -341,7 +341,7 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) & trans=trans_,scale='U',choice=psb_none_,work=aux) end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if @@ -355,8 +355,8 @@ subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) call psb_halo(y,desc_data,info,data=psb_comm_mov_) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -389,6 +389,7 @@ subroutine psb_c_bjac_precinit(prec,info) info = psb_success_ call psb_realloc(psb_ifpsz,prec%iprcparm,info) + call psb_realloc(psb_rfpsz,prec%rprcparm,info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_Errpush(info,name) @@ -399,6 +400,11 @@ subroutine psb_c_bjac_precinit(prec,info) prec%iprcparm(psb_p_type_) = psb_bjac_ prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ prec%iprcparm(psb_ilu_fill_in_) = 0 + prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ + prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + + prec%rprcparm(:) = 0 + prec%rprcparm(psb_fact_eps_) = 1E-1_psb_spk_ call psb_erractionrestore(err_act) @@ -413,7 +419,7 @@ end subroutine psb_c_bjac_precinit subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - use psb_prec_mod, only : psb_ilu_fct + use psb_c_ilu_fact_mod use psb_c_bjacprec, psb_protect_name => psb_c_bjac_precbld Implicit None @@ -425,12 +431,13 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! .. Local Scalars .. - integer(psb_ipk_) :: i, m + ! .. Local Scalars .. + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale integer(psb_ipk_) :: ierr(5) character :: trans, unitd - type(psb_c_csr_sparse_mat), allocatable :: lf, uf + type(psb_cspmat_type), allocatable :: lf, uf complex(psb_spk_), allocatable :: dd(:) + real(psb_spk_) :: fact_eps integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='c_bjac_precbld' @@ -458,19 +465,214 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) trans = 'N' unitd = 'U' + ! We check if all the information contained in the preconditioner structure + ! are meaningful, otherwise we give an error and get out of the build + ! procedure + ialg = prec%iprcparm(psb_ilu_ialg_) + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_ialg_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + iscale = prec%iprcparm(psb_ilu_scale_) + if ((iscale == psb_ilu_scale_none_).or.& + (iscale == psb_ilu_scale_maxval_).or.& + (iscale == psb_ilu_scale_diag_).or.& + (iscale == psb_ilu_scale_arwsum_).or.& + (iscale == psb_ilu_scale_aclsum_).or.& + (iscale == psb_ilu_scale_arcsum_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_scale_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fact_eps = prec%rprcparm(psb_fact_eps_) + if(fact_eps > 1 ) then + info=psb_err_from_subroutine_ + ch_err='psb_fact_eps_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + + ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_) + ! ILU(0) Factorization: the total number of nonzeros of the factorized matrix + ! is equal to the one of the input matrix + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + ! call psb_ilu_fct(a,lf,uf,dd,info) + call psb_ilu0_fact(ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu0_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_k_) + ! ILU(N) Incomplete LU-factorization with N levels of fill-in. Depending on + ! the type of the variant of the algorithm the may be forgotten or added to + ! the diagonal (MILU) + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if - if (allocated(prec%av)) then - if (size(prec%av) < psb_bp_ilu_avsz) then - do i=1,size(prec%av) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + call psb_iluk_fact(fill_in,ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_iluk_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_t_) + ! ILU(N,E) Incomplete LU factorization with thresholding and level of fill + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) endif end if - if (.not.allocated(prec%av)) then + if (.not.allocated(prec%av)) then allocate(prec%av(psb_max_avsz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_dealloc_,name) @@ -497,27 +699,27 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) end if allocate(dd(n_row),stat=info) - if (info == psb_success_) then + if (info == psb_success_) then allocate(prec%dv, stat=info) - if (info == 0) then - if (present(vmold)) then - allocate(prec%dv%v,mold=vmold,stat=info) - else - allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_c_base_vect_type :: prec%dv%v,stat=info) end if end if end if - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 endif - ! This is where we have no renumbering, thus no need - call psb_ilu_fct(a,lf,uf,dd,info) + ! This is where we have no renumbering, thus no need + call psb_ilut_fact(fill_in,fact_eps,a,lf,uf,dd,info,iscale=iscale) if(info == psb_success_) then - call prec%av(psb_l_pr_)%mv_from(lf) - call prec%av(psb_u_pr_)%mv_from(uf) + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) call prec%av(psb_l_pr_)%set_asb() call prec%av(psb_u_pr_)%set_asb() call prec%av(psb_l_pr_)%trim() @@ -526,12 +728,12 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! call move_alloc(dd,prec%d) else info=psb_err_from_subroutine_ - ch_err='psb_ilu_fct' + ch_err='psb_ilut_fact' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - case(psb_f_none_) + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' call psb_errpush(info,name,a_err=ch_err) @@ -544,7 +746,7 @@ subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end select - if (present(amold)) then + if (present(amold)) then call prec%av(psb_l_pr_)%cscnv(info,mold=amold) call prec%av(psb_u_pr_)%cscnv(info,mold=amold) end if @@ -563,8 +765,8 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) Implicit None class(psb_c_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, nrow character(len=20) :: name='c_bjac_precset' @@ -572,15 +774,15 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) call psb_erractionsave(err_act) info = psb_success_ - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if select case(what) - case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -588,9 +790,10 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_f_type_) = val - case (psb_ilu_fill_in_) + case (psb_ilu_fill_in_) if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -598,6 +801,24 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_ilu_fill_in_) = val + case (psb_ilu_ialg_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_ialg_) = val + + case (psb_ilu_scale_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_scale_) = val + case default write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' @@ -609,3 +830,57 @@ subroutine psb_c_bjac_precseti(prec,what,val,info) 9999 call psb_error_handler(err_act) return end subroutine psb_c_bjac_precseti + +subroutine psb_c_bjac_precsetr(prec,what,val,info) + + use psb_base_mod + use psb_c_bjacprec, psb_protect_name => psb_c_bjac_precsetr + Implicit None + + class(psb_c_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, nrow + character(len=20) :: name='c_bjac_precset' + + call psb_erractionsave(err_act) + + info = psb_success_ + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + select case(what) + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_f_type_) = val + + case (psb_fact_eps_) + if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%rprcparm(psb_fact_eps_) = val + + case default + write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_c_bjac_precsetr diff --git a/prec/impl/psb_d_bjacprec_impl.f90 b/prec/impl/psb_d_bjacprec_impl.f90 index ef8c52c3..34a2b63a 100644 --- a/prec/impl/psb_d_bjacprec_impl.f90 +++ b/prec/impl/psb_d_bjacprec_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,13 +27,13 @@ ! 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 psb_d_bjac_dump(prec,info,prefix,head) use psb_base_mod use psb_d_bjacprec, psb_protect_name => psb_d_bjac_dump - implicit none + implicit none class(psb_d_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head @@ -42,13 +42,13 @@ subroutine psb_d_bjac_dump(prec,info,prefix,head) character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - ! len of prefix_ + ! len of prefix_ info = 0 ictxt = prec%get_ctxt() call psb_info(ictxt,iam,np) - if (present(prefix)) then + if (present(prefix)) then prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) else prefix_ = "dump_fact_d" @@ -73,7 +73,7 @@ end subroutine psb_d_bjac_dump subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_d_bjacprec, psb_protect_name => psb_d_bjac_apply_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_d_bjac_prec_type), intent(inout) :: prec real(psb_dpk_),intent(in) :: alpha,beta @@ -116,12 +116,12 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then + if (x%get_nrows() < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (y%get_nrows() < n_row) then + if (y%get_nrows() < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -138,9 +138,9 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then aux => work(n_col+1:) else allocate(aux(4*n_col),stat=info) @@ -150,19 +150,19 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) allocate(ww(n_col),aux(4*n_col),stat=info) endif - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if do_alloc_wrk = .not.prec%is_allocated_wrk() if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) - + select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) + select case(trans_) case('N') call psb_spsm(done,prec%av(psb_l_pr_),x,dzero,wv,desc_data,info,& @@ -170,31 +170,31 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + case('T') call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + case('C') - + call psb_spsm(done,prec%av(psb_u_pr_),x,dzero,wv,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + call wv1%mlt(done,prec%dv,wv,dzero,info,conjgx=trans_) - + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if - + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -202,12 +202,12 @@ subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end select end associate - + call psb_halo(y,desc_data,info,data=psb_comm_mov_) if (do_alloc_wrk) call prec%free_wrk(info) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -229,7 +229,7 @@ end subroutine psb_d_bjac_apply_vect subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_d_bjacprec, psb_protect_name => psb_d_bjac_apply - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_d_bjac_prec_type), intent(inout) :: prec real(psb_dpk_),intent(in) :: alpha,beta @@ -270,12 +270,12 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (size(x) < n_row) then + if (size(x) < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (size(y) < n_row) then + if (size(y) < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -292,29 +292,29 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + 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 + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif else allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) select case(trans_) case('N') @@ -341,7 +341,7 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) & trans=trans_,scale='U',choice=psb_none_,work=aux) end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if @@ -355,8 +355,8 @@ subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) call psb_halo(y,desc_data,info,data=psb_comm_mov_) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -389,6 +389,7 @@ subroutine psb_d_bjac_precinit(prec,info) info = psb_success_ call psb_realloc(psb_ifpsz,prec%iprcparm,info) + call psb_realloc(psb_rfpsz,prec%rprcparm,info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_Errpush(info,name) @@ -399,6 +400,11 @@ subroutine psb_d_bjac_precinit(prec,info) prec%iprcparm(psb_p_type_) = psb_bjac_ prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ prec%iprcparm(psb_ilu_fill_in_) = 0 + prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ + prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + + prec%rprcparm(:) = 0 + prec%rprcparm(psb_fact_eps_) = 1E-1_psb_dpk_ call psb_erractionrestore(err_act) @@ -413,7 +419,7 @@ end subroutine psb_d_bjac_precinit subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - use psb_prec_mod, only : psb_ilu_fct + use psb_d_ilu_fact_mod use psb_d_bjacprec, psb_protect_name => psb_d_bjac_precbld Implicit None @@ -425,12 +431,13 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! .. Local Scalars .. - integer(psb_ipk_) :: i, m + ! .. Local Scalars .. + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale integer(psb_ipk_) :: ierr(5) character :: trans, unitd - type(psb_d_csr_sparse_mat), allocatable :: lf, uf + type(psb_dspmat_type), allocatable :: lf, uf real(psb_dpk_), allocatable :: dd(:) + real(psb_dpk_) :: fact_eps integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='d_bjac_precbld' @@ -458,19 +465,214 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) trans = 'N' unitd = 'U' + ! We check if all the information contained in the preconditioner structure + ! are meaningful, otherwise we give an error and get out of the build + ! procedure + ialg = prec%iprcparm(psb_ilu_ialg_) + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_ialg_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + iscale = prec%iprcparm(psb_ilu_scale_) + if ((iscale == psb_ilu_scale_none_).or.& + (iscale == psb_ilu_scale_maxval_).or.& + (iscale == psb_ilu_scale_diag_).or.& + (iscale == psb_ilu_scale_arwsum_).or.& + (iscale == psb_ilu_scale_aclsum_).or.& + (iscale == psb_ilu_scale_arcsum_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_scale_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fact_eps = prec%rprcparm(psb_fact_eps_) + if(fact_eps > 1 ) then + info=psb_err_from_subroutine_ + ch_err='psb_fact_eps_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + + ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_) + ! ILU(0) Factorization: the total number of nonzeros of the factorized matrix + ! is equal to the one of the input matrix + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + ! call psb_ilu_fct(a,lf,uf,dd,info) + call psb_ilu0_fact(ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu0_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_k_) + ! ILU(N) Incomplete LU-factorization with N levels of fill-in. Depending on + ! the type of the variant of the algorithm the may be forgotten or added to + ! the diagonal (MILU) + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if - if (allocated(prec%av)) then - if (size(prec%av) < psb_bp_ilu_avsz) then - do i=1,size(prec%av) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + call psb_iluk_fact(fill_in,ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_iluk_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_t_) + ! ILU(N,E) Incomplete LU factorization with thresholding and level of fill + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) endif end if - if (.not.allocated(prec%av)) then + if (.not.allocated(prec%av)) then allocate(prec%av(psb_max_avsz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_dealloc_,name) @@ -497,27 +699,27 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) end if allocate(dd(n_row),stat=info) - if (info == psb_success_) then + if (info == psb_success_) then allocate(prec%dv, stat=info) - if (info == 0) then - if (present(vmold)) then - allocate(prec%dv%v,mold=vmold,stat=info) - else - allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_d_base_vect_type :: prec%dv%v,stat=info) end if end if end if - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 endif - ! This is where we have no renumbering, thus no need - call psb_ilu_fct(a,lf,uf,dd,info) + ! This is where we have no renumbering, thus no need + call psb_ilut_fact(fill_in,fact_eps,a,lf,uf,dd,info,iscale=iscale) if(info == psb_success_) then - call prec%av(psb_l_pr_)%mv_from(lf) - call prec%av(psb_u_pr_)%mv_from(uf) + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) call prec%av(psb_l_pr_)%set_asb() call prec%av(psb_u_pr_)%set_asb() call prec%av(psb_l_pr_)%trim() @@ -526,12 +728,12 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! call move_alloc(dd,prec%d) else info=psb_err_from_subroutine_ - ch_err='psb_ilu_fct' + ch_err='psb_ilut_fact' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - case(psb_f_none_) + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' call psb_errpush(info,name,a_err=ch_err) @@ -544,7 +746,7 @@ subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end select - if (present(amold)) then + if (present(amold)) then call prec%av(psb_l_pr_)%cscnv(info,mold=amold) call prec%av(psb_u_pr_)%cscnv(info,mold=amold) end if @@ -563,8 +765,8 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) Implicit None class(psb_d_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, nrow character(len=20) :: name='d_bjac_precset' @@ -572,15 +774,15 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) call psb_erractionsave(err_act) info = psb_success_ - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if select case(what) - case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -588,9 +790,10 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_f_type_) = val - case (psb_ilu_fill_in_) + case (psb_ilu_fill_in_) if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -598,6 +801,24 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_ilu_fill_in_) = val + case (psb_ilu_ialg_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_ialg_) = val + + case (psb_ilu_scale_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_scale_) = val + case default write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' @@ -609,3 +830,57 @@ subroutine psb_d_bjac_precseti(prec,what,val,info) 9999 call psb_error_handler(err_act) return end subroutine psb_d_bjac_precseti + +subroutine psb_d_bjac_precsetr(prec,what,val,info) + + use psb_base_mod + use psb_d_bjacprec, psb_protect_name => psb_d_bjac_precsetr + Implicit None + + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, nrow + character(len=20) :: name='d_bjac_precset' + + call psb_erractionsave(err_act) + + info = psb_success_ + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + select case(what) + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_f_type_) = val + + case (psb_fact_eps_) + if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%rprcparm(psb_fact_eps_) = val + + case default + write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_d_bjac_precsetr diff --git a/prec/impl/psb_s_bjacprec_impl.f90 b/prec/impl/psb_s_bjacprec_impl.f90 index 3a9cfce2..1448e43f 100644 --- a/prec/impl/psb_s_bjacprec_impl.f90 +++ b/prec/impl/psb_s_bjacprec_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,13 +27,13 @@ ! 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 psb_s_bjac_dump(prec,info,prefix,head) use psb_base_mod use psb_s_bjacprec, psb_protect_name => psb_s_bjac_dump - implicit none + implicit none class(psb_s_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head @@ -42,13 +42,13 @@ subroutine psb_s_bjac_dump(prec,info,prefix,head) character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - ! len of prefix_ + ! len of prefix_ info = 0 ictxt = prec%get_ctxt() call psb_info(ictxt,iam,np) - if (present(prefix)) then + if (present(prefix)) then prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) else prefix_ = "dump_fact_s" @@ -73,7 +73,7 @@ end subroutine psb_s_bjac_dump subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_s_bjacprec, psb_protect_name => psb_s_bjac_apply_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_s_bjac_prec_type), intent(inout) :: prec real(psb_spk_),intent(in) :: alpha,beta @@ -116,12 +116,12 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then + if (x%get_nrows() < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (y%get_nrows() < n_row) then + if (y%get_nrows() < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -138,9 +138,9 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then aux => work(n_col+1:) else allocate(aux(4*n_col),stat=info) @@ -150,19 +150,19 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) allocate(ww(n_col),aux(4*n_col),stat=info) endif - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if do_alloc_wrk = .not.prec%is_allocated_wrk() if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) - + select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) + select case(trans_) case('N') call psb_spsm(sone,prec%av(psb_l_pr_),x,szero,wv,desc_data,info,& @@ -170,31 +170,31 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + case('T') call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + case('C') - + call psb_spsm(sone,prec%av(psb_u_pr_),x,szero,wv,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + call wv1%mlt(sone,prec%dv,wv,szero,info,conjgx=trans_) - + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if - + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -202,12 +202,12 @@ subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end select end associate - + call psb_halo(y,desc_data,info,data=psb_comm_mov_) if (do_alloc_wrk) call prec%free_wrk(info) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -229,7 +229,7 @@ end subroutine psb_s_bjac_apply_vect subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_s_bjacprec, psb_protect_name => psb_s_bjac_apply - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_s_bjac_prec_type), intent(inout) :: prec real(psb_spk_),intent(in) :: alpha,beta @@ -270,12 +270,12 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (size(x) < n_row) then + if (size(x) < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (size(y) < n_row) then + if (size(y) < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -292,29 +292,29 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + 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 + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif else allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) select case(trans_) case('N') @@ -341,7 +341,7 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) & trans=trans_,scale='U',choice=psb_none_,work=aux) end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if @@ -355,8 +355,8 @@ subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) call psb_halo(y,desc_data,info,data=psb_comm_mov_) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -389,6 +389,7 @@ subroutine psb_s_bjac_precinit(prec,info) info = psb_success_ call psb_realloc(psb_ifpsz,prec%iprcparm,info) + call psb_realloc(psb_rfpsz,prec%rprcparm,info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_Errpush(info,name) @@ -399,6 +400,11 @@ subroutine psb_s_bjac_precinit(prec,info) prec%iprcparm(psb_p_type_) = psb_bjac_ prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ prec%iprcparm(psb_ilu_fill_in_) = 0 + prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ + prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + + prec%rprcparm(:) = 0 + prec%rprcparm(psb_fact_eps_) = 1E-1_psb_spk_ call psb_erractionrestore(err_act) @@ -413,7 +419,7 @@ end subroutine psb_s_bjac_precinit subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - use psb_prec_mod, only : psb_ilu_fct + use psb_s_ilu_fact_mod use psb_s_bjacprec, psb_protect_name => psb_s_bjac_precbld Implicit None @@ -425,12 +431,13 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! .. Local Scalars .. - integer(psb_ipk_) :: i, m + ! .. Local Scalars .. + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale integer(psb_ipk_) :: ierr(5) character :: trans, unitd - type(psb_s_csr_sparse_mat), allocatable :: lf, uf + type(psb_sspmat_type), allocatable :: lf, uf real(psb_spk_), allocatable :: dd(:) + real(psb_spk_) :: fact_eps integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='s_bjac_precbld' @@ -458,19 +465,214 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) trans = 'N' unitd = 'U' + ! We check if all the information contained in the preconditioner structure + ! are meaningful, otherwise we give an error and get out of the build + ! procedure + ialg = prec%iprcparm(psb_ilu_ialg_) + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_ialg_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + iscale = prec%iprcparm(psb_ilu_scale_) + if ((iscale == psb_ilu_scale_none_).or.& + (iscale == psb_ilu_scale_maxval_).or.& + (iscale == psb_ilu_scale_diag_).or.& + (iscale == psb_ilu_scale_arwsum_).or.& + (iscale == psb_ilu_scale_aclsum_).or.& + (iscale == psb_ilu_scale_arcsum_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_scale_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fact_eps = prec%rprcparm(psb_fact_eps_) + if(fact_eps > 1 ) then + info=psb_err_from_subroutine_ + ch_err='psb_fact_eps_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + + ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_) + ! ILU(0) Factorization: the total number of nonzeros of the factorized matrix + ! is equal to the one of the input matrix + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + ! call psb_ilu_fct(a,lf,uf,dd,info) + call psb_ilu0_fact(ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu0_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_k_) + ! ILU(N) Incomplete LU-factorization with N levels of fill-in. Depending on + ! the type of the variant of the algorithm the may be forgotten or added to + ! the diagonal (MILU) + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if - if (allocated(prec%av)) then - if (size(prec%av) < psb_bp_ilu_avsz) then - do i=1,size(prec%av) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + call psb_iluk_fact(fill_in,ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_iluk_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_t_) + ! ILU(N,E) Incomplete LU factorization with thresholding and level of fill + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) endif end if - if (.not.allocated(prec%av)) then + if (.not.allocated(prec%av)) then allocate(prec%av(psb_max_avsz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_dealloc_,name) @@ -497,27 +699,27 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) end if allocate(dd(n_row),stat=info) - if (info == psb_success_) then + if (info == psb_success_) then allocate(prec%dv, stat=info) - if (info == 0) then - if (present(vmold)) then - allocate(prec%dv%v,mold=vmold,stat=info) - else - allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_s_base_vect_type :: prec%dv%v,stat=info) end if end if end if - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 endif - ! This is where we have no renumbering, thus no need - call psb_ilu_fct(a,lf,uf,dd,info) + ! This is where we have no renumbering, thus no need + call psb_ilut_fact(fill_in,fact_eps,a,lf,uf,dd,info,iscale=iscale) if(info == psb_success_) then - call prec%av(psb_l_pr_)%mv_from(lf) - call prec%av(psb_u_pr_)%mv_from(uf) + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) call prec%av(psb_l_pr_)%set_asb() call prec%av(psb_u_pr_)%set_asb() call prec%av(psb_l_pr_)%trim() @@ -526,12 +728,12 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! call move_alloc(dd,prec%d) else info=psb_err_from_subroutine_ - ch_err='psb_ilu_fct' + ch_err='psb_ilut_fact' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - case(psb_f_none_) + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' call psb_errpush(info,name,a_err=ch_err) @@ -544,7 +746,7 @@ subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end select - if (present(amold)) then + if (present(amold)) then call prec%av(psb_l_pr_)%cscnv(info,mold=amold) call prec%av(psb_u_pr_)%cscnv(info,mold=amold) end if @@ -563,8 +765,8 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) Implicit None class(psb_s_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, nrow character(len=20) :: name='s_bjac_precset' @@ -572,15 +774,15 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) call psb_erractionsave(err_act) info = psb_success_ - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if select case(what) - case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -588,9 +790,10 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_f_type_) = val - case (psb_ilu_fill_in_) + case (psb_ilu_fill_in_) if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -598,6 +801,24 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_ilu_fill_in_) = val + case (psb_ilu_ialg_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_ialg_) = val + + case (psb_ilu_scale_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_scale_) = val + case default write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' @@ -609,3 +830,57 @@ subroutine psb_s_bjac_precseti(prec,what,val,info) 9999 call psb_error_handler(err_act) return end subroutine psb_s_bjac_precseti + +subroutine psb_s_bjac_precsetr(prec,what,val,info) + + use psb_base_mod + use psb_s_bjacprec, psb_protect_name => psb_s_bjac_precsetr + Implicit None + + class(psb_s_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, nrow + character(len=20) :: name='s_bjac_precset' + + call psb_erractionsave(err_act) + + info = psb_success_ + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + select case(what) + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_f_type_) = val + + case (psb_fact_eps_) + if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%rprcparm(psb_fact_eps_) = val + + case default + write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_s_bjac_precsetr diff --git a/prec/impl/psb_z_bjacprec_impl.f90 b/prec/impl/psb_z_bjacprec_impl.f90 index b70018f4..a1859e87 100644 --- a/prec/impl/psb_z_bjacprec_impl.f90 +++ b/prec/impl/psb_z_bjacprec_impl.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,13 +27,13 @@ ! 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 psb_z_bjac_dump(prec,info,prefix,head) use psb_base_mod use psb_z_bjacprec, psb_protect_name => psb_z_bjac_dump - implicit none + implicit none class(psb_z_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(out) :: info character(len=*), intent(in), optional :: prefix,head @@ -42,13 +42,13 @@ subroutine psb_z_bjac_dump(prec,info,prefix,head) character(len=80) :: prefix_ character(len=120) :: fname ! len should be at least 20 more than - ! len of prefix_ + ! len of prefix_ info = 0 ictxt = prec%get_ctxt() call psb_info(ictxt,iam,np) - if (present(prefix)) then + if (present(prefix)) then prefix_ = trim(prefix(1:min(len(prefix),len(prefix_)))) else prefix_ = "dump_fact_z" @@ -73,7 +73,7 @@ end subroutine psb_z_bjac_dump subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_z_bjacprec, psb_protect_name => psb_z_bjac_apply_vect - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_z_bjac_prec_type), intent(inout) :: prec complex(psb_dpk_),intent(in) :: alpha,beta @@ -116,12 +116,12 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (x%get_nrows() < n_row) then + if (x%get_nrows() < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (y%get_nrows() < n_row) then + if (y%get_nrows() < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -138,9 +138,9 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then aux => work(n_col+1:) else allocate(aux(4*n_col),stat=info) @@ -150,19 +150,19 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) allocate(ww(n_col),aux(4*n_col),stat=info) endif - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if do_alloc_wrk = .not.prec%is_allocated_wrk() if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v) associate (wv => prec%wrk(1), wv1 => prec%wrk(2)) - + select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) - + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) + select case(trans_) case('N') call psb_spsm(zone,prec%av(psb_l_pr_),x,zzero,wv,desc_data,info,& @@ -170,31 +170,31 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_u_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + case('T') call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& & trans=trans_,scale='L',diag=prec%dv,choice=psb_none_, work=aux) if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + case('C') - + call psb_spsm(zone,prec%av(psb_u_pr_),x,zzero,wv,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_, work=aux) - + call wv1%mlt(zone,prec%dv,wv,zzero,info,conjgx=trans_) - + if(info == psb_success_) call psb_spsm(alpha,prec%av(psb_l_pr_),wv1,& & beta,y,desc_data,info,& & trans=trans_,scale='U',choice=psb_none_,work=aux) - + end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if - + case default info = psb_err_internal_error_ call psb_errpush(info,name,a_err='Invalid factorization') @@ -202,12 +202,12 @@ subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) end select end associate - + call psb_halo(y,desc_data,info,data=psb_comm_mov_) if (do_alloc_wrk) call prec%free_wrk(info) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -229,7 +229,7 @@ end subroutine psb_z_bjac_apply_vect subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) use psb_base_mod use psb_z_bjacprec, psb_protect_name => psb_z_bjac_apply - implicit none + implicit none type(psb_desc_type),intent(in) :: desc_data class(psb_z_bjac_prec_type), intent(inout) :: prec complex(psb_dpk_),intent(in) :: alpha,beta @@ -270,12 +270,12 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) n_row = desc_data%get_local_rows() n_col = desc_data%get_local_cols() - if (size(x) < n_row) then + if (size(x) < n_row) then info = 36; ierr(1) = 2; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 end if - if (size(y) < n_row) then + if (size(y) < n_row) then info = 36; ierr(1) = 3; ierr(2) = n_row; call psb_errpush(info,name,i_err=ierr) goto 9999 @@ -292,29 +292,29 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) end if - if (n_col <= size(work)) then + if (n_col <= size(work)) then ww => work(1:n_col) - if ((4*n_col+n_col) <= size(work)) then + 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 + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif else allocate(ww(n_col),aux(4*n_col),stat=info) - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 end if endif select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_,psb_f_ilu_k_,psb_f_ilu_t_) select case(trans_) case('N') @@ -341,7 +341,7 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) & trans=trans_,scale='U',choice=psb_none_,work=aux) end select - if (info /= psb_success_) then + if (info /= psb_success_) then ch_err="psb_spsm" goto 9999 end if @@ -355,8 +355,8 @@ subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) call psb_halo(y,desc_data,info,data=psb_comm_mov_) - if (n_col <= size(work)) then - if ((4*n_col+n_col) <= size(work)) then + if (n_col <= size(work)) then + if ((4*n_col+n_col) <= size(work)) then else deallocate(aux) endif @@ -389,6 +389,7 @@ subroutine psb_z_bjac_precinit(prec,info) info = psb_success_ call psb_realloc(psb_ifpsz,prec%iprcparm,info) + call psb_realloc(psb_rfpsz,prec%rprcparm,info) if (info /= psb_success_) then info = psb_err_alloc_dealloc_ call psb_Errpush(info,name) @@ -399,6 +400,11 @@ subroutine psb_z_bjac_precinit(prec,info) prec%iprcparm(psb_p_type_) = psb_bjac_ prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ prec%iprcparm(psb_ilu_fill_in_) = 0 + prec%iprcparm(psb_ilu_ialg_) = psb_ilu_n_ + prec%iprcparm(psb_ilu_scale_) = psb_ilu_scale_none_ + + prec%rprcparm(:) = 0 + prec%rprcparm(psb_fact_eps_) = 1E-1_psb_dpk_ call psb_erractionrestore(err_act) @@ -413,7 +419,7 @@ end subroutine psb_z_bjac_precinit subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) use psb_base_mod - use psb_prec_mod, only : psb_ilu_fct + use psb_z_ilu_fact_mod use psb_z_bjacprec, psb_protect_name => psb_z_bjac_precbld Implicit None @@ -425,12 +431,13 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_i_base_vect_type), intent(in), optional :: imold - ! .. Local Scalars .. - integer(psb_ipk_) :: i, m + ! .. Local Scalars .. + integer(psb_ipk_) :: i, m, ialg, fill_in, iscale integer(psb_ipk_) :: ierr(5) character :: trans, unitd - type(psb_z_csr_sparse_mat), allocatable :: lf, uf + type(psb_zspmat_type), allocatable :: lf, uf complex(psb_dpk_), allocatable :: dd(:) + real(psb_dpk_) :: fact_eps integer(psb_ipk_) :: nztota, err_act, n_row, nrow_a,n_col, nhalo integer(psb_ipk_) :: ictxt,np,me character(len=20) :: name='z_bjac_precbld' @@ -458,19 +465,214 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) trans = 'N' unitd = 'U' + ! We check if all the information contained in the preconditioner structure + ! are meaningful, otherwise we give an error and get out of the build + ! procedure + ialg = prec%iprcparm(psb_ilu_ialg_) + if ((ialg == psb_ilu_n_).or.(ialg == psb_milu_n_).or.(ialg == psb_ilu_t_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_ialg_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + iscale = prec%iprcparm(psb_ilu_scale_) + if ((iscale == psb_ilu_scale_none_).or.& + (iscale == psb_ilu_scale_maxval_).or.& + (iscale == psb_ilu_scale_diag_).or.& + (iscale == psb_ilu_scale_arwsum_).or.& + (iscale == psb_ilu_scale_aclsum_).or.& + (iscale == psb_ilu_scale_arcsum_)) then + ! Do nothing: admissible request + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu_scale_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fact_eps = prec%rprcparm(psb_fact_eps_) + if(fact_eps > 1 ) then + info=psb_err_from_subroutine_ + ch_err='psb_fact_eps_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + fill_in = prec%iprcparm(psb_ilu_fill_in_) + if(fill_in < 0) then + info=psb_err_from_subroutine_ + ch_err='psb_ilu_fill_in_' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + else if (fill_in == 0) then + ! If the requested level of fill is equal to zero, we default to the + ! specialized ILU(0) routine + prec%iprcparm(psb_f_type_) = psb_f_ilu_n_ + end if + + ! Select on the type of factorization to be used select case(prec%iprcparm(psb_f_type_)) - case(psb_f_ilu_n_) + case(psb_f_ilu_n_) + ! ILU(0) Factorization: the total number of nonzeros of the factorized matrix + ! is equal to the one of the input matrix + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + ! call psb_ilu_fct(a,lf,uf,dd,info) + call psb_ilu0_fact(ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_ilu0_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_k_) + ! ILU(N) Incomplete LU-factorization with N levels of fill-in. Depending on + ! the type of the variant of the algorithm the may be forgotten or added to + ! the diagonal (MILU) + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) + call prec%av(i)%free() + enddo + deallocate(prec%av,stat=info) + endif + end if + if (.not.allocated(prec%av)) then + allocate(prec%av(psb_max_avsz),stat=info) + if (info /= psb_success_) then + call psb_errpush(psb_err_alloc_dealloc_,name) + goto 9999 + end if + endif + + nrow_a = desc_a%get_local_rows() + nztota = a%get_nzeros() + + n_col = desc_a%get_local_cols() + nhalo = n_col-nrow_a + n_row = nrow_a + + allocate(lf,uf,stat=info) + if (info == psb_success_) call lf%allocate(n_row,n_row,nztota) + if (info == psb_success_) call uf%allocate(n_row,n_row,nztota) + + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_sp_all' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(dd(n_row),stat=info) + if (info == psb_success_) then + allocate(prec%dv, stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + end if + end if + end if - if (allocated(prec%av)) then - if (size(prec%av) < psb_bp_ilu_avsz) then - do i=1,size(prec%av) + if (info /= psb_success_) then + call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') + goto 9999 + endif + ! This is where we have no renumbering, thus no need + call psb_iluk_fact(fill_in,ialg,a,lf,uf,dd,info) + + if(info == psb_success_) then + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) + call prec%av(psb_l_pr_)%set_asb() + call prec%av(psb_u_pr_)%set_asb() + call prec%av(psb_l_pr_)%trim() + call prec%av(psb_u_pr_)%trim() + call prec%dv%bld(dd) + ! call move_alloc(dd,prec%d) + else + info=psb_err_from_subroutine_ + ch_err='psb_iluk_fact' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + case(psb_f_ilu_t_) + ! ILU(N,E) Incomplete LU factorization with thresholding and level of fill + + if (allocated(prec%av)) then + if (size(prec%av) < psb_bp_ilu_avsz) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) endif end if - if (.not.allocated(prec%av)) then + if (.not.allocated(prec%av)) then allocate(prec%av(psb_max_avsz),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_alloc_dealloc_,name) @@ -497,27 +699,27 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) end if allocate(dd(n_row),stat=info) - if (info == psb_success_) then + if (info == psb_success_) then allocate(prec%dv, stat=info) - if (info == 0) then - if (present(vmold)) then - allocate(prec%dv%v,mold=vmold,stat=info) - else - allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) + if (info == 0) then + if (present(vmold)) then + allocate(prec%dv%v,mold=vmold,stat=info) + else + allocate(psb_z_base_vect_type :: prec%dv%v,stat=info) end if end if end if - if (info /= psb_success_) then + if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') - goto 9999 + goto 9999 endif - ! This is where we have no renumbering, thus no need - call psb_ilu_fct(a,lf,uf,dd,info) + ! This is where we have no renumbering, thus no need + call psb_ilut_fact(fill_in,fact_eps,a,lf,uf,dd,info,iscale=iscale) if(info == psb_success_) then - call prec%av(psb_l_pr_)%mv_from(lf) - call prec%av(psb_u_pr_)%mv_from(uf) + call prec%av(psb_l_pr_)%mv_from(lf%a) + call prec%av(psb_u_pr_)%mv_from(uf%a) call prec%av(psb_l_pr_)%set_asb() call prec%av(psb_u_pr_)%set_asb() call prec%av(psb_l_pr_)%trim() @@ -526,12 +728,12 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) ! call move_alloc(dd,prec%d) else info=psb_err_from_subroutine_ - ch_err='psb_ilu_fct' + ch_err='psb_ilut_fact' call psb_errpush(info,name,a_err=ch_err) goto 9999 end if - case(psb_f_none_) + case(psb_f_none_) info=psb_err_from_subroutine_ ch_err='Inconsistent prec psb_f_none_' call psb_errpush(info,name,a_err=ch_err) @@ -544,7 +746,7 @@ subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) goto 9999 end select - if (present(amold)) then + if (present(amold)) then call prec%av(psb_l_pr_)%cscnv(info,mold=amold) call prec%av(psb_u_pr_)%cscnv(info,mold=amold) end if @@ -563,8 +765,8 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) Implicit None class(psb_z_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: err_act, nrow character(len=20) :: name='z_bjac_precset' @@ -572,15 +774,15 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) call psb_erractionsave(err_act) info = psb_success_ - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 end if select case(what) - case (psb_f_type_) - if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -588,9 +790,10 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_f_type_) = val - case (psb_ilu_fill_in_) + case (psb_ilu_fill_in_) if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& - & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + & ((prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_t_))) then write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& & prec%iprcparm(psb_p_type_),& & 'ignoring user specification' @@ -598,6 +801,24 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) endif prec%iprcparm(psb_ilu_fill_in_) = val + case (psb_ilu_ialg_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_ialg_) = val + + case (psb_ilu_scale_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_ilu_scale_) = val + case default write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' @@ -609,3 +830,57 @@ subroutine psb_z_bjac_precseti(prec,what,val,info) 9999 call psb_error_handler(err_act) return end subroutine psb_z_bjac_precseti + +subroutine psb_z_bjac_precsetr(prec,what,val,info) + + use psb_base_mod + use psb_z_bjacprec, psb_protect_name => psb_z_bjac_precsetr + Implicit None + + class(psb_z_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + integer(psb_ipk_) :: err_act, nrow + character(len=20) :: name='z_bjac_precset' + + call psb_erractionsave(err_act) + + info = psb_success_ + if (.not.allocated(prec%iprcparm)) then + info = 1124 + call psb_errpush(info,name,a_err="preconditioner") + goto 9999 + end if + + select case(what) + case (psb_f_type_) + if (prec%iprcparm(psb_p_type_) /= psb_bjac_) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%iprcparm(psb_f_type_) = val + + case (psb_fact_eps_) + if ((prec%iprcparm(psb_p_type_) /= psb_bjac_).or.& + & (prec%iprcparm(psb_f_type_) /= psb_f_ilu_n_)) then + write(psb_err_unit,*) 'WHAT is invalid for current preconditioner ',& + & prec%iprcparm(psb_p_type_),& + & 'ignoring user specification' + return + endif + prec%rprcparm(psb_fact_eps_) = val + + case default + write(psb_err_unit,*) 'WHAT is invalid, ignoring user specification' + + end select + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return +end subroutine psb_z_bjac_precsetr diff --git a/prec/psb_c_bjacprec.f90 b/prec/psb_c_bjacprec.f90 index 2a46a6df..7a232c46 100644 --- a/prec/psb_c_bjacprec.f90 +++ b/prec/psb_c_bjacprec.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,15 +27,16 @@ ! 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 psb_c_bjacprec use psb_c_base_prec_mod use psb_c_ilu_fact_mod - + type, extends(psb_c_base_prec_type) :: psb_c_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) type(psb_cspmat_type), allocatable :: av(:) type(psb_c_vect_type), allocatable :: dv, wrk(:) contains @@ -44,6 +45,7 @@ module psb_c_bjacprec procedure, pass(prec) :: precbld => psb_c_bjac_precbld procedure, pass(prec) :: precinit => psb_c_bjac_precinit procedure, pass(prec) :: precseti => psb_c_bjac_precseti + procedure, pass(prec) :: precsetr => psb_c_bjac_precsetr procedure, pass(prec) :: precdescr => psb_c_bjac_precdescr procedure, pass(prec) :: dump => psb_c_bjac_dump procedure, pass(prec) :: clone => psb_c_bjac_clone @@ -56,14 +58,14 @@ module psb_c_bjacprec end type psb_c_bjac_prec_type private :: psb_c_bjac_sizeof, psb_c_bjac_precdescr, psb_c_bjac_get_nzeros - + character(len=15), parameter, private :: & & fact_names(0:2)=(/'None ','ILU(n) ',& & 'ILU(eps) '/) - - interface + + interface subroutine psb_c_bjac_dump(prec,info,prefix,head) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ class(psb_c_bjac_prec_type), intent(in) :: prec @@ -72,7 +74,7 @@ module psb_c_bjacprec end subroutine psb_c_bjac_dump end interface - interface + interface subroutine psb_c_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -89,7 +91,7 @@ module psb_c_bjacprec interface subroutine psb_c_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ - + type(psb_desc_type),intent(in) :: desc_data class(psb_c_bjac_prec_type), intent(inout) :: prec complex(psb_spk_),intent(in) :: alpha,beta @@ -100,7 +102,7 @@ module psb_c_bjacprec complex(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_c_bjac_apply end interface - + interface subroutine psb_c_bjac_precinit(prec,info) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ @@ -108,7 +110,7 @@ module psb_c_bjacprec integer(psb_ipk_), intent(out) :: info end subroutine psb_c_bjac_precinit end interface - + interface subroutine psb_c_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_, & @@ -123,24 +125,34 @@ module psb_c_bjacprec class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine psb_c_bjac_precbld end interface - + interface subroutine psb_c_bjac_precseti(prec,what,val,info) import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ class(psb_c_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info end subroutine psb_c_bjac_precseti end interface - + + interface + subroutine psb_c_bjac_precsetr(prec,what,val,info) + import :: psb_ipk_, psb_desc_type, psb_c_bjac_prec_type, psb_c_vect_type, psb_spk_ + class(psb_c_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + end subroutine psb_c_bjac_precsetr + end interface + contains subroutine psb_c_bjac_precdescr(prec,iout,root) use psb_penv_mod use psb_error_mod - implicit none + implicit none class(psb_c_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout @@ -153,19 +165,19 @@ contains call psb_erractionsave(err_act) info = psb_success_ - - if (present(iout)) then + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (present(root)) then + if (present(root)) then root_ = root else root_ = psb_root_ end if - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -179,26 +191,26 @@ contains & write(iout_,*) trim(prec%desc_prefix()),' ',& & 'Block Jacobi with: ',& & fact_names(prec%iprcparm(psb_f_type_)) - + call psb_erractionrestore(err_act) return 9999 call psb_error_handler(err_act) return - + end subroutine psb_c_bjac_precdescr function psb_c_bjac_sizeof(prec) result(val) class(psb_c_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + (2*psb_sizeof_sp) * prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%sizeof() val = val + prec%av(psb_u_pr_)%sizeof() endif @@ -209,12 +221,12 @@ contains class(psb_c_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%get_nzeros() val = val + prec%av(psb_u_pr_)%get_nzeros() endif @@ -235,18 +247,18 @@ contains call psb_erractionsave(err_act) info = psb_success_ - if (allocated(prec%av)) then - do i=1,size(prec%av) + if (allocated(prec%av)) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then call prec%dv%free(info) if (info == 0) deallocate(prec%dv,stat=info) end if - if (allocated(prec%iprcparm)) then + if (allocated(prec%iprcparm)) then deallocate(prec%iprcparm,stat=info) end if call psb_erractionrestore(err_act) @@ -282,19 +294,19 @@ contains & allocate(psb_c_bjac_prec_type :: precout, stat=info) if (info /= 0) goto 9999 select type(pout => precout) - type is (psb_c_bjac_prec_type) + type is (psb_c_bjac_prec_type) call pout%set_ctxt(prec%get_ctxt()) - if (allocated(prec%av)) then + if (allocated(prec%av)) then allocate(pout%av(size(prec%av)),stat=info) - do i=1,size(prec%av) + do i=1,size(prec%av) if (info /= psb_success_) exit call prec%av(i)%clone(pout%av(i),info) enddo if (info /= psb_success_) goto 9999 end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then allocate(pout%dv,stat=info) if (info == 0) call prec%dv%clone(pout%dv,info) end if @@ -315,7 +327,7 @@ contains subroutine psb_c_bjac_allocate_wrk(prec,info,vmold,desc) use psb_base_mod implicit none - + ! Arguments class(psb_c_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -325,11 +337,11 @@ contains ! Local variables integer(psb_ipk_) :: err_act, i character(len=20) :: name - + info=psb_success_ name = 'psb_c_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 if (allocated(prec%wrk)) then if (size(prec%wrk)<2) then @@ -342,11 +354,11 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + if (.not.allocated(prec%wrk)) then - if (.not.present(desc)) then + if (.not.present(desc)) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 - end if + end if allocate(prec%wrk(2),stat=info) do i=1, 2 if (info == 0) call psb_geall(prec%wrk(i),desc,info) @@ -356,19 +368,19 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_c_bjac_allocate_wrk subroutine psb_c_bjac_free_wrk(prec,info) use psb_base_mod implicit none - + ! Arguments class(psb_c_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -377,14 +389,14 @@ contains integer(psb_ipk_) :: err_act integer(psb_ipk_) :: i character(len=20) :: name - + info=psb_success_ name = 'psb_c_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 - info = psb_success_ + info = psb_success_ if (allocated(prec%wrk)) then do i=1,size(prec%wrk) if (info == 0) call prec%wrk(i)%free(info) @@ -394,29 +406,29 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_c_bjac_free_wrk function psb_c_bjac_is_allocated_wrk(prec) result(res) use psb_base_mod implicit none - + ! Arguments class(psb_c_bjac_prec_type), intent(in) :: prec logical :: res - ! In the base version we can say yes, because + ! In the base version we can say yes, because ! there is nothing to allocate res = allocated(prec%wrk) - + end function psb_c_bjac_is_allocated_wrk - + end module psb_c_bjacprec diff --git a/prec/psb_d_bjacprec.f90 b/prec/psb_d_bjacprec.f90 index 06279ae1..e8c8e47a 100644 --- a/prec/psb_d_bjacprec.f90 +++ b/prec/psb_d_bjacprec.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,15 +27,16 @@ ! 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 psb_d_bjacprec use psb_d_base_prec_mod use psb_d_ilu_fact_mod - + type, extends(psb_d_base_prec_type) :: psb_d_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) type(psb_dspmat_type), allocatable :: av(:) type(psb_d_vect_type), allocatable :: dv, wrk(:) contains @@ -44,6 +45,7 @@ module psb_d_bjacprec procedure, pass(prec) :: precbld => psb_d_bjac_precbld procedure, pass(prec) :: precinit => psb_d_bjac_precinit procedure, pass(prec) :: precseti => psb_d_bjac_precseti + procedure, pass(prec) :: precsetr => psb_d_bjac_precsetr procedure, pass(prec) :: precdescr => psb_d_bjac_precdescr procedure, pass(prec) :: dump => psb_d_bjac_dump procedure, pass(prec) :: clone => psb_d_bjac_clone @@ -56,14 +58,14 @@ module psb_d_bjacprec end type psb_d_bjac_prec_type private :: psb_d_bjac_sizeof, psb_d_bjac_precdescr, psb_d_bjac_get_nzeros - + character(len=15), parameter, private :: & & fact_names(0:2)=(/'None ','ILU(n) ',& & 'ILU(eps) '/) - - interface + + interface subroutine psb_d_bjac_dump(prec,info,prefix,head) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ class(psb_d_bjac_prec_type), intent(in) :: prec @@ -72,7 +74,7 @@ module psb_d_bjacprec end subroutine psb_d_bjac_dump end interface - interface + interface subroutine psb_d_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -89,7 +91,7 @@ module psb_d_bjacprec interface subroutine psb_d_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ - + type(psb_desc_type),intent(in) :: desc_data class(psb_d_bjac_prec_type), intent(inout) :: prec real(psb_dpk_),intent(in) :: alpha,beta @@ -100,7 +102,7 @@ module psb_d_bjacprec real(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_d_bjac_apply end interface - + interface subroutine psb_d_bjac_precinit(prec,info) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ @@ -108,7 +110,7 @@ module psb_d_bjacprec integer(psb_ipk_), intent(out) :: info end subroutine psb_d_bjac_precinit end interface - + interface subroutine psb_d_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_, & @@ -123,24 +125,34 @@ module psb_d_bjacprec class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine psb_d_bjac_precbld end interface - + interface subroutine psb_d_bjac_precseti(prec,what,val,info) import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ class(psb_d_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info end subroutine psb_d_bjac_precseti end interface - + + interface + subroutine psb_d_bjac_precsetr(prec,what,val,info) + import :: psb_ipk_, psb_desc_type, psb_d_bjac_prec_type, psb_d_vect_type, psb_dpk_ + class(psb_d_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + end subroutine psb_d_bjac_precsetr + end interface + contains subroutine psb_d_bjac_precdescr(prec,iout,root) use psb_penv_mod use psb_error_mod - implicit none + implicit none class(psb_d_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout @@ -153,19 +165,19 @@ contains call psb_erractionsave(err_act) info = psb_success_ - - if (present(iout)) then + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (present(root)) then + if (present(root)) then root_ = root else root_ = psb_root_ end if - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -179,26 +191,26 @@ contains & write(iout_,*) trim(prec%desc_prefix()),' ',& & 'Block Jacobi with: ',& & fact_names(prec%iprcparm(psb_f_type_)) - + call psb_erractionrestore(err_act) return 9999 call psb_error_handler(err_act) return - + end subroutine psb_d_bjac_precdescr function psb_d_bjac_sizeof(prec) result(val) class(psb_d_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + psb_sizeof_dp * prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%sizeof() val = val + prec%av(psb_u_pr_)%sizeof() endif @@ -209,12 +221,12 @@ contains class(psb_d_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%get_nzeros() val = val + prec%av(psb_u_pr_)%get_nzeros() endif @@ -235,18 +247,18 @@ contains call psb_erractionsave(err_act) info = psb_success_ - if (allocated(prec%av)) then - do i=1,size(prec%av) + if (allocated(prec%av)) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then call prec%dv%free(info) if (info == 0) deallocate(prec%dv,stat=info) end if - if (allocated(prec%iprcparm)) then + if (allocated(prec%iprcparm)) then deallocate(prec%iprcparm,stat=info) end if call psb_erractionrestore(err_act) @@ -282,19 +294,19 @@ contains & allocate(psb_d_bjac_prec_type :: precout, stat=info) if (info /= 0) goto 9999 select type(pout => precout) - type is (psb_d_bjac_prec_type) + type is (psb_d_bjac_prec_type) call pout%set_ctxt(prec%get_ctxt()) - if (allocated(prec%av)) then + if (allocated(prec%av)) then allocate(pout%av(size(prec%av)),stat=info) - do i=1,size(prec%av) + do i=1,size(prec%av) if (info /= psb_success_) exit call prec%av(i)%clone(pout%av(i),info) enddo if (info /= psb_success_) goto 9999 end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then allocate(pout%dv,stat=info) if (info == 0) call prec%dv%clone(pout%dv,info) end if @@ -315,7 +327,7 @@ contains subroutine psb_d_bjac_allocate_wrk(prec,info,vmold,desc) use psb_base_mod implicit none - + ! Arguments class(psb_d_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -325,11 +337,11 @@ contains ! Local variables integer(psb_ipk_) :: err_act, i character(len=20) :: name - + info=psb_success_ name = 'psb_d_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 if (allocated(prec%wrk)) then if (size(prec%wrk)<2) then @@ -342,11 +354,11 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + if (.not.allocated(prec%wrk)) then - if (.not.present(desc)) then + if (.not.present(desc)) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 - end if + end if allocate(prec%wrk(2),stat=info) do i=1, 2 if (info == 0) call psb_geall(prec%wrk(i),desc,info) @@ -356,19 +368,19 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_d_bjac_allocate_wrk subroutine psb_d_bjac_free_wrk(prec,info) use psb_base_mod implicit none - + ! Arguments class(psb_d_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -377,14 +389,14 @@ contains integer(psb_ipk_) :: err_act integer(psb_ipk_) :: i character(len=20) :: name - + info=psb_success_ name = 'psb_d_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 - info = psb_success_ + info = psb_success_ if (allocated(prec%wrk)) then do i=1,size(prec%wrk) if (info == 0) call prec%wrk(i)%free(info) @@ -394,29 +406,29 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_d_bjac_free_wrk function psb_d_bjac_is_allocated_wrk(prec) result(res) use psb_base_mod implicit none - + ! Arguments class(psb_d_bjac_prec_type), intent(in) :: prec logical :: res - ! In the base version we can say yes, because + ! In the base version we can say yes, because ! there is nothing to allocate res = allocated(prec%wrk) - + end function psb_d_bjac_is_allocated_wrk - + end module psb_d_bjacprec diff --git a/prec/psb_prec_const_mod.f90 b/prec/psb_prec_const_mod.f90 index f7b32d2f..0e0e019b 100644 --- a/prec/psb_prec_const_mod.f90 +++ b/prec/psb_prec_const_mod.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,8 +27,8 @@ ! 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 to define PREC_DATA, !! !! structure for preconditioning. !! @@ -43,19 +43,21 @@ module psb_prec_const_mod ! Entries in iprcparm: preconditioner type, factorization type, ! prolongation type, restriction type, renumbering algorithm, - ! number of overlap layers, pointer to SuperLU factors, - ! levels of fill in for ILU(N), + ! number of overlap layers, pointer to SuperLU factors, + ! levels of fill in for ILU(N), integer(psb_ipk_), parameter :: psb_p_type_=1, psb_f_type_=2 integer(psb_ipk_), parameter :: psb_ilu_fill_in_=8 + integer(psb_ipk_), parameter :: psb_ilu_ialg_=9 !Renumbering. SEE BELOW integer(psb_ipk_), parameter :: psb_renum_none_=0, psb_renum_glb_=1, psb_renum_gps_=2 integer(psb_ipk_), parameter :: psb_ifpsz=10 ! Entries in rprcparm: ILU(E) epsilon, smoother omega + integer(psb_ipk_), parameter :: psb_ilu_scale_=7 integer(psb_ipk_), parameter :: psb_fact_eps_=1 integer(psb_ipk_), parameter :: psb_rfpsz=4 - ! Factorization types: none, ILU(N), ILU(E) - integer(psb_ipk_), parameter :: psb_f_none_=0,psb_f_ilu_n_=1 - ! Fields for sparse matrices ensembles: + ! Factorization types: none, ILU(0), ILU(N), ILU(N,E) + integer(psb_ipk_), parameter :: psb_f_none_=0,psb_f_ilu_n_=1,psb_f_ilu_k_=2,psb_f_ilu_t_=3 + ! Fields for sparse matrices ensembles: integer(psb_ipk_), parameter :: psb_l_pr_=1, psb_u_pr_=2, psb_bp_ilu_avsz=2 integer(psb_ipk_), parameter :: psb_max_avsz=psb_bp_ilu_avsz @@ -65,11 +67,11 @@ module psb_prec_const_mod integer(psb_ipk_), parameter :: psb_ilu_scale_none_ = 0 integer(psb_ipk_), parameter :: psb_ilu_scale_maxval_ = 1 integer(psb_ipk_), parameter :: psb_ilu_scale_diag_ = 2 - integer(psb_ipk_), parameter :: psb_ilu_scale_arwsum_ = 3 + integer(psb_ipk_), parameter :: psb_ilu_scale_arwsum_ = 3 integer(psb_ipk_), parameter :: psb_ilu_scale_aclsum_ = 4 integer(psb_ipk_), parameter :: psb_ilu_scale_arcsum_ = 5 - + interface psb_check_def module procedure psb_icheck_def, psb_scheck_def, psb_dcheck_def @@ -87,9 +89,9 @@ contains select case(iprec) case(psb_noprec_) pr_to_str='NOPREC' - case(psb_diag_) + case(psb_diag_) pr_to_str='DIAG' - case(psb_bjac_) + case(psb_bjac_) pr_to_str='BJAC' case default pr_to_str='???' @@ -125,7 +127,7 @@ contains integer(psb_ipk_), intent(inout) :: ip integer(psb_ipk_), intent(in) :: id character(len=*), intent(in) :: name - interface + interface function is_legal(i) import :: psb_ipk_ integer(psb_ipk_), intent(in) :: i @@ -133,7 +135,7 @@ contains end function is_legal end interface - if (.not.is_legal(ip)) then + if (.not.is_legal(ip)) then write(psb_err_unit,*) 'Illegal value for ',name,' :',ip, '. defaulting to ',id ip = id end if @@ -143,7 +145,7 @@ contains real(psb_spk_), intent(inout) :: ip real(psb_spk_), intent(in) :: id character(len=*), intent(in) :: name - interface + interface function is_legal(i) import :: psb_ipk_, psb_spk_ real(psb_spk_), intent(in) :: i @@ -151,7 +153,7 @@ contains end function is_legal end interface - if (.not.is_legal(ip)) then + if (.not.is_legal(ip)) then write(psb_err_unit,*) 'Illegal value for ',name,' :',ip, '. defaulting to ',id ip = id end if @@ -161,7 +163,7 @@ contains real(psb_dpk_), intent(inout) :: ip real(psb_dpk_), intent(in) :: id character(len=*), intent(in) :: name - interface + interface function is_legal(i) import :: psb_ipk_, psb_spk_, psb_dpk_ real(psb_dpk_), intent(in) :: i @@ -169,7 +171,7 @@ contains end function is_legal end interface - if (.not.is_legal(ip)) then + if (.not.is_legal(ip)) then write(psb_err_unit,*) 'Illegal value for ',name,' :',ip, '. defaulting to ',id ip = id end if diff --git a/prec/psb_s_bjacprec.f90 b/prec/psb_s_bjacprec.f90 index 25ad7642..2bcf2e02 100644 --- a/prec/psb_s_bjacprec.f90 +++ b/prec/psb_s_bjacprec.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,15 +27,16 @@ ! 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 psb_s_bjacprec use psb_s_base_prec_mod use psb_s_ilu_fact_mod - + type, extends(psb_s_base_prec_type) :: psb_s_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) + real(psb_spk_), allocatable :: rprcparm(:) type(psb_sspmat_type), allocatable :: av(:) type(psb_s_vect_type), allocatable :: dv, wrk(:) contains @@ -44,6 +45,7 @@ module psb_s_bjacprec procedure, pass(prec) :: precbld => psb_s_bjac_precbld procedure, pass(prec) :: precinit => psb_s_bjac_precinit procedure, pass(prec) :: precseti => psb_s_bjac_precseti + procedure, pass(prec) :: precsetr => psb_s_bjac_precsetr procedure, pass(prec) :: precdescr => psb_s_bjac_precdescr procedure, pass(prec) :: dump => psb_s_bjac_dump procedure, pass(prec) :: clone => psb_s_bjac_clone @@ -56,14 +58,14 @@ module psb_s_bjacprec end type psb_s_bjac_prec_type private :: psb_s_bjac_sizeof, psb_s_bjac_precdescr, psb_s_bjac_get_nzeros - + character(len=15), parameter, private :: & & fact_names(0:2)=(/'None ','ILU(n) ',& & 'ILU(eps) '/) - - interface + + interface subroutine psb_s_bjac_dump(prec,info,prefix,head) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ class(psb_s_bjac_prec_type), intent(in) :: prec @@ -72,7 +74,7 @@ module psb_s_bjacprec end subroutine psb_s_bjac_dump end interface - interface + interface subroutine psb_s_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ type(psb_desc_type),intent(in) :: desc_data @@ -89,7 +91,7 @@ module psb_s_bjacprec interface subroutine psb_s_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ - + type(psb_desc_type),intent(in) :: desc_data class(psb_s_bjac_prec_type), intent(inout) :: prec real(psb_spk_),intent(in) :: alpha,beta @@ -100,7 +102,7 @@ module psb_s_bjacprec real(psb_spk_),intent(inout), optional, target :: work(:) end subroutine psb_s_bjac_apply end interface - + interface subroutine psb_s_bjac_precinit(prec,info) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ @@ -108,7 +110,7 @@ module psb_s_bjacprec integer(psb_ipk_), intent(out) :: info end subroutine psb_s_bjac_precinit end interface - + interface subroutine psb_s_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_, & @@ -123,24 +125,34 @@ module psb_s_bjacprec class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine psb_s_bjac_precbld end interface - + interface subroutine psb_s_bjac_precseti(prec,what,val,info) import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ class(psb_s_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info end subroutine psb_s_bjac_precseti end interface - + + interface + subroutine psb_s_bjac_precsetr(prec,what,val,info) + import :: psb_ipk_, psb_desc_type, psb_s_bjac_prec_type, psb_s_vect_type, psb_spk_ + class(psb_s_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_spk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + end subroutine psb_s_bjac_precsetr + end interface + contains subroutine psb_s_bjac_precdescr(prec,iout,root) use psb_penv_mod use psb_error_mod - implicit none + implicit none class(psb_s_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout @@ -153,19 +165,19 @@ contains call psb_erractionsave(err_act) info = psb_success_ - - if (present(iout)) then + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (present(root)) then + if (present(root)) then root_ = root else root_ = psb_root_ end if - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -179,26 +191,26 @@ contains & write(iout_,*) trim(prec%desc_prefix()),' ',& & 'Block Jacobi with: ',& & fact_names(prec%iprcparm(psb_f_type_)) - + call psb_erractionrestore(err_act) return 9999 call psb_error_handler(err_act) return - + end subroutine psb_s_bjac_precdescr function psb_s_bjac_sizeof(prec) result(val) class(psb_s_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + psb_sizeof_sp * prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%sizeof() val = val + prec%av(psb_u_pr_)%sizeof() endif @@ -209,12 +221,12 @@ contains class(psb_s_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%get_nzeros() val = val + prec%av(psb_u_pr_)%get_nzeros() endif @@ -235,18 +247,18 @@ contains call psb_erractionsave(err_act) info = psb_success_ - if (allocated(prec%av)) then - do i=1,size(prec%av) + if (allocated(prec%av)) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then call prec%dv%free(info) if (info == 0) deallocate(prec%dv,stat=info) end if - if (allocated(prec%iprcparm)) then + if (allocated(prec%iprcparm)) then deallocate(prec%iprcparm,stat=info) end if call psb_erractionrestore(err_act) @@ -282,19 +294,19 @@ contains & allocate(psb_s_bjac_prec_type :: precout, stat=info) if (info /= 0) goto 9999 select type(pout => precout) - type is (psb_s_bjac_prec_type) + type is (psb_s_bjac_prec_type) call pout%set_ctxt(prec%get_ctxt()) - if (allocated(prec%av)) then + if (allocated(prec%av)) then allocate(pout%av(size(prec%av)),stat=info) - do i=1,size(prec%av) + do i=1,size(prec%av) if (info /= psb_success_) exit call prec%av(i)%clone(pout%av(i),info) enddo if (info /= psb_success_) goto 9999 end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then allocate(pout%dv,stat=info) if (info == 0) call prec%dv%clone(pout%dv,info) end if @@ -315,7 +327,7 @@ contains subroutine psb_s_bjac_allocate_wrk(prec,info,vmold,desc) use psb_base_mod implicit none - + ! Arguments class(psb_s_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -325,11 +337,11 @@ contains ! Local variables integer(psb_ipk_) :: err_act, i character(len=20) :: name - + info=psb_success_ name = 'psb_s_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 if (allocated(prec%wrk)) then if (size(prec%wrk)<2) then @@ -342,11 +354,11 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + if (.not.allocated(prec%wrk)) then - if (.not.present(desc)) then + if (.not.present(desc)) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 - end if + end if allocate(prec%wrk(2),stat=info) do i=1, 2 if (info == 0) call psb_geall(prec%wrk(i),desc,info) @@ -356,19 +368,19 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_s_bjac_allocate_wrk subroutine psb_s_bjac_free_wrk(prec,info) use psb_base_mod implicit none - + ! Arguments class(psb_s_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -377,14 +389,14 @@ contains integer(psb_ipk_) :: err_act integer(psb_ipk_) :: i character(len=20) :: name - + info=psb_success_ name = 'psb_s_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 - info = psb_success_ + info = psb_success_ if (allocated(prec%wrk)) then do i=1,size(prec%wrk) if (info == 0) call prec%wrk(i)%free(info) @@ -394,29 +406,29 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_s_bjac_free_wrk function psb_s_bjac_is_allocated_wrk(prec) result(res) use psb_base_mod implicit none - + ! Arguments class(psb_s_bjac_prec_type), intent(in) :: prec logical :: res - ! In the base version we can say yes, because + ! In the base version we can say yes, because ! there is nothing to allocate res = allocated(prec%wrk) - + end function psb_s_bjac_is_allocated_wrk - + end module psb_s_bjacprec diff --git a/prec/psb_z_bjacprec.f90 b/prec/psb_z_bjacprec.f90 index 8ca5616a..de9b3518 100644 --- a/prec/psb_z_bjacprec.f90 +++ b/prec/psb_z_bjacprec.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,15 +27,16 @@ ! 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 psb_z_bjacprec use psb_z_base_prec_mod use psb_z_ilu_fact_mod - + type, extends(psb_z_base_prec_type) :: psb_z_bjac_prec_type integer(psb_ipk_), allocatable :: iprcparm(:) + real(psb_dpk_), allocatable :: rprcparm(:) type(psb_zspmat_type), allocatable :: av(:) type(psb_z_vect_type), allocatable :: dv, wrk(:) contains @@ -44,6 +45,7 @@ module psb_z_bjacprec procedure, pass(prec) :: precbld => psb_z_bjac_precbld procedure, pass(prec) :: precinit => psb_z_bjac_precinit procedure, pass(prec) :: precseti => psb_z_bjac_precseti + procedure, pass(prec) :: precsetr => psb_z_bjac_precsetr procedure, pass(prec) :: precdescr => psb_z_bjac_precdescr procedure, pass(prec) :: dump => psb_z_bjac_dump procedure, pass(prec) :: clone => psb_z_bjac_clone @@ -56,14 +58,14 @@ module psb_z_bjacprec end type psb_z_bjac_prec_type private :: psb_z_bjac_sizeof, psb_z_bjac_precdescr, psb_z_bjac_get_nzeros - + character(len=15), parameter, private :: & & fact_names(0:2)=(/'None ','ILU(n) ',& & 'ILU(eps) '/) - - interface + + interface subroutine psb_z_bjac_dump(prec,info,prefix,head) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ class(psb_z_bjac_prec_type), intent(in) :: prec @@ -72,7 +74,7 @@ module psb_z_bjacprec end subroutine psb_z_bjac_dump end interface - interface + interface subroutine psb_z_bjac_apply_vect(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ type(psb_desc_type),intent(in) :: desc_data @@ -89,7 +91,7 @@ module psb_z_bjacprec interface subroutine psb_z_bjac_apply(alpha,prec,x,beta,y,desc_data,info,trans,work) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ - + type(psb_desc_type),intent(in) :: desc_data class(psb_z_bjac_prec_type), intent(inout) :: prec complex(psb_dpk_),intent(in) :: alpha,beta @@ -100,7 +102,7 @@ module psb_z_bjacprec complex(psb_dpk_),intent(inout), optional, target :: work(:) end subroutine psb_z_bjac_apply end interface - + interface subroutine psb_z_bjac_precinit(prec,info) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ @@ -108,7 +110,7 @@ module psb_z_bjacprec integer(psb_ipk_), intent(out) :: info end subroutine psb_z_bjac_precinit end interface - + interface subroutine psb_z_bjac_precbld(a,desc_a,prec,info,amold,vmold,imold) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_, & @@ -123,24 +125,34 @@ module psb_z_bjacprec class(psb_i_base_vect_type), intent(in), optional :: imold end subroutine psb_z_bjac_precbld end interface - + interface subroutine psb_z_bjac_precseti(prec,what,val,info) import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ class(psb_z_bjac_prec_type),intent(inout) :: prec - integer(psb_ipk_), intent(in) :: what - integer(psb_ipk_), intent(in) :: val + integer(psb_ipk_), intent(in) :: what + integer(psb_ipk_), intent(in) :: val integer(psb_ipk_), intent(out) :: info end subroutine psb_z_bjac_precseti end interface - + + interface + subroutine psb_z_bjac_precsetr(prec,what,val,info) + import :: psb_ipk_, psb_desc_type, psb_z_bjac_prec_type, psb_z_vect_type, psb_dpk_ + class(psb_z_bjac_prec_type),intent(inout) :: prec + integer(psb_ipk_), intent(in) :: what + real(psb_dpk_), intent(in) :: val + integer(psb_ipk_), intent(out) :: info + end subroutine psb_z_bjac_precsetr + end interface + contains subroutine psb_z_bjac_precdescr(prec,iout,root) use psb_penv_mod use psb_error_mod - implicit none + implicit none class(psb_z_bjac_prec_type), intent(in) :: prec integer(psb_ipk_), intent(in), optional :: iout @@ -153,19 +165,19 @@ contains call psb_erractionsave(err_act) info = psb_success_ - - if (present(iout)) then + + if (present(iout)) then iout_ = iout else - iout_ = 6 + iout_ = 6 end if - if (present(root)) then + if (present(root)) then root_ = root else root_ = psb_root_ end if - if (.not.allocated(prec%iprcparm)) then + if (.not.allocated(prec%iprcparm)) then info = 1124 call psb_errpush(info,name,a_err="preconditioner") goto 9999 @@ -179,26 +191,26 @@ contains & write(iout_,*) trim(prec%desc_prefix()),' ',& & 'Block Jacobi with: ',& & fact_names(prec%iprcparm(psb_f_type_)) - + call psb_erractionrestore(err_act) return 9999 call psb_error_handler(err_act) return - + end subroutine psb_z_bjac_precdescr function psb_z_bjac_sizeof(prec) result(val) class(psb_z_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + (2*psb_sizeof_dp) * prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%sizeof() val = val + prec%av(psb_u_pr_)%sizeof() endif @@ -209,12 +221,12 @@ contains class(psb_z_bjac_prec_type), intent(in) :: prec integer(psb_epk_) :: val - + val = 0 - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then val = val + prec%dv%get_nrows() endif - if (allocated(prec%av)) then + if (allocated(prec%av)) then val = val + prec%av(psb_l_pr_)%get_nzeros() val = val + prec%av(psb_u_pr_)%get_nzeros() endif @@ -235,18 +247,18 @@ contains call psb_erractionsave(err_act) info = psb_success_ - if (allocated(prec%av)) then - do i=1,size(prec%av) + if (allocated(prec%av)) then + do i=1,size(prec%av) call prec%av(i)%free() enddo deallocate(prec%av,stat=info) end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then call prec%dv%free(info) if (info == 0) deallocate(prec%dv,stat=info) end if - if (allocated(prec%iprcparm)) then + if (allocated(prec%iprcparm)) then deallocate(prec%iprcparm,stat=info) end if call psb_erractionrestore(err_act) @@ -282,19 +294,19 @@ contains & allocate(psb_z_bjac_prec_type :: precout, stat=info) if (info /= 0) goto 9999 select type(pout => precout) - type is (psb_z_bjac_prec_type) + type is (psb_z_bjac_prec_type) call pout%set_ctxt(prec%get_ctxt()) - if (allocated(prec%av)) then + if (allocated(prec%av)) then allocate(pout%av(size(prec%av)),stat=info) - do i=1,size(prec%av) + do i=1,size(prec%av) if (info /= psb_success_) exit call prec%av(i)%clone(pout%av(i),info) enddo if (info /= psb_success_) goto 9999 end if - if (allocated(prec%dv)) then + if (allocated(prec%dv)) then allocate(pout%dv,stat=info) if (info == 0) call prec%dv%clone(pout%dv,info) end if @@ -315,7 +327,7 @@ contains subroutine psb_z_bjac_allocate_wrk(prec,info,vmold,desc) use psb_base_mod implicit none - + ! Arguments class(psb_z_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -325,11 +337,11 @@ contains ! Local variables integer(psb_ipk_) :: err_act, i character(len=20) :: name - + info=psb_success_ name = 'psb_z_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 if (allocated(prec%wrk)) then if (size(prec%wrk)<2) then @@ -342,11 +354,11 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + if (.not.allocated(prec%wrk)) then - if (.not.present(desc)) then + if (.not.present(desc)) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="no desc?"); goto 9999 - end if + end if allocate(prec%wrk(2),stat=info) do i=1, 2 if (info == 0) call psb_geall(prec%wrk(i),desc,info) @@ -356,19 +368,19 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="allocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_z_bjac_allocate_wrk subroutine psb_z_bjac_free_wrk(prec,info) use psb_base_mod implicit none - + ! Arguments class(psb_z_bjac_prec_type), intent(inout) :: prec integer(psb_ipk_), intent(out) :: info @@ -377,14 +389,14 @@ contains integer(psb_ipk_) :: err_act integer(psb_ipk_) :: i character(len=20) :: name - + info=psb_success_ name = 'psb_z_allocate_wrk' call psb_erractionsave(err_act) - + if (psb_get_errstatus().ne.0) goto 9999 - info = psb_success_ + info = psb_success_ if (allocated(prec%wrk)) then do i=1,size(prec%wrk) if (info == 0) call prec%wrk(i)%free(info) @@ -394,29 +406,29 @@ contains if (info /= 0) then info = psb_err_internal_error_; call psb_errpush(info,name,a_err="deallocate"); goto 9999 end if - + call psb_erractionrestore(err_act) return - + 9999 call psb_error_handler(err_act) return - + end subroutine psb_z_bjac_free_wrk function psb_z_bjac_is_allocated_wrk(prec) result(res) use psb_base_mod implicit none - + ! Arguments class(psb_z_bjac_prec_type), intent(in) :: prec logical :: res - ! In the base version we can say yes, because + ! In the base version we can say yes, because ! there is nothing to allocate res = allocated(prec%wrk) - + end function psb_z_bjac_is_allocated_wrk - + end module psb_z_bjacprec diff --git a/test/pargen/psb_d_pde3d.f90 b/test/pargen/psb_d_pde3d.f90 index 429e9a0e..1a07fffd 100644 --- a/test/pargen/psb_d_pde3d.f90 +++ b/test/pargen/psb_d_pde3d.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,23 +27,23 @@ ! 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. -! -! +! +! ! File: psb_d_pde3d.f90 ! ! Program: psb_d_pde3d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 3d ! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f -! dxdx dydy dzdz dx dy dz +! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -60,37 +60,37 @@ ! module psb_d_pde3d_mod - + use psb_base_mod, only : psb_dpk_, psb_ipk_, psb_lpk_, psb_desc_type,& & psb_dspmat_type, psb_d_vect_type, dzero,& & psb_d_base_sparse_mat, psb_d_base_vect_type, & & psb_i_base_vect_type, psb_l_base_vect_type - interface + interface function d_func_3d(x,y,z) result(val) import :: psb_dpk_ real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_) :: val end function d_func_3d - end interface + end interface interface psb_gen_pde3d module procedure psb_d_gen_pde3d end interface psb_gen_pde3d - + contains function d_null_func_3d(x,y,z) result(val) real(psb_dpk_), intent(in) :: x,y,z real(psb_dpk_) :: val - + val = dzero end function d_null_func_3d ! - ! functions parametrizing the differential equation - ! + ! functions parametrizing the differential equation + ! ! ! Note: b1, b2 and b3 are the coefficients of the first @@ -103,70 +103,70 @@ contains ! function b1(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: b1 real(psb_dpk_), intent(in) :: x,y,z b1=dzero end function b1 function b2(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: b2 real(psb_dpk_), intent(in) :: x,y,z b2=dzero end function b2 function b3(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: b3 - real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_), intent(in) :: x,y,z b3=dzero end function b3 function c(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: c - real(psb_dpk_), intent(in) :: x,y,z + real(psb_dpk_), intent(in) :: x,y,z c=dzero end function c function a1(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none - real(psb_dpk_) :: a1 + implicit none + real(psb_dpk_) :: a1 real(psb_dpk_), intent(in) :: x,y,z a1=done/80 end function a1 function a2(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: a2 real(psb_dpk_), intent(in) :: x,y,z a2=done/80 end function a2 function a3(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: a3 real(psb_dpk_), intent(in) :: x,y,z a3=done/80 end function a3 function g(x,y,z) use psb_base_mod, only : psb_dpk_, done, dzero - implicit none + implicit none real(psb_dpk_) :: g real(psb_dpk_), intent(in) :: x,y,z g = dzero if (x == done) then g = done - else if (x == dzero) then + else if (x == dzero) then g = exp(y**2-z**2) end if end function g - + ! ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. + ! the rhs. ! subroutine psb_d_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,info,& & f,amold,vmold,imold,partition,nrl,iv) @@ -174,13 +174,13 @@ contains use psb_util_mod ! ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) + ! + ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f - ! dxdx dydy dzdz dx dy dz + ! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions - ! u = g + ! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -196,7 +196,7 @@ contains character(len=*) :: afmt procedure(d_func_3d), optional :: f class(psb_d_base_sparse_mat), optional :: amold - class(psb_d_base_vect_type), optional :: vmold + class(psb_d_base_vect_type), optional :: vmold class(psb_i_base_vect_type), optional :: imold integer(psb_ipk_), optional :: partition, nrl,iv(:) @@ -237,7 +237,7 @@ contains call psb_info(ictxt, iam, np) - if (present(f)) then + if (present(f)) then f_ => f else f_ => d_null_func_3d @@ -257,10 +257,10 @@ contains else partition_ = 3 end if - + ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - + ! estimate of the number of non zeroes + m = (1_psb_lpk_*idim)*idim*idim n = m nnz = ((n*7)/(np)) @@ -268,8 +268,8 @@ contains t0 = psb_wtime() select case(partition_) case(1) - ! A BLOCK partition - if (present(nrl)) then + ! A BLOCK partition + if (present(nrl)) then nr = nrl else ! @@ -280,46 +280,46 @@ contains end if nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) then + call psb_sum(ictxt,nt) + if (nt /= m) then write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! First example of use of CDALL: specify for each process a number of ! contiguous rows - ! + ! call psb_cdall(ictxt,desc_a,info,nl=nr) myidx = desc_a%get_global_indices() nlr = size(myidx) case(2) ! A partition defined by the user through IV - - if (present(iv)) then + + if (present(iv)) then if (size(iv) /= m) then write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if else write(psb_err_unit,*) iam, 'Initialization error: IV not present' info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! + ! process that owns it + ! call psb_cdall(ictxt,desc_a,info,vg=iv) myidx = desc_a%get_global_indices() nlr = size(myidx) @@ -335,7 +335,7 @@ contains npz = npdims(3) allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. + ! We can reuse idx2ijk for process indices as well. call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) ! Now let's split the 3D cube in hexahedra call dist1Didx(bndx,idim,npx) @@ -345,7 +345,7 @@ contains call dist1Didx(bndz,idim,npz) mynz = bndz(iamz+1)-bndz(iamz) - ! How many indices do I own? + ! How many indices do I own? nlr = mynx*myny*mynz allocate(myidx(nlr)) ! Now, let's generate the list of indices I own @@ -369,9 +369,9 @@ contains ! ! Third example of use of CDALL: specify for each process ! the set of global indices it owns. - ! + ! call psb_cdall(ictxt,desc_a,info,vl=myidx) - + case default write(psb_err_unit,*) iam, 'Initialization error: should not get here' info = -1 @@ -380,9 +380,9 @@ contains return end select - + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess + ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) @@ -397,12 +397,12 @@ contains end if ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! allocate(val(20*nb),irow(20*nb),& &icol(20*nb),stat=info) - if (info /= psb_success_ ) then + if (info /= psb_success_ ) then info=psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 @@ -415,11 +415,11 @@ contains call psb_barrier(ictxt) t1 = psb_wtime() do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) + ib = min(nb,nlr-ii+1) icoeff = 1 do k=1,ib i=ii+k-1 - ! local matrix pointer + ! local matrix pointer glob_row=myidx(i) ! compute gridpoint coordinates call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) @@ -429,11 +429,11 @@ contains z = (iz-1)*deltah zt(k) = f_(x,y,z) ! internal point: build discretization - ! + ! ! term depending on (x-1,y,z) ! val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 - if (ix == 1) then + if (ix == 1) then zt(k) = g(dzero,y,z)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) @@ -442,19 +442,19 @@ contains endif ! term depending on (x,y-1,z) val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 - if (iy == 1) then + if (iy == 1) then zt(k) = g(x,dzero,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x,y,z-1) val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 - if (iz == 1) then + if (iz == 1) then zt(k) = g(x,y,dzero)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif @@ -462,33 +462,33 @@ contains ! term depending on (x,y,z) val(icoeff)=(2*done)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & & + c(x,y,z) - call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) irow(icoeff) = glob_row - icoeff = icoeff+1 + icoeff = icoeff+1 ! term depending on (x,y,z+1) val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 - if (iz == idim) then + if (iz == idim) then zt(k) = g(x,y,done)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x,y+1,z) val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 - if (iy == idim) then + if (iy == idim) then zt(k) = g(x,done,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x+1,y,z) val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 - if (ix==idim) then + if (ix==idim) then zt(k) = g(done,y,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif @@ -519,8 +519,8 @@ contains tcdasb = psb_wtime()-t1 call psb_barrier(ictxt) t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then + if (info == psb_success_) then + if (present(amold)) then call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) else call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) @@ -543,7 +543,7 @@ contains end if tasb = psb_wtime()-t1 call psb_barrier(ictxt) - ttot = psb_wtime() - t0 + ttot = psb_wtime() - t0 call psb_amx(ictxt,talc) call psb_amx(ictxt,tgen) @@ -585,9 +585,9 @@ program psb_d_pde3d integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_dpk_), parameter :: one = done - real(psb_dpk_) :: t1, t2, tprec + real(psb_dpk_) :: t1, t2, tprec ! sparse matrix and preconditioner type(psb_dspmat_type) :: a @@ -604,6 +604,14 @@ program psb_d_pde3d integer(psb_epk_) :: amatsize, precsize, descsize, d2size real(psb_dpk_) :: err, eps + ! Parameters for solvers in Block-Jacobi preconditioner + type ainvparms + character(len=12) :: alg, orth_alg + integer(psb_ipk_) :: fill, inv_fill + real(psb_dpk_) :: thresh, inv_thresh + end type ainvparms + type(ainvparms) :: parms + ! other variables integer(psb_ipk_) :: info, i character(len=20) :: name,ch_err @@ -615,7 +623,7 @@ program psb_d_pde3d call psb_init(ictxt) call psb_info(ictxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ictxt) stop @@ -626,21 +634,21 @@ program psb_d_pde3d ! ! Hello world ! - if (iam == psb_root_) then + if (iam == psb_root_) then write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) + call psb_gen_pde3d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -653,7 +661,7 @@ program psb_d_pde3d if (iam == psb_root_) write(psb_out_unit,'(" ")') ! ! prepare the preconditioner. - ! + ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) @@ -675,14 +683,14 @@ program psb_d_pde3d if (iam == psb_root_) write(psb_out_unit,'(" ")') call prec%descr() ! - ! iterative method parameters + ! iterative method parameters ! if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd call psb_barrier(ictxt) - t1 = psb_wtime() + t1 = psb_wtime() eps = 1.d-6 - call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) if(info /= psb_success_) then info=psb_err_from_subroutine_ @@ -712,14 +720,14 @@ program psb_d_pde3d write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err write(psb_out_unit,'("Info on exit : ",i12)')info write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize - write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize write(psb_out_unit,'("Storage format for A: ",a)') a%get_fmt() write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt() end if - ! + ! ! cleanup storage and exit ! call psb_gefree(bv,desc_a,info) @@ -745,13 +753,15 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,& + itmax,itrace,irst,ipart,parms) integer(psb_ipk_) :: ictxt character(len=*) :: kmethd, ptype, afmt integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart integer(psb_ipk_) :: np, iam integer(psb_ipk_) :: ip, inp_unit character(len=1024) :: filename + type(ainvparms) :: parms call psb_info(ictxt, iam, np) @@ -780,12 +790,12 @@ contains if (ip >= 4) then read(inp_unit,*) ipart else - ipart = 3 + ipart = 3 endif if (ip >= 5) then read(inp_unit,*) istopc else - istopc=1 + istopc=1 endif if (ip >= 6) then read(inp_unit,*) itmax @@ -802,8 +812,23 @@ contains else irst=1 endif + if (ip >= 9) then + read(psb_inp_unit,*) parms%alg + read(psb_inp_unit,*) parms%fill + read(psb_inp_unit,*) parms%inv_fill + read(psb_inp_unit,*) parms%thresh + read(psb_inp_unit,*) parms%inv_thresh + read(psb_inp_unit,*) parms%orth_alg + else + parms%alg = 'ILU' ! AINV variant: ILU,ILUT,MILU,INVK,AINVT,AORTH + parms%fill = 0 ! Fill in for forward factorization + parms%inv_fill = 1 ! Fill in for inverse factorization + parms%thresh = 1E-1_psb_dpk_ ! Threshold for forward factorization + parms%inv_thresh = 1E-1_psb_dpk_ ! Threshold for inverse factorization + parms%orth_alg = 'LLK' ! What orthogonalization algorithm? + endif - write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Solving matrix : ell1")') write(psb_out_unit,& & '("Grid dimensions : ",i4," x ",i4," x ",i4)') & & idim,idim,idim @@ -818,11 +843,28 @@ contains write(psb_out_unit,'("Unknown data distrbution, defaulting to 3D")') end select write(psb_out_unit,'("Preconditioner : ",a)') ptype + if( psb_toupper(ptype) == "BJAC" ) then + write(psb_out_unit,'("Block subsolver : ",a)') parms%alg + select case (psb_toupper(parms%alg)) + case ('ILU','ILUT','MILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + case ('INVK') + write(psb_out_unit,'("Fill : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + write(psb_out_unit,'("Invese Fill : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + case ('AINVT','AORTH') + write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + case default + write(psb_out_unit,'("Unknown diagonal solver")') + end select + end if write(psb_out_unit,'("Iterative method : ",a)') kmethd write(psb_out_unit,'(" ")') else ! wrong number of parameter, print an error message and exit - call pr_usage(izero) + call pr_usage(izero) call psb_abort(ictxt) stop 1 endif @@ -841,20 +883,26 @@ contains call psb_bcast(ictxt,itmax) call psb_bcast(ictxt,itrace) call psb_bcast(ictxt,irst) + call psb_bcast(ictxt,parms%alg) + call psb_bcast(ictxt,parms%fill) + call psb_bcast(ictxt,parms%inv_fill) + call psb_bcast(ictxt,parms%thresh) + call psb_bcast(ictxt,parms%inv_thresh) + call psb_bcast(ictxt,parms%orth_alg) return end subroutine get_parms ! - ! print an error message - ! + ! print an error message + ! subroutine pr_usage(iout) integer(psb_ipk_) :: iout write(iout,*)'incorrect parameter(s) found' write(iout,*)' usage: pde3d90 methd prec dim & - &[istop itmax itrace]' + &[istop itmax itrace]' write(iout,*)' where:' - write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' write(iout,*)' prec : bjac diag none' write(iout,*)' dim number of points along each axis' write(iout,*)' the size of the resulting linear ' @@ -862,11 +910,9 @@ contains write(iout,*)' ipart data partition 1 3 ' write(iout,*)' istop stopping criterion 1, 2 ' write(iout,*)' itmax maximum number of iterations [500] ' - write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' itrace <=0 (no tracing, default) or ' write(iout,*)' >= 1 do tracing every itrace' - write(iout,*)' iterations ' + write(iout,*)' iterations ' end subroutine pr_usage end program psb_d_pde3d - - diff --git a/test/pargen/psb_s_pde3d.f90 b/test/pargen/psb_s_pde3d.f90 index e7c7725e..d4ad1492 100644 --- a/test/pargen/psb_s_pde3d.f90 +++ b/test/pargen/psb_s_pde3d.f90 @@ -1,9 +1,9 @@ -! +! ! Parallel Sparse BLAS version 3.5 ! (C) Copyright 2006-2018 -! Salvatore Filippone -! Alfredo Buttari -! +! Salvatore Filippone +! Alfredo Buttari +! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the following conditions ! are met: @@ -15,7 +15,7 @@ ! 3. The name of the PSBLAS 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 @@ -27,23 +27,23 @@ ! 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. -! -! +! +! ! File: psb_s_pde3d.f90 ! ! Program: psb_s_pde3d ! This sample program solves a linear system obtained by discretizing a -! PDE with Dirichlet BCs. -! +! PDE with Dirichlet BCs. +! ! ! The PDE is a general second order equation in 3d ! -! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) +! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f -! dxdx dydy dzdz dx dy dz +! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions -! u = g +! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -60,37 +60,37 @@ ! module psb_s_pde3d_mod - + use psb_base_mod, only : psb_spk_, psb_ipk_, psb_lpk_, psb_desc_type,& & psb_sspmat_type, psb_s_vect_type, szero,& & psb_s_base_sparse_mat, psb_s_base_vect_type, & & psb_i_base_vect_type, psb_l_base_vect_type - interface + interface function s_func_3d(x,y,z) result(val) import :: psb_spk_ real(psb_spk_), intent(in) :: x,y,z real(psb_spk_) :: val end function s_func_3d - end interface + end interface interface psb_gen_pde3d module procedure psb_s_gen_pde3d end interface psb_gen_pde3d - + contains function s_null_func_3d(x,y,z) result(val) real(psb_spk_), intent(in) :: x,y,z real(psb_spk_) :: val - + val = szero end function s_null_func_3d ! - ! functions parametrizing the differential equation - ! + ! functions parametrizing the differential equation + ! ! ! Note: b1, b2 and b3 are the coefficients of the first @@ -103,70 +103,70 @@ contains ! function b1(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: b1 real(psb_spk_), intent(in) :: x,y,z b1=szero end function b1 function b2(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: b2 real(psb_spk_), intent(in) :: x,y,z b2=szero end function b2 function b3(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: b3 - real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_), intent(in) :: x,y,z b3=szero end function b3 function c(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: c - real(psb_spk_), intent(in) :: x,y,z + real(psb_spk_), intent(in) :: x,y,z c=szero end function c function a1(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none - real(psb_spk_) :: a1 + implicit none + real(psb_spk_) :: a1 real(psb_spk_), intent(in) :: x,y,z a1=sone/80 end function a1 function a2(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: a2 real(psb_spk_), intent(in) :: x,y,z a2=sone/80 end function a2 function a3(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: a3 real(psb_spk_), intent(in) :: x,y,z a3=sone/80 end function a3 function g(x,y,z) use psb_base_mod, only : psb_spk_, sone, szero - implicit none + implicit none real(psb_spk_) :: g real(psb_spk_), intent(in) :: x,y,z g = szero if (x == sone) then g = sone - else if (x == szero) then + else if (x == szero) then g = exp(y**2-z**2) end if end function g - + ! ! subroutine to allocate and fill in the coefficient matrix and - ! the rhs. + ! the rhs. ! subroutine psb_s_gen_pde3d(ictxt,idim,a,bv,xv,desc_a,afmt,info,& & f,amold,vmold,imold,partition,nrl,iv) @@ -174,13 +174,13 @@ contains use psb_util_mod ! ! Discretizes the partial differential equation - ! - ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) + ! + ! a1 dd(u) a2 dd(u) a3 dd(u) b1 d(u) b2 d(u) b3 d(u) ! - ------ - ------ - ------ + ----- + ------ + ------ + c u = f - ! dxdx dydy dzdz dx dy dz + ! dxdx dydy dzdz dx dy dz ! ! with Dirichlet boundary conditions - ! u = g + ! u = g ! ! on the unit cube 0<=x,y,z<=1. ! @@ -196,7 +196,7 @@ contains character(len=*) :: afmt procedure(s_func_3d), optional :: f class(psb_s_base_sparse_mat), optional :: amold - class(psb_s_base_vect_type), optional :: vmold + class(psb_s_base_vect_type), optional :: vmold class(psb_i_base_vect_type), optional :: imold integer(psb_ipk_), optional :: partition, nrl,iv(:) @@ -237,7 +237,7 @@ contains call psb_info(ictxt, iam, np) - if (present(f)) then + if (present(f)) then f_ => f else f_ => s_null_func_3d @@ -257,10 +257,10 @@ contains else partition_ = 3 end if - + ! initialize array descriptor and sparse matrix storage. provide an - ! estimate of the number of non zeroes - + ! estimate of the number of non zeroes + m = (1_psb_lpk_*idim)*idim*idim n = m nnz = ((n*7)/(np)) @@ -268,8 +268,8 @@ contains t0 = psb_wtime() select case(partition_) case(1) - ! A BLOCK partition - if (present(nrl)) then + ! A BLOCK partition + if (present(nrl)) then nr = nrl else ! @@ -280,46 +280,46 @@ contains end if nt = nr - call psb_sum(ictxt,nt) - if (nt /= m) then + call psb_sum(ictxt,nt) + if (nt /= m) then write(psb_err_unit,*) iam, 'Initialization error ',nr,nt,m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! First example of use of CDALL: specify for each process a number of ! contiguous rows - ! + ! call psb_cdall(ictxt,desc_a,info,nl=nr) myidx = desc_a%get_global_indices() nlr = size(myidx) case(2) ! A partition defined by the user through IV - - if (present(iv)) then + + if (present(iv)) then if (size(iv) /= m) then write(psb_err_unit,*) iam, 'Initialization error: wrong IV size',size(iv),m info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if else write(psb_err_unit,*) iam, 'Initialization error: IV not present' info = -1 call psb_barrier(ictxt) call psb_abort(ictxt) - return + return end if ! ! Second example of use of CDALL: specify for each row the - ! process that owns it - ! + ! process that owns it + ! call psb_cdall(ictxt,desc_a,info,vg=iv) myidx = desc_a%get_global_indices() nlr = size(myidx) @@ -335,7 +335,7 @@ contains npz = npdims(3) allocate(bndx(0:npx),bndy(0:npy),bndz(0:npz)) - ! We can reuse idx2ijk for process indices as well. + ! We can reuse idx2ijk for process indices as well. call idx2ijk(iamx,iamy,iamz,iam,npx,npy,npz,base=0) ! Now let's split the 3D cube in hexahedra call dist1Didx(bndx,idim,npx) @@ -345,7 +345,7 @@ contains call dist1Didx(bndz,idim,npz) mynz = bndz(iamz+1)-bndz(iamz) - ! How many indices do I own? + ! How many indices do I own? nlr = mynx*myny*mynz allocate(myidx(nlr)) ! Now, let's generate the list of indices I own @@ -369,9 +369,9 @@ contains ! ! Third example of use of CDALL: specify for each process ! the set of global indices it owns. - ! + ! call psb_cdall(ictxt,desc_a,info,vl=myidx) - + case default write(psb_err_unit,*) iam, 'Initialization error: should not get here' info = -1 @@ -380,9 +380,9 @@ contains return end select - + if (info == psb_success_) call psb_spall(a,desc_a,info,nnz=nnz) - ! define rhs from boundary conditions; also build initial guess + ! define rhs from boundary conditions; also build initial guess if (info == psb_success_) call psb_geall(xv,desc_a,info) if (info == psb_success_) call psb_geall(bv,desc_a,info) @@ -397,12 +397,12 @@ contains end if ! we build an auxiliary matrix consisting of one row at a - ! time; just a small matrix. might be extended to generate - ! a bunch of rows per call. - ! + ! time; just a small matrix. might be extended to generate + ! a bunch of rows per call. + ! allocate(val(20*nb),irow(20*nb),& &icol(20*nb),stat=info) - if (info /= psb_success_ ) then + if (info /= psb_success_ ) then info=psb_err_alloc_dealloc_ call psb_errpush(info,name) goto 9999 @@ -415,11 +415,11 @@ contains call psb_barrier(ictxt) t1 = psb_wtime() do ii=1, nlr,nb - ib = min(nb,nlr-ii+1) + ib = min(nb,nlr-ii+1) icoeff = 1 do k=1,ib i=ii+k-1 - ! local matrix pointer + ! local matrix pointer glob_row=myidx(i) ! compute gridpoint coordinates call idx2ijk(ix,iy,iz,glob_row,idim,idim,idim) @@ -429,11 +429,11 @@ contains z = (iz-1)*deltah zt(k) = f_(x,y,z) ! internal point: build discretization - ! + ! ! term depending on (x-1,y,z) ! val(icoeff) = -a1(x,y,z)/sqdeltah-b1(x,y,z)/deltah2 - if (ix == 1) then + if (ix == 1) then zt(k) = g(szero,y,z)*(-val(icoeff)) + zt(k) else call ijk2idx(icol(icoeff),ix-1,iy,iz,idim,idim,idim) @@ -442,19 +442,19 @@ contains endif ! term depending on (x,y-1,z) val(icoeff) = -a2(x,y,z)/sqdeltah-b2(x,y,z)/deltah2 - if (iy == 1) then + if (iy == 1) then zt(k) = g(x,szero,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy-1,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x,y,z-1) val(icoeff)=-a3(x,y,z)/sqdeltah-b3(x,y,z)/deltah2 - if (iz == 1) then + if (iz == 1) then zt(k) = g(x,y,szero)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz-1,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif @@ -462,33 +462,33 @@ contains ! term depending on (x,y,z) val(icoeff)=(2*sone)*(a1(x,y,z)+a2(x,y,z)+a3(x,y,z))/sqdeltah & & + c(x,y,z) - call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz,idim,idim,idim) irow(icoeff) = glob_row - icoeff = icoeff+1 + icoeff = icoeff+1 ! term depending on (x,y,z+1) val(icoeff)=-a3(x,y,z)/sqdeltah+b3(x,y,z)/deltah2 - if (iz == idim) then + if (iz == idim) then zt(k) = g(x,y,sone)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy,iz+1,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x,y+1,z) val(icoeff)=-a2(x,y,z)/sqdeltah+b2(x,y,z)/deltah2 - if (iy == idim) then + if (iy == idim) then zt(k) = g(x,sone,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix,iy+1,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif ! term depending on (x+1,y,z) val(icoeff)=-a1(x,y,z)/sqdeltah+b1(x,y,z)/deltah2 - if (ix==idim) then + if (ix==idim) then zt(k) = g(sone,y,z)*(-val(icoeff)) + zt(k) else - call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) + call ijk2idx(icol(icoeff),ix+1,iy,iz,idim,idim,idim) irow(icoeff) = glob_row icoeff = icoeff+1 endif @@ -519,8 +519,8 @@ contains tcdasb = psb_wtime()-t1 call psb_barrier(ictxt) t1 = psb_wtime() - if (info == psb_success_) then - if (present(amold)) then + if (info == psb_success_) then + if (present(amold)) then call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,mold=amold) else call psb_spasb(a,desc_a,info,dupl=psb_dupl_err_,afmt=afmt) @@ -543,7 +543,7 @@ contains end if tasb = psb_wtime()-t1 call psb_barrier(ictxt) - ttot = psb_wtime() - t0 + ttot = psb_wtime() - t0 call psb_amx(ictxt,talc) call psb_amx(ictxt,tgen) @@ -585,9 +585,9 @@ program psb_s_pde3d integer(psb_ipk_) :: idim integer(psb_epk_) :: system_size - ! miscellaneous + ! miscellaneous real(psb_spk_), parameter :: one = sone - real(psb_dpk_) :: t1, t2, tprec + real(psb_dpk_) :: t1, t2, tprec ! sparse matrix and preconditioner type(psb_sspmat_type) :: a @@ -604,6 +604,14 @@ program psb_s_pde3d integer(psb_epk_) :: amatsize, precsize, descsize, d2size real(psb_spk_) :: err, eps + ! Parameters for solvers in Block-Jacobi preconditioner + type ainvparms + character(len=12) :: alg, orth_alg + integer(psb_ipk_) :: fill, inv_fill + real(psb_dpk_) :: thresh, inv_thresh + end type ainvparms + type(ainvparms) :: parms + ! other variables integer(psb_ipk_) :: info, i character(len=20) :: name,ch_err @@ -615,7 +623,7 @@ program psb_s_pde3d call psb_init(ictxt) call psb_info(ictxt,iam,np) - if (iam < 0) then + if (iam < 0) then ! This should not happen, but just in case call psb_exit(ictxt) stop @@ -626,21 +634,21 @@ program psb_s_pde3d ! ! Hello world ! - if (iam == psb_root_) then + if (iam == psb_root_) then write(*,*) 'Welcome to PSBLAS version: ',psb_version_string_ write(*,*) 'This is the ',trim(name),' sample program' end if ! ! get parameters ! - call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + call get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart,parms) ! - ! allocate and fill in the coefficient matrix, rhs and initial guess + ! allocate and fill in the coefficient matrix, rhs and initial guess ! call psb_barrier(ictxt) t1 = psb_wtime() - call psb_gen_pde3d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) + call psb_gen_pde3d(ictxt,idim,a,bv,xxv,desc_a,afmt,info,partition=ipart) call psb_barrier(ictxt) t2 = psb_wtime() - t1 if(info /= psb_success_) then @@ -653,7 +661,7 @@ program psb_s_pde3d if (iam == psb_root_) write(psb_out_unit,'(" ")') ! ! prepare the preconditioner. - ! + ! if(iam == psb_root_) write(psb_out_unit,'("Setting preconditioner to : ",a)')ptype call prec%init(ictxt,ptype,info) @@ -675,14 +683,14 @@ program psb_s_pde3d if (iam == psb_root_) write(psb_out_unit,'(" ")') call prec%descr() ! - ! iterative method parameters + ! iterative method parameters ! if(iam == psb_root_) write(psb_out_unit,'("Calling iterative method ",a)')kmethd call psb_barrier(ictxt) - t1 = psb_wtime() + t1 = psb_wtime() eps = 1.d-6 - call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& - & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) + call psb_krylov(kmethd,a,prec,bv,xxv,eps,desc_a,info,& + & itmax=itmax,iter=iter,err=err,itrace=itrace,istop=istopc,irst=irst) if(info /= psb_success_) then info=psb_err_from_subroutine_ @@ -712,14 +720,14 @@ program psb_s_pde3d write(psb_out_unit,'("Convergence indicator on exit : ",es12.5)')err write(psb_out_unit,'("Info on exit : ",i12)')info write(psb_out_unit,'("Total memory occupation for A: ",i12)')amatsize - write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize + write(psb_out_unit,'("Total memory occupation for PREC: ",i12)')precsize write(psb_out_unit,'("Total memory occupation for DESC_A: ",i12)')descsize write(psb_out_unit,'("Storage format for A: ",a)') a%get_fmt() write(psb_out_unit,'("Storage format for DESC_A: ",a)') desc_a%get_fmt() end if - ! + ! ! cleanup storage and exit ! call psb_gefree(bv,desc_a,info) @@ -745,13 +753,15 @@ contains ! ! get iteration parameters from standard input ! - subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,itmax,itrace,irst,ipart) + subroutine get_parms(ictxt,kmethd,ptype,afmt,idim,istopc,& + itmax,itrace,irst,ipart,parms) integer(psb_ipk_) :: ictxt character(len=*) :: kmethd, ptype, afmt integer(psb_ipk_) :: idim, istopc,itmax,itrace,irst,ipart integer(psb_ipk_) :: np, iam integer(psb_ipk_) :: ip, inp_unit character(len=1024) :: filename + type(ainvparms) :: parms call psb_info(ictxt, iam, np) @@ -780,12 +790,12 @@ contains if (ip >= 4) then read(inp_unit,*) ipart else - ipart = 3 + ipart = 3 endif if (ip >= 5) then read(inp_unit,*) istopc else - istopc=1 + istopc=1 endif if (ip >= 6) then read(inp_unit,*) itmax @@ -802,8 +812,23 @@ contains else irst=1 endif + if (ip >= 9) then + read(psb_inp_unit,*) parms%alg + read(psb_inp_unit,*) parms%fill + read(psb_inp_unit,*) parms%inv_fill + read(psb_inp_unit,*) parms%thresh + read(psb_inp_unit,*) parms%inv_thresh + read(psb_inp_unit,*) parms%orth_alg + else + parms%alg = 'ILU' ! AINV variant: ILU,ILUT,MILU,INVK,AINVT,AORTH + parms%fill = 0 ! Fill in for forward factorization + parms%inv_fill = 1 ! Fill in for inverse factorization + parms%thresh = 1E-1_psb_spk_ ! Threshold for forward factorization + parms%inv_thresh = 1E-1_psb_spk_ ! Threshold for inverse factorization + parms%orth_alg = 'LLK' ! What orthogonalization algorithm? + endif - write(psb_out_unit,'("Solving matrix : ell1")') + write(psb_out_unit,'("Solving matrix : ell1")') write(psb_out_unit,& & '("Grid dimensions : ",i4," x ",i4," x ",i4)') & & idim,idim,idim @@ -818,11 +843,28 @@ contains write(psb_out_unit,'("Unknown data distrbution, defaulting to 3D")') end select write(psb_out_unit,'("Preconditioner : ",a)') ptype + if( psb_toupper(ptype) == "BJAC" ) then + write(psb_out_unit,'("Block subsolver : ",a)') parms%alg + select case (psb_toupper(parms%alg)) + case ('ILU','ILUT','MILU') + write(psb_out_unit,'("Fill in : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + case ('INVK') + write(psb_out_unit,'("Fill : ",i0)') parms%fill + write(psb_out_unit,'("Threshold : ",e2.2)') parms%thresh + write(psb_out_unit,'("Invese Fill : ",i0)') parms%inv_fill + write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + case ('AINVT','AORTH') + write(psb_out_unit,'("Inverse Threshold : ",e2.2)') parms%inv_thresh + case default + write(psb_out_unit,'("Unknown diagonal solver")') + end select + end if write(psb_out_unit,'("Iterative method : ",a)') kmethd write(psb_out_unit,'(" ")') else ! wrong number of parameter, print an error message and exit - call pr_usage(izero) + call pr_usage(izero) call psb_abort(ictxt) stop 1 endif @@ -841,20 +883,26 @@ contains call psb_bcast(ictxt,itmax) call psb_bcast(ictxt,itrace) call psb_bcast(ictxt,irst) + call psb_bcast(ictxt,parms%alg) + call psb_bcast(ictxt,parms%fill) + call psb_bcast(ictxt,parms%inv_fill) + call psb_bcast(ictxt,parms%thresh) + call psb_bcast(ictxt,parms%inv_thresh) + call psb_bcast(ictxt,parms%orth_alg) return end subroutine get_parms ! - ! print an error message - ! + ! print an error message + ! subroutine pr_usage(iout) integer(psb_ipk_) :: iout write(iout,*)'incorrect parameter(s) found' write(iout,*)' usage: pde3d90 methd prec dim & - &[istop itmax itrace]' + &[istop itmax itrace]' write(iout,*)' where:' - write(iout,*)' methd: cgstab cgs rgmres bicgstabl' + write(iout,*)' methd: cgstab cgs rgmres bicgstabl' write(iout,*)' prec : bjac diag none' write(iout,*)' dim number of points along each axis' write(iout,*)' the size of the resulting linear ' @@ -862,11 +910,9 @@ contains write(iout,*)' ipart data partition 1 3 ' write(iout,*)' istop stopping criterion 1, 2 ' write(iout,*)' itmax maximum number of iterations [500] ' - write(iout,*)' itrace <=0 (no tracing, default) or ' + write(iout,*)' itrace <=0 (no tracing, default) or ' write(iout,*)' >= 1 do tracing every itrace' - write(iout,*)' iterations ' + write(iout,*)' iterations ' end subroutine pr_usage end program psb_s_pde3d - - diff --git a/test/pargen/runs/ppde.inp b/test/pargen/runs/ppde.inp index f6fe33eb..d54ce5b6 100644 --- a/test/pargen/runs/ppde.inp +++ b/test/pargen/runs/ppde.inp @@ -8,5 +8,11 @@ CSR Storage format for matrix A: CSR COO 0100 MAXIT 05 ITRACE 002 IRST restart for RGMRES and BiCGSTABL +ILU Factorization variant: ILU,ILUT,MILU,INVK,AINVT,AORTH +0 Fill in for forward factorization +1 Fill in for inverse factorization (ignored if not INVK) +1E-1 Threshold for forward factorization (ignored if ILU) +1E-1 Threshold for inverse factorization (ignored if ILU,ILUT,MILU) +LLK What orthogonalization algorithm? (ignored if ILU,ILUT,MILU,INVK)