fdurastante/amg4psblas
1
0
Fork 0
mirror of https://github.com/sfilippone/amg4psblas.git synced 2026-10-16 19:02:10 +00:00
Code Issues Packages Projects Releases Wiki Activity

Changes for new L1-Jacobi smoother/solver.

Browse Source
This commit is contained in:
Salvatore Filippone
2019-05-07 17:24:01 +01:00
parent b2047f95e0
commit 88976ec108
33 changed files with 1000 additions and 62 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/impl/solver/mld_c_diag_solver_bld.f90
+78
View File
@@ -111,3 +111,81 @@ subroutine mld_c_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
return
end subroutine mld_c_diag_solver_bld
subroutine mld_c_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_c_l1_diag_solver, mld_protect_name => mld_c_l1_diag_solver_bld
Implicit None
! Arguments
type(psb_cspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(inout) :: desc_a
class(mld_c_l1_diag_solver_type), intent(inout) :: sv
integer(psb_ipk_), intent(out) :: info
type(psb_cspmat_type), intent(in), target, optional :: b
class(psb_c_base_sparse_mat), intent(in), optional :: amold
class(psb_c_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
complex(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_spk_), allocatable :: tdb(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='c_l1_diag_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
nrow_a = a%get_nrows()
sv%d = a%arwsum(info)
if (info == psb_success_) call psb_realloc(n_row,sv%d,info)
if (present(b)) then
tdb=b%arwsum(info)
if (size(tdb)+nrow_a > n_row) call psb_realloc(nrow_a+size(tdb),sv%d,info)
if (info == psb_success_) sv%d(nrow_a+1:nrow_a+size(tdb)) = tdb(:)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='arwsum')
goto 9999
end if
do i=1,n_row
if (sv%d(i) == czero) then
sv%d(i) = cone
else
sv%d(i) = cone/sv%d(i)
end if
end do
allocate(sv%dv,stat=info)
if (info == psb_success_) then
call sv%dv%bld(sv%d)
if (present(vmold)) call sv%dv%cnv(vmold)
call sv%dv%sync()
else
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='Allocate sv%dv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_l1_diag_solver_bld
mlprec/impl/solver/mld_c_diag_solver_clone.f90
+2 -2
View File
@@ -57,14 +57,14 @@ subroutine mld_c_diag_solver_clone(sv,svout,info)
if (info == psb_success_) deallocate(svout, stat=info)
end if
if (info == psb_success_) &
& allocate(mld_c_diag_solver_type :: svout, stat=info)
& allocate(svout, mold=sv, stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
goto 9999
end if
select type(svo => svout)
type is (mld_c_diag_solver_type)
class is (mld_c_diag_solver_type)
call psb_safe_ab_cpy(sv%d,svo%d,info)
if (info == psb_success_) &
& call sv%dv%clone(svo%dv,info)
mlprec/impl/solver/mld_c_diag_solver_dmp.f90
+46
View File
@@ -81,3 +81,49 @@ subroutine mld_c_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
end if
end subroutine mld_c_diag_solver_dmp
subroutine mld_c_l1_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_base_mod
use mld_c_l1_diag_solver, mld_protect_name => mld_c_l1_diag_solver_dmp
implicit none
class(mld_c_l1_diag_solver_type), intent(in) :: sv
integer(psb_ipk_), intent(in) :: ictxt,level
integer(psb_ipk_), intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
if (solver_) then
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_c"
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_l1_diag.mtx'
if (allocated(sv%d)) &
& call psb_geprt(fname,sv%d,head=head)
end if
end subroutine mld_c_l1_diag_solver_dmp
mlprec/impl/solver/mld_d_diag_solver_bld.f90
+78
View File
@@ -111,3 +111,81 @@ subroutine mld_d_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
return
end subroutine mld_d_diag_solver_bld
subroutine mld_d_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_d_l1_diag_solver, mld_protect_name => mld_d_l1_diag_solver_bld
Implicit None
! Arguments
type(psb_dspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(inout) :: desc_a
class(mld_d_l1_diag_solver_type), intent(inout) :: sv
integer(psb_ipk_), intent(out) :: info
type(psb_dspmat_type), intent(in), target, optional :: b
class(psb_d_base_sparse_mat), intent(in), optional :: amold
class(psb_d_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
real(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_dpk_), allocatable :: tdb(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='d_l1_diag_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
nrow_a = a%get_nrows()
sv%d = a%arwsum(info)
if (info == psb_success_) call psb_realloc(n_row,sv%d,info)
if (present(b)) then
tdb=b%arwsum(info)
if (size(tdb)+nrow_a > n_row) call psb_realloc(nrow_a+size(tdb),sv%d,info)
if (info == psb_success_) sv%d(nrow_a+1:nrow_a+size(tdb)) = tdb(:)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='arwsum')
goto 9999
end if
do i=1,n_row
if (sv%d(i) == dzero) then
sv%d(i) = done
else
sv%d(i) = done/sv%d(i)
end if
end do
allocate(sv%dv,stat=info)
if (info == psb_success_) then
call sv%dv%bld(sv%d)
if (present(vmold)) call sv%dv%cnv(vmold)
call sv%dv%sync()
else
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='Allocate sv%dv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_l1_diag_solver_bld
mlprec/impl/solver/mld_d_diag_solver_clone.f90
+2 -2
View File
@@ -57,14 +57,14 @@ subroutine mld_d_diag_solver_clone(sv,svout,info)
if (info == psb_success_) deallocate(svout, stat=info)
end if
if (info == psb_success_) &
& allocate(mld_d_diag_solver_type :: svout, stat=info)
& allocate(svout, mold=sv, stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
goto 9999
end if
select type(svo => svout)
type is (mld_d_diag_solver_type)
class is (mld_d_diag_solver_type)
call psb_safe_ab_cpy(sv%d,svo%d,info)
if (info == psb_success_) &
& call sv%dv%clone(svo%dv,info)
mlprec/impl/solver/mld_d_diag_solver_dmp.f90
+46
View File
@@ -81,3 +81,49 @@ subroutine mld_d_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
end if
end subroutine mld_d_diag_solver_dmp
subroutine mld_d_l1_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_base_mod
use mld_d_l1_diag_solver, mld_protect_name => mld_d_l1_diag_solver_dmp
implicit none
class(mld_d_l1_diag_solver_type), intent(in) :: sv
integer(psb_ipk_), intent(in) :: ictxt,level
integer(psb_ipk_), intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
if (solver_) then
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_d"
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_l1_diag.mtx'
if (allocated(sv%d)) &
& call psb_geprt(fname,sv%d,head=head)
end if
end subroutine mld_d_l1_diag_solver_dmp
mlprec/impl/solver/mld_s_diag_solver_bld.f90
+78
View File
@@ -111,3 +111,81 @@ subroutine mld_s_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
return
end subroutine mld_s_diag_solver_bld
subroutine mld_s_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_s_l1_diag_solver, mld_protect_name => mld_s_l1_diag_solver_bld
Implicit None
! Arguments
type(psb_sspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(inout) :: desc_a
class(mld_s_l1_diag_solver_type), intent(inout) :: sv
integer(psb_ipk_), intent(out) :: info
type(psb_sspmat_type), intent(in), target, optional :: b
class(psb_s_base_sparse_mat), intent(in), optional :: amold
class(psb_s_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
real(psb_spk_), pointer :: ww(:), aux(:), tx(:),ty(:)
real(psb_spk_), allocatable :: tdb(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='s_l1_diag_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
nrow_a = a%get_nrows()
sv%d = a%arwsum(info)
if (info == psb_success_) call psb_realloc(n_row,sv%d,info)
if (present(b)) then
tdb=b%arwsum(info)
if (size(tdb)+nrow_a > n_row) call psb_realloc(nrow_a+size(tdb),sv%d,info)
if (info == psb_success_) sv%d(nrow_a+1:nrow_a+size(tdb)) = tdb(:)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='arwsum')
goto 9999
end if
do i=1,n_row
if (sv%d(i) == szero) then
sv%d(i) = sone
else
sv%d(i) = sone/sv%d(i)
end if
end do
allocate(sv%dv,stat=info)
if (info == psb_success_) then
call sv%dv%bld(sv%d)
if (present(vmold)) call sv%dv%cnv(vmold)
call sv%dv%sync()
else
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='Allocate sv%dv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_l1_diag_solver_bld
mlprec/impl/solver/mld_s_diag_solver_clone.f90
+2 -2
View File
@@ -57,14 +57,14 @@ subroutine mld_s_diag_solver_clone(sv,svout,info)
if (info == psb_success_) deallocate(svout, stat=info)
end if
if (info == psb_success_) &
& allocate(mld_s_diag_solver_type :: svout, stat=info)
& allocate(svout, mold=sv, stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
goto 9999
end if
select type(svo => svout)
type is (mld_s_diag_solver_type)
class is (mld_s_diag_solver_type)
call psb_safe_ab_cpy(sv%d,svo%d,info)
if (info == psb_success_) &
& call sv%dv%clone(svo%dv,info)
mlprec/impl/solver/mld_s_diag_solver_dmp.f90
+46
View File
@@ -81,3 +81,49 @@ subroutine mld_s_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
end if
end subroutine mld_s_diag_solver_dmp
subroutine mld_s_l1_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_base_mod
use mld_s_l1_diag_solver, mld_protect_name => mld_s_l1_diag_solver_dmp
implicit none
class(mld_s_l1_diag_solver_type), intent(in) :: sv
integer(psb_ipk_), intent(in) :: ictxt,level
integer(psb_ipk_), intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
if (solver_) then
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_s"
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_l1_diag.mtx'
if (allocated(sv%d)) &
& call psb_geprt(fname,sv%d,head=head)
end if
end subroutine mld_s_l1_diag_solver_dmp
mlprec/impl/solver/mld_z_diag_solver_bld.f90
+78
View File
@@ -111,3 +111,81 @@ subroutine mld_z_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
return
end subroutine mld_z_diag_solver_bld
subroutine mld_z_l1_diag_solver_bld(a,desc_a,sv,info,b,amold,vmold,imold)
use psb_base_mod
use mld_z_l1_diag_solver, mld_protect_name => mld_z_l1_diag_solver_bld
Implicit None
! Arguments
type(psb_zspmat_type), intent(in), target :: a
Type(psb_desc_type), Intent(inout) :: desc_a
class(mld_z_l1_diag_solver_type), intent(inout) :: sv
integer(psb_ipk_), intent(out) :: info
type(psb_zspmat_type), intent(in), target, optional :: b
class(psb_z_base_sparse_mat), intent(in), optional :: amold
class(psb_z_base_vect_type), intent(in), optional :: vmold
class(psb_i_base_vect_type), intent(in), optional :: imold
! Local variables
integer(psb_ipk_) :: n_row,n_col, nrow_a, nztota
complex(psb_dpk_), pointer :: ww(:), aux(:), tx(:),ty(:)
complex(psb_dpk_), allocatable :: tdb(:)
integer(psb_ipk_) :: ictxt,np,me,i, err_act, debug_unit, debug_level
character(len=20) :: name='z_l1_diag_solver_bld', ch_err
info=psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
ictxt = desc_a%get_context()
call psb_info(ictxt, me, np)
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' start'
n_row = desc_a%get_local_rows()
nrow_a = a%get_nrows()
sv%d = a%arwsum(info)
if (info == psb_success_) call psb_realloc(n_row,sv%d,info)
if (present(b)) then
tdb=b%arwsum(info)
if (size(tdb)+nrow_a > n_row) call psb_realloc(nrow_a+size(tdb),sv%d,info)
if (info == psb_success_) sv%d(nrow_a+1:nrow_a+size(tdb)) = tdb(:)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='arwsum')
goto 9999
end if
do i=1,n_row
if (sv%d(i) == zzero) then
sv%d(i) = zone
else
sv%d(i) = zone/sv%d(i)
end if
end do
allocate(sv%dv,stat=info)
if (info == psb_success_) then
call sv%dv%bld(sv%d)
if (present(vmold)) call sv%dv%cnv(vmold)
call sv%dv%sync()
else
call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='Allocate sv%dv')
goto 9999
end if
if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),' end'
call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_l1_diag_solver_bld
mlprec/impl/solver/mld_z_diag_solver_clone.f90
+2 -2
View File
@@ -57,14 +57,14 @@ subroutine mld_z_diag_solver_clone(sv,svout,info)
if (info == psb_success_) deallocate(svout, stat=info)
end if
if (info == psb_success_) &
& allocate(mld_z_diag_solver_type :: svout, stat=info)
& allocate(svout, mold=sv, stat=info)
if (info /= 0) then
info = psb_err_alloc_dealloc_
goto 9999
end if
select type(svo => svout)
type is (mld_z_diag_solver_type)
class is (mld_z_diag_solver_type)
call psb_safe_ab_cpy(sv%d,svo%d,info)
if (info == psb_success_) &
& call sv%dv%clone(svo%dv,info)
mlprec/impl/solver/mld_z_diag_solver_dmp.f90
+46
View File
@@ -81,3 +81,49 @@ subroutine mld_z_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
end if
end subroutine mld_z_diag_solver_dmp
subroutine mld_z_l1_diag_solver_dmp(sv,ictxt,level,info,prefix,head,solver)
use psb_base_mod
use mld_z_l1_diag_solver, mld_protect_name => mld_z_l1_diag_solver_dmp
implicit none
class(mld_z_l1_diag_solver_type), intent(in) :: sv
integer(psb_ipk_), intent(in) :: ictxt,level
integer(psb_ipk_), intent(out) :: info
character(len=*), intent(in), optional :: prefix, head
logical, optional, intent(in) :: solver
integer(psb_ipk_) :: i, j, il1, iln, lname, lev
integer(psb_ipk_) :: icontxt,iam, np
character(len=80) :: prefix_
character(len=120) :: fname ! len should be at least 20 more than
logical :: solver_
! len of prefix_
info = 0
call psb_info(ictxt,iam,np)
if (present(solver)) then
solver_ = solver
else
solver_ = .false.
end if
if (solver_) then
if (present(prefix)) then
prefix_ = trim(prefix(1:min(len(prefix),len(prefix_))))
else
prefix_ = "dump_slv_z"
end if
lname = len_trim(prefix_)
fname = trim(prefix_)
write(fname(lname+1:lname+5),'(a,i3.3)') '_p',iam
lname = lname + 5
write(fname(lname+1:),'(a,i3.3,a)')'_l',level,'_l1_diag.mtx'
if (allocated(sv%d)) &
& call psb_geprt(fname,sv%d,head=head)
end if
end subroutine mld_z_l1_diag_solver_dmp