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

Added implementation in BJAC and test for ILU-type factorizations

Browse Source
implement-ainv
Cirdans-Home 4 years ago
parent f0bb949192
commit fbf23c3959
20 changed files with 2230 additions and 818 deletions
Show Stats Download Patch File Download Diff File

14
base/modules/serial/psb_c_mat_mod.F90
Unescape Escape View File

@ -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

14
base/modules/serial/psb_d_mat_mod.F90
Unescape Escape View File

@ -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

14
base/modules/serial/psb_s_mat_mod.F90
Unescape Escape View File

@ -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

14
base/modules/serial/psb_z_mat_mod.F90
Unescape Escape View File

@ -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

43
base/serial/impl/psb_c_mat_impl.F90
Unescape Escape View File

@ -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)

43
base/serial/impl/psb_d_mat_impl.F90
Unescape Escape View File

@ -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)

43
base/serial/impl/psb_s_mat_impl.F90
Unescape Escape View File

@ -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)

43
base/serial/impl/psb_z_mat_impl.F90
Unescape Escape View File

@ -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)

425
prec/impl/psb_c_bjacprec_impl.f90
Unescape Escape View File

@ -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

425
prec/impl/psb_d_bjacprec_impl.f90
Unescape Escape View File

@ -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

425
prec/impl/psb_s_bjacprec_impl.f90
Unescape Escape View File

@ -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

425
prec/impl/psb_z_bjacprec_impl.f90
Unescape Escape View File

@ -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

134
prec/psb_c_bjacprec.f90
Unescape Escape View File

@ -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

134
prec/psb_d_bjacprec.f90
Unescape Escape View File

@ -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

46
prec/psb_prec_const_mod.f90
Unescape Escape View File

@ -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

134
prec/psb_s_bjacprec.f90
Unescape Escape View File

@ -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

134
prec/psb_z_bjacprec.f90
Unescape Escape View File

@ -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

266
test/pargen/psb_d_pde3d.f90
Unescape Escape View File

@ -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

266
test/pargen/psb_s_pde3d.f90
Unescape Escape View File

@ -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

6
test/pargen/runs/ppde.inp
Unescape Escape View File

@ -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)

Write Preview
Loading…
Cancel
Save