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

mld2p4-2:

Browse Source 
mlprec/impl/level/mld_c_base_onelev_build.f90
 mlprec/impl/level/mld_d_base_onelev_build.f90
 mlprec/impl/level/mld_s_base_onelev_build.f90
 mlprec/impl/level/mld_z_base_onelev_build.f90
 mlprec/impl/mld_cmlprec_aply.f90
 mlprec/impl/mld_dmlprec_aply.f90
 mlprec/impl/mld_smlprec_aply.f90
 mlprec/impl/mld_zmlprec_aply.f90
 mlprec/mld_base_prec_type.F90

First major tep in restructuring MLPREC_APLY.
stopcriterion
Salvatore Filippone 9 years ago
parent a9b9ea958d
commit 436f3e49f7
5 changed files with 2684 additions and 2717 deletions
Show Stats Download Patch File Download Diff File

647
mlprec/impl/mld_cmlprec_aply.f90
Unescape Escape View File

@ -39,7 +39,7 @@
! File: mld_cmlprec_aply.f90 ! File: mld_cmlprec_aply.f90
! !
! Subroutine: mld_cmlprec_aply ! Subroutine: mld_cmlprec_aply
! Version: complex ! Version: real
! !
! This routine computes ! This routine computes
! !
@ -84,7 +84,7 @@
! The multilevel preconditioner data structure containing the ! The multilevel preconditioner data structure containing the
! local part of the preconditioner to be applied. ! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! p%precv(ilev)%prec - type(psb_dbaseprec_type) ! p%precv(ilev)%prec - type(psb_cbaseprec_type)
! The 'base preconditioner' for the current level ! The 'base preconditioner' for the current level
! p%precv(ilev)%ac - type(psb_cspmat_type) ! p%precv(ilev)%ac - type(psb_cspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
@ -333,7 +333,7 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_cmlprec_aply' name='mld_cmlprec_aply'
info = psb_success_ info = psb_success_
@ -1016,11 +1016,12 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act
character(len=20) :: name character(len=20) :: name
character :: trans_ character :: trans_
complex(psb_spk_) :: par
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_cmlprec_aply' name='mld_cmlprec_aply'
info = psb_success_ info = psb_success_
@ -1036,7 +1037,6 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1065,12 +1065,12 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999 goto 9999
end if end if
end do end do
level = 1 level = 1
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call mlprec_wrk(level)%vy2l%zero() call mlprec_wrk(level)%vy2l%zero()
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1082,6 +1082,7 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
do level = 1, nlev do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info) call mlprec_wrk(level)%vtx%free(info)
@ -1118,29 +1119,30 @@ contains
! Arguments ! Arguments
integer(psb_ipk_) :: level integer(psb_ipk_) :: level
type(mld_cprec_type), target, intent(inout) :: p type(mld_cprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans character, intent(in) :: trans
complex(psb_spk_),target :: work(:) complex(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_c_vect_type),intent(inout), optional :: u type(psb_c_vect_type),intent(inout), optional :: u
type(psb_c_vect_type) :: res type(psb_c_vect_type) :: res
type(psb_c_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables ! Local variables
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name character(len=20) :: name
name = 'inner_ml_aply' name = 'inner_ml_aply'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1167,57 +1169,9 @@ contains
goto 9999 goto 9999
case(mld_add_ml_) case(mld_add_ml_)
!
! Additive multilevel
!
if (level > 1) then call mld_c_inner_add(p, mlprec_wrk, level, trans, work)
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,&
& czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply')
goto 9999
end if
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -1230,164 +1184,108 @@ contains
select case(p%precv(level)%parms%smoother_pos) select case(p%precv(level)%parms%smoother_pos)
case(mld_post_smooth_) case(mld_post_smooth_)
p%precv(level)%parms%sweeps_pre = 0
call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
select case (trans_)
case('N')
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,&
& czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
! This is one step of post-smoothing case(mld_pre_smooth_)
p%precv(level)%parms%sweeps_post = 0
call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
if (level < nlev) then case(mld_twoside_smooth_)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
! case default
! Apply the prolongator info = psb_err_from_subroutine_ai_
! call psb_errpush(info,name,a_err='invalid smooth_pos',&
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/))
& czero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999 goto 9999
end if
! end select
! Compute the residual
!
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
sweeps = p%precv(level)%parms%sweeps_post case(mld_mult_dev_ml_)
call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C') case(mld_vcycle_ml_, mld_wcycle_ml_)
! Post-smoothing transpose is pre-smoothing call mld_c_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
if (level > 1) then call mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u)
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,&
& czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then case default
call psb_errpush(psb_err_internal_error_,name,& info = psb_err_from_subroutine_ai_
& a_err='Error during restriction') call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999 goto 9999
end if
end select
call psb_erractionrestore(err_act)
return
end if 9999 call psb_error_handler(err_act)
return
! end subroutine inner_ml_aply
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
! recursive subroutine mld_c_inner_add(p, mlprec_wrk, level, trans, work)
! Compute the residual (at all levels but the coarsest one) use psb_base_mod
! use mld_prec_mod
if (level < nlev) then
call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) implicit none
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& cone,mlprec_wrk(level)%vy2l,& integer(psb_ipk_), intent(in) :: level
& p%precv(level+1)%map,info,work=work) character, intent(in) :: trans
if (info /= psb_success_) then complex(psb_spk_),target :: work(:)
type(psb_c_vect_type) :: res
type(psb_c_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='wrong call level to inner_add')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if end if
case default if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end select end if
case(mld_pre_smooth_)
select case (trans_)
case('N')
! One step of pre-smoothing
if (level > 1) then if (level > 1) then
@ -1404,78 +1302,17 @@ contains
end if end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps sweeps = p%precv(level)%parms%sweeps
end if call p%precv(level)%sm%apply(cone,&
call p%precv(level)%sm%apply(cone,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,&
& p%precv(level)%base_desc, trans,& & sweeps,work,info)
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case('T','C')
! pre-smooth transpose is post-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,&
& czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during ADD smoother_apply')
goto 9999 goto 9999
end if end if
end if
if (level < nlev) then if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1488,7 +1325,7 @@ contains
! Apply the prolongator ! Apply the prolongator
! !
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& czero,mlprec_wrk(level)%vy2l,& & cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1496,56 +1333,86 @@ contains
goto 9999 goto 9999
end if end if
!
! Compute the residual
!
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if end if
call psb_erractionrestore(err_act)
return
sweeps = p%precv(level)%parms%sweeps_pre 9999 call psb_error_handler(err_act)
call p%precv(level)%sm%apply(cone,& return
& mlprec_wrk(level)%vx2l,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& end subroutine mld_c_inner_add
& sweeps,work,info)
if (info /= psb_success_) then recursive subroutine mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
call psb_errpush(psb_err_internal_error_,name,& use psb_base_mod
& a_err='Error during PRE smoother_apply') use mld_prec_mod
goto 9999
end if implicit none
else
sweeps = p%precv(level)%parms%sweeps !Input/Oputput variables
call p%precv(level)%sm%apply(cone,& type(mld_cprec_type), intent(inout) :: p
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& sweeps,work,info) integer(psb_ipk_), intent(in) :: level
if (info /= psb_success_) then character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
type(psb_c_vect_type) :: res
type(psb_c_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_mult'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during smoother_apply') & a_err='wrong call level to inner_mult')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level
end if end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_twoside_smooth_) if ((level < nlev).or.(nlev == 1)) then
sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre
else
sweeps_post = p%precv(level-1)%parms%sweeps_post
sweeps_pre = p%precv(level-1)%parms%sweeps_pre
endif
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
if (level > 1) then if (level > 1) then
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vty,& if (pre) then
current => mlprec_wrk(level-1)%vty
else
current => mlprec_wrk(level-1)%vx2l
endif
call psb_map_X2Y(cone,current,&
& czero,mlprec_wrk(level)%vx2l,& & czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work) & p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -1553,13 +1420,14 @@ contains
end if end if
end if end if
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vtx,& if (level < nlev) then
& p%precv(level)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (level < nlev) then
if (pre) then
if (trans == 'N') then if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
@ -1573,26 +1441,18 @@ contains
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info)
end if end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply') & a_err='Error during 2-PRE smoother_apply')
goto 9999 goto 9999
end if end if
endif
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual and call recursively
! and call recursively
! !
if(level < nlev) then if (pre) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,& & czero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
@ -1605,6 +1465,7 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1617,10 +1478,16 @@ contains
! !
! Apply the prolongator ! Apply the prolongator
! !
if (pre) then
par = cone
else
par = czero
endif
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,& & par,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -1630,6 +1497,10 @@ contains
! !
! Compute the residual ! Compute the residual
! !
if (post) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& & cone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
@ -1660,19 +1531,88 @@ contains
& a_err='Error during 2-POST smoother_apply') & a_err='Error during 2-POST smoother_apply')
goto 9999 goto 9999
end if end if
endif endif
case default else if (level == nlev) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid smooth_pos',& sweeps = p%precv(level)%parms%sweeps
& i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end if
end select call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_inner_mult
recursive subroutine mld_c_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
type(psb_c_vect_type),intent(inout), optional :: u
type(psb_c_vect_type) :: res
type(psb_c_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle !V/W cycle
if (level > 1) then if (level > 1) then
@ -1703,7 +1643,7 @@ contains
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(czero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1766,7 +1706,7 @@ contains
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, u=mlprec_wrk(level+1)%vy2l)
endif endif
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1826,8 +1766,69 @@ contains
endif endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_c_inner_vw_cycle
recursive subroutine mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
type(psb_c_vect_type),intent(inout), optional :: u
type(psb_c_vect_type) :: res
type(psb_c_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
!K cycle !K cycle
@ -1844,7 +1845,7 @@ contains
& czero,mlprec_wrk(level)%vy2l,& & czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(czero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1984,23 +1985,13 @@ contains
endif endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine inner_ml_aply end subroutine mld_c_inner_k_cycle
recursive subroutine mld_cinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) recursive subroutine mld_cinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)

643
mlprec/impl/mld_dmlprec_aply.f90
Unescape Escape View File

@ -333,7 +333,7 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_dmlprec_aply' name='mld_dmlprec_aply'
info = psb_success_ info = psb_success_
@ -1016,11 +1016,12 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act
character(len=20) :: name character(len=20) :: name
character :: trans_ character :: trans_
real(psb_dpk_) :: par
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_dmlprec_aply' name='mld_dmlprec_aply'
info = psb_success_ info = psb_success_
@ -1036,7 +1037,6 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1065,12 +1065,12 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999 goto 9999
end if end if
end do end do
level = 1 level = 1
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call mlprec_wrk(level)%vy2l%zero() call mlprec_wrk(level)%vy2l%zero()
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1082,6 +1082,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
do level = 1, nlev do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info) call mlprec_wrk(level)%vtx%free(info)
@ -1118,29 +1119,30 @@ contains
! Arguments ! Arguments
integer(psb_ipk_) :: level integer(psb_ipk_) :: level
type(mld_dprec_type), target, intent(inout) :: p type(mld_dprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans character, intent(in) :: trans
real(psb_dpk_),target :: work(:) real(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_d_vect_type),intent(inout), optional :: u type(psb_d_vect_type),intent(inout), optional :: u
type(psb_d_vect_type) :: res type(psb_d_vect_type) :: res
type(psb_d_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables ! Local variables
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name character(len=20) :: name
name = 'inner_ml_aply' name = 'inner_ml_aply'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1167,57 +1169,9 @@ contains
goto 9999 goto 9999
case(mld_add_ml_) case(mld_add_ml_)
!
! Additive multilevel
!
if (level > 1) then call mld_d_inner_add(p, mlprec_wrk, level, trans, work)
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply')
goto 9999
end if
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -1230,164 +1184,108 @@ contains
select case(p%precv(level)%parms%smoother_pos) select case(p%precv(level)%parms%smoother_pos)
case(mld_post_smooth_) case(mld_post_smooth_)
p%precv(level)%parms%sweeps_pre = 0
call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
select case (trans_)
case('N')
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
! This is one step of post-smoothing case(mld_pre_smooth_)
p%precv(level)%parms%sweeps_post = 0
call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
if (level < nlev) then case(mld_twoside_smooth_)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
! case default
! Apply the prolongator info = psb_err_from_subroutine_ai_
! call psb_errpush(info,name,a_err='invalid smooth_pos',&
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/))
& dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999 goto 9999
end if
! end select
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
sweeps = p%precv(level)%parms%sweeps_post case(mld_mult_dev_ml_)
call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C') case(mld_vcycle_ml_, mld_wcycle_ml_)
! Post-smoothing transpose is pre-smoothing call mld_d_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
if (level > 1) then call mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u)
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then case default
call psb_errpush(psb_err_internal_error_,name,& info = psb_err_from_subroutine_ai_
& a_err='Error during restriction') call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999 goto 9999
end if
end select
call psb_erractionrestore(err_act)
return
end if 9999 call psb_error_handler(err_act)
return
! end subroutine inner_ml_aply
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
! recursive subroutine mld_d_inner_add(p, mlprec_wrk, level, trans, work)
! Compute the residual (at all levels but the coarsest one) use psb_base_mod
! use mld_prec_mod
if (level < nlev) then
call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) implicit none
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& done,mlprec_wrk(level)%vy2l,& integer(psb_ipk_), intent(in) :: level
& p%precv(level+1)%map,info,work=work) character, intent(in) :: trans
if (info /= psb_success_) then real(psb_dpk_),target :: work(:)
type(psb_d_vect_type) :: res
type(psb_d_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='wrong call level to inner_add')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if end if
case default if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end select end if
case(mld_pre_smooth_)
select case (trans_)
case('N')
! One step of pre-smoothing
if (level > 1) then if (level > 1) then
@ -1404,78 +1302,17 @@ contains
end if end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps sweeps = p%precv(level)%parms%sweeps
end if call p%precv(level)%sm%apply(done,&
call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,&
& p%precv(level)%base_desc, trans,& & sweeps,work,info)
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case('T','C')
! pre-smooth transpose is post-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during ADD smoother_apply')
goto 9999 goto 9999
end if end if
end if
if (level < nlev) then if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1488,7 +1325,7 @@ contains
! Apply the prolongator ! Apply the prolongator
! !
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& dzero,mlprec_wrk(level)%vy2l,& & done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1496,56 +1333,86 @@ contains
goto 9999 goto 9999
end if end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if end if
call psb_erractionrestore(err_act)
return
sweeps = p%precv(level)%parms%sweeps_pre 9999 call psb_error_handler(err_act)
call p%precv(level)%sm%apply(done,& return
& mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& end subroutine mld_d_inner_add
& sweeps,work,info)
if (info /= psb_success_) then recursive subroutine mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
call psb_errpush(psb_err_internal_error_,name,& use psb_base_mod
& a_err='Error during PRE smoother_apply') use mld_prec_mod
goto 9999
end if implicit none
else
sweeps = p%precv(level)%parms%sweeps !Input/Oputput variables
call p%precv(level)%sm%apply(done,& type(mld_dprec_type), intent(inout) :: p
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& sweeps,work,info) integer(psb_ipk_), intent(in) :: level
if (info /= psb_success_) then character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
type(psb_d_vect_type) :: res
type(psb_d_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_mult'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during smoother_apply') & a_err='wrong call level to inner_mult')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level
end if end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_twoside_smooth_) if ((level < nlev).or.(nlev == 1)) then
sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre
else
sweeps_post = p%precv(level-1)%parms%sweeps_post
sweeps_pre = p%precv(level-1)%parms%sweeps_pre
endif
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
if (level > 1) then if (level > 1) then
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vty,& if (pre) then
current => mlprec_wrk(level-1)%vty
else
current => mlprec_wrk(level-1)%vx2l
endif
call psb_map_X2Y(done,current,&
& dzero,mlprec_wrk(level)%vx2l,& & dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work) & p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -1553,13 +1420,14 @@ contains
end if end if
end if end if
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,& if (level < nlev) then
& p%precv(level)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (level < nlev) then
if (pre) then
if (trans == 'N') then if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
@ -1573,26 +1441,18 @@ contains
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info)
end if end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply') & a_err='Error during 2-PRE smoother_apply')
goto 9999 goto 9999
end if end if
endif
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual and call recursively
! and call recursively
! !
if(level < nlev) then if (pre) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,& & dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
@ -1605,6 +1465,7 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1617,10 +1478,16 @@ contains
! !
! Apply the prolongator ! Apply the prolongator
! !
if (pre) then
par = done
else
par = dzero
endif
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,& & par,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -1630,6 +1497,10 @@ contains
! !
! Compute the residual ! Compute the residual
! !
if (post) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& & done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
@ -1660,19 +1531,88 @@ contains
& a_err='Error during 2-POST smoother_apply') & a_err='Error during 2-POST smoother_apply')
goto 9999 goto 9999
end if end if
endif endif
case default else if (level == nlev) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid smooth_pos',& sweeps = p%precv(level)%parms%sweeps
& i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end if
end select call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_inner_mult
recursive subroutine mld_d_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
type(psb_d_vect_type),intent(inout), optional :: u
type(psb_d_vect_type) :: res
type(psb_d_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle !V/W cycle
if (level > 1) then if (level > 1) then
@ -1703,7 +1643,7 @@ contains
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(dzero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1766,7 +1706,7 @@ contains
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, u=mlprec_wrk(level+1)%vy2l)
endif endif
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1826,8 +1766,69 @@ contains
endif endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_d_inner_vw_cycle
recursive subroutine mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
type(psb_d_vect_type),intent(inout), optional :: u
type(psb_d_vect_type) :: res
type(psb_d_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
!K cycle !K cycle
@ -1844,7 +1845,7 @@ contains
& dzero,mlprec_wrk(level)%vy2l,& & dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(dzero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1984,23 +1985,13 @@ contains
endif endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine inner_ml_aply end subroutine mld_d_inner_k_cycle
recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)

645
mlprec/impl/mld_smlprec_aply.f90
Unescape Escape View File

@ -84,7 +84,7 @@
! The multilevel preconditioner data structure containing the ! The multilevel preconditioner data structure containing the
! local part of the preconditioner to be applied. ! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! p%precv(ilev)%prec - type(psb_dbaseprec_type) ! p%precv(ilev)%prec - type(psb_sbaseprec_type)
! The 'base preconditioner' for the current level ! The 'base preconditioner' for the current level
! p%precv(ilev)%ac - type(psb_sspmat_type) ! p%precv(ilev)%ac - type(psb_sspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
@ -333,7 +333,7 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_smlprec_aply' name='mld_smlprec_aply'
info = psb_success_ info = psb_success_
@ -1016,11 +1016,12 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act
character(len=20) :: name character(len=20) :: name
character :: trans_ character :: trans_
real(psb_spk_) :: par
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_smlprec_aply' name='mld_smlprec_aply'
info = psb_success_ info = psb_success_
@ -1036,7 +1037,6 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1065,12 +1065,12 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999 goto 9999
end if end if
end do end do
level = 1 level = 1
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call mlprec_wrk(level)%vy2l%zero() call mlprec_wrk(level)%vy2l%zero()
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1082,6 +1082,7 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
do level = 1, nlev do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info) call mlprec_wrk(level)%vtx%free(info)
@ -1118,29 +1119,30 @@ contains
! Arguments ! Arguments
integer(psb_ipk_) :: level integer(psb_ipk_) :: level
type(mld_sprec_type), target, intent(inout) :: p type(mld_sprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans character, intent(in) :: trans
real(psb_spk_),target :: work(:) real(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_s_vect_type),intent(inout), optional :: u type(psb_s_vect_type),intent(inout), optional :: u
type(psb_s_vect_type) :: res type(psb_s_vect_type) :: res
type(psb_s_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables ! Local variables
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name character(len=20) :: name
name = 'inner_ml_aply' name = 'inner_ml_aply'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1167,57 +1169,9 @@ contains
goto 9999 goto 9999
case(mld_add_ml_) case(mld_add_ml_)
!
! Additive multilevel
!
if (level > 1) then call mld_s_inner_add(p, mlprec_wrk, level, trans, work)
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,&
& szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply')
goto 9999
end if
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -1230,164 +1184,108 @@ contains
select case(p%precv(level)%parms%smoother_pos) select case(p%precv(level)%parms%smoother_pos)
case(mld_post_smooth_) case(mld_post_smooth_)
p%precv(level)%parms%sweeps_pre = 0
call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
select case (trans_)
case('N')
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,&
& szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
! This is one step of post-smoothing case(mld_pre_smooth_)
p%precv(level)%parms%sweeps_post = 0
call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
if (level < nlev) then case(mld_twoside_smooth_)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
! case default
! Apply the prolongator info = psb_err_from_subroutine_ai_
! call psb_errpush(info,name,a_err='invalid smooth_pos',&
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/))
& szero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999 goto 9999
end if
! end select
! Compute the residual
!
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
sweeps = p%precv(level)%parms%sweeps_post case(mld_mult_dev_ml_)
call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C') case(mld_vcycle_ml_, mld_wcycle_ml_)
! Post-smoothing transpose is pre-smoothing call mld_s_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
if (level > 1) then call mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u)
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,&
& szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then case default
call psb_errpush(psb_err_internal_error_,name,& info = psb_err_from_subroutine_ai_
& a_err='Error during restriction') call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999 goto 9999
end if
end select
call psb_erractionrestore(err_act)
return
end if 9999 call psb_error_handler(err_act)
return
! end subroutine inner_ml_aply
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
! recursive subroutine mld_s_inner_add(p, mlprec_wrk, level, trans, work)
! Compute the residual (at all levels but the coarsest one) use psb_base_mod
! use mld_prec_mod
if (level < nlev) then
call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) implicit none
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& sone,mlprec_wrk(level)%vy2l,& integer(psb_ipk_), intent(in) :: level
& p%precv(level+1)%map,info,work=work) character, intent(in) :: trans
if (info /= psb_success_) then real(psb_spk_),target :: work(:)
type(psb_s_vect_type) :: res
type(psb_s_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='wrong call level to inner_add')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if end if
case default if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end select end if
case(mld_pre_smooth_)
select case (trans_)
case('N')
! One step of pre-smoothing
if (level > 1) then if (level > 1) then
@ -1404,78 +1302,17 @@ contains
end if end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps sweeps = p%precv(level)%parms%sweeps
end if call p%precv(level)%sm%apply(sone,&
call p%precv(level)%sm%apply(sone,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,&
& p%precv(level)%base_desc, trans,& & sweeps,work,info)
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case('T','C')
! pre-smooth transpose is post-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,&
& szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during ADD smoother_apply')
goto 9999 goto 9999
end if end if
end if
if (level < nlev) then if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1488,7 +1325,7 @@ contains
! Apply the prolongator ! Apply the prolongator
! !
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& szero,mlprec_wrk(level)%vy2l,& & sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1496,56 +1333,86 @@ contains
goto 9999 goto 9999
end if end if
!
! Compute the residual
!
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if end if
call psb_erractionrestore(err_act)
return
sweeps = p%precv(level)%parms%sweeps_pre 9999 call psb_error_handler(err_act)
call p%precv(level)%sm%apply(sone,& return
& mlprec_wrk(level)%vx2l,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& end subroutine mld_s_inner_add
& sweeps,work,info)
if (info /= psb_success_) then recursive subroutine mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
call psb_errpush(psb_err_internal_error_,name,& use psb_base_mod
& a_err='Error during PRE smoother_apply') use mld_prec_mod
goto 9999
end if implicit none
else
sweeps = p%precv(level)%parms%sweeps !Input/Oputput variables
call p%precv(level)%sm%apply(sone,& type(mld_sprec_type), intent(inout) :: p
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& sweeps,work,info) integer(psb_ipk_), intent(in) :: level
if (info /= psb_success_) then character, intent(in) :: trans
real(psb_spk_),target :: work(:)
type(psb_s_vect_type) :: res
type(psb_s_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_mult'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during smoother_apply') & a_err='wrong call level to inner_mult')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level
end if end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_twoside_smooth_) if ((level < nlev).or.(nlev == 1)) then
sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre
else
sweeps_post = p%precv(level-1)%parms%sweeps_post
sweeps_pre = p%precv(level-1)%parms%sweeps_pre
endif
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
if (level > 1) then if (level > 1) then
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vty,& if (pre) then
current => mlprec_wrk(level-1)%vty
else
current => mlprec_wrk(level-1)%vx2l
endif
call psb_map_X2Y(sone,current,&
& szero,mlprec_wrk(level)%vx2l,& & szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work) & p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -1553,13 +1420,14 @@ contains
end if end if
end if end if
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vtx,& if (level < nlev) then
& p%precv(level)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (level < nlev) then
if (pre) then
if (trans == 'N') then if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
@ -1573,26 +1441,18 @@ contains
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info)
end if end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply') & a_err='Error during 2-PRE smoother_apply')
goto 9999 goto 9999
end if end if
endif
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual and call recursively
! and call recursively
! !
if(level < nlev) then if (pre) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,& & szero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
@ -1605,6 +1465,7 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1617,10 +1478,16 @@ contains
! !
! Apply the prolongator ! Apply the prolongator
! !
if (pre) then
par = sone
else
par = szero
endif
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,& & par,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -1630,6 +1497,10 @@ contains
! !
! Compute the residual ! Compute the residual
! !
if (post) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& & sone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
@ -1660,19 +1531,88 @@ contains
& a_err='Error during 2-POST smoother_apply') & a_err='Error during 2-POST smoother_apply')
goto 9999 goto 9999
end if end if
endif endif
case default else if (level == nlev) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid smooth_pos',& sweeps = p%precv(level)%parms%sweeps
& i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end if
end select call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_inner_mult
recursive subroutine mld_s_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
type(psb_s_vect_type),intent(inout), optional :: u
type(psb_s_vect_type) :: res
type(psb_s_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle !V/W cycle
if (level > 1) then if (level > 1) then
@ -1703,7 +1643,7 @@ contains
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(szero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1766,7 +1706,7 @@ contains
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, u=mlprec_wrk(level+1)%vy2l)
endif endif
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1826,8 +1766,69 @@ contains
endif endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_s_inner_vw_cycle
recursive subroutine mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
type(psb_s_vect_type),intent(inout), optional :: u
type(psb_s_vect_type) :: res
type(psb_s_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
!K cycle !K cycle
@ -1844,7 +1845,7 @@ contains
& szero,mlprec_wrk(level)%vy2l,& & szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(szero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1984,23 +1985,13 @@ contains
endif endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine inner_ml_aply end subroutine mld_s_inner_k_cycle
recursive subroutine mld_sinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) recursive subroutine mld_sinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)

647
mlprec/impl/mld_zmlprec_aply.f90
Unescape Escape View File

@ -39,7 +39,7 @@
! File: mld_zmlprec_aply.f90 ! File: mld_zmlprec_aply.f90
! !
! Subroutine: mld_zmlprec_aply ! Subroutine: mld_zmlprec_aply
! Version: complex ! Version: real
! !
! This routine computes ! This routine computes
! !
@ -84,7 +84,7 @@
! The multilevel preconditioner data structure containing the ! The multilevel preconditioner data structure containing the
! local part of the preconditioner to be applied. ! local part of the preconditioner to be applied.
! Note that nlev = size(p%precv) = number of levels. ! Note that nlev = size(p%precv) = number of levels.
! p%precv(ilev)%prec - type(psb_dbaseprec_type) ! p%precv(ilev)%prec - type(psb_zbaseprec_type)
! The 'base preconditioner' for the current level ! The 'base preconditioner' for the current level
! p%precv(ilev)%ac - type(psb_zspmat_type) ! p%precv(ilev)%ac - type(psb_zspmat_type)
! The local part of the matrix A(ilev). ! The local part of the matrix A(ilev).
@ -333,7 +333,7 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_zmlprec_aply' name='mld_zmlprec_aply'
info = psb_success_ info = psb_success_
@ -1016,11 +1016,12 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act
character(len=20) :: name character(len=20) :: name
character :: trans_ character :: trans_
complex(psb_dpk_) :: par
type mld_mlprec_wrk_type type mld_mlprec_wrk_type
complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_zmlprec_aply' name='mld_zmlprec_aply'
info = psb_success_ info = psb_success_
@ -1036,7 +1037,6 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv) & ' Entry ', size(p%precv)
trans_ = psb_toupper(trans) trans_ = psb_toupper(trans)
nlev = size(p%precv) nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info) allocate(mlprec_wrk(nlev),stat=info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1065,12 +1065,12 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999 goto 9999
end if end if
end do end do
level = 1 level = 1
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call mlprec_wrk(level)%vy2l%zero() call mlprec_wrk(level)%vy2l%zero()
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1082,6 +1082,7 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
do level = 1, nlev do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info) call mlprec_wrk(level)%vtx%free(info)
@ -1118,29 +1119,30 @@ contains
! Arguments ! Arguments
integer(psb_ipk_) :: level integer(psb_ipk_) :: level
type(mld_zprec_type), target, intent(inout) :: p type(mld_zprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans character, intent(in) :: trans
complex(psb_dpk_),target :: work(:) complex(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
type(psb_z_vect_type),intent(inout), optional :: u type(psb_z_vect_type),intent(inout), optional :: u
type(psb_z_vect_type) :: res type(psb_z_vect_type) :: res
type(psb_z_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables ! Local variables
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name character(len=20) :: name
name = 'inner_ml_aply' name = 'inner_ml_aply'
info = psb_success_ info = psb_success_
call psb_erractionsave(err_act) call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit() debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level() debug_level = psb_get_debug_level()
nlev = size(p%precv) nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1167,57 +1169,9 @@ contains
goto 9999 goto 9999
case(mld_add_ml_) case(mld_add_ml_)
!
! Additive multilevel
!
if (level > 1) then call mld_z_inner_add(p, mlprec_wrk, level, trans, work)
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,&
& zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply')
goto 9999
end if
if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case(mld_mult_ml_) case(mld_mult_ml_)
! !
@ -1230,164 +1184,108 @@ contains
select case(p%precv(level)%parms%smoother_pos) select case(p%precv(level)%parms%smoother_pos)
case(mld_post_smooth_) case(mld_post_smooth_)
p%precv(level)%parms%sweeps_pre = 0
call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
select case (trans_)
case('N')
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,&
& zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
! This is one step of post-smoothing case(mld_pre_smooth_)
p%precv(level)%parms%sweeps_post = 0
call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
if (level < nlev) then case(mld_twoside_smooth_)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
! case default
! Apply the prolongator info = psb_err_from_subroutine_ai_
! call psb_errpush(info,name,a_err='invalid smooth_pos',&
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/))
& zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999 goto 9999
end if
! end select
! Compute the residual
!
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
sweeps = p%precv(level)%parms%sweeps_post case(mld_mult_dev_ml_)
call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
else call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C') case(mld_vcycle_ml_, mld_wcycle_ml_)
! Post-smoothing transpose is pre-smoothing call mld_z_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
if (level > 1) then call mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u)
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,&
& zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then case default
call psb_errpush(psb_err_internal_error_,name,& info = psb_err_from_subroutine_ai_
& a_err='Error during restriction') call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999 goto 9999
end if
end select
call psb_erractionrestore(err_act)
return
end if 9999 call psb_error_handler(err_act)
return
! end subroutine inner_ml_aply
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_post
else
sweeps = p%precv(level)%parms%sweeps
end if
call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during POST smoother_apply')
goto 9999
end if
! recursive subroutine mld_z_inner_add(p, mlprec_wrk, level, trans, work)
! Compute the residual (at all levels but the coarsest one) use psb_base_mod
! use mld_prec_mod
if (level < nlev) then
call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) implicit none
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& zone,mlprec_wrk(level)%vy2l,& integer(psb_ipk_), intent(in) :: level
& p%precv(level+1)%map,info,work=work) character, intent(in) :: trans
if (info /= psb_success_) then complex(psb_dpk_),target :: work(:)
type(psb_z_vect_type) :: res
type(psb_z_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='wrong call level to inner_add')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if end if
case default if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_ info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans') call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end select end if
case(mld_pre_smooth_)
select case (trans_)
case('N')
! One step of pre-smoothing
if (level > 1) then if (level > 1) then
@ -1404,78 +1302,17 @@ contains
end if end if
!
! Apply the base preconditioner
!
if (level < nlev) then
sweeps = p%precv(level)%parms%sweeps_pre
else
sweeps = p%precv(level)%parms%sweeps sweeps = p%precv(level)%parms%sweeps
end if call p%precv(level)%sm%apply(zone,&
call p%precv(level)%sm%apply(zone,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,&
& p%precv(level)%base_desc, trans,& & sweeps,work,info)
& sweeps,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
!
if (level < nlev) then
call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
end if
case('T','C')
! pre-smooth transpose is post-smoothing
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,&
& zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during ADD smoother_apply')
goto 9999 goto 9999
end if end if
end if
if (level < nlev) then if (level < nlev) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1488,7 +1325,7 @@ contains
! Apply the prolongator ! Apply the prolongator
! !
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zzero,mlprec_wrk(level)%vy2l,& & zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1496,56 +1333,86 @@ contains
goto 9999 goto 9999
end if end if
!
! Compute the residual
!
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if end if
call psb_erractionrestore(err_act)
return
sweeps = p%precv(level)%parms%sweeps_pre 9999 call psb_error_handler(err_act)
call p%precv(level)%sm%apply(zone,& return
& mlprec_wrk(level)%vx2l,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& end subroutine mld_z_inner_add
& sweeps,work,info)
if (info /= psb_success_) then recursive subroutine mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
call psb_errpush(psb_err_internal_error_,name,& use psb_base_mod
& a_err='Error during PRE smoother_apply') use mld_prec_mod
goto 9999
end if implicit none
else
sweeps = p%precv(level)%parms%sweeps !Input/Oputput variables
call p%precv(level)%sm%apply(zone,& type(mld_zprec_type), intent(inout) :: p
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
& sweeps,work,info) integer(psb_ipk_), intent(in) :: level
if (info /= psb_success_) then character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
type(psb_z_vect_type) :: res
type(psb_z_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_mult'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during smoother_apply') & a_err='wrong call level to inner_mult')
goto 9999 goto 9999
end if end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_mult at level ',level
end if end if
case default
info = psb_err_internal_error_
call psb_errpush(info,name,a_err='invalid trans')
goto 9999
end select
case(mld_twoside_smooth_) if ((level < nlev).or.(nlev == 1)) then
sweeps_post = p%precv(level)%parms%sweeps_post
sweeps_pre = p%precv(level)%parms%sweeps_pre
else
sweeps_post = p%precv(level-1)%parms%sweeps_post
sweeps_pre = p%precv(level-1)%parms%sweeps_pre
endif
pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N'))
post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N'))
if (level > 1) then if (level > 1) then
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vty,& if (pre) then
current => mlprec_wrk(level-1)%vty
else
current => mlprec_wrk(level-1)%vx2l
endif
call psb_map_X2Y(zone,current,&
& zzero,mlprec_wrk(level)%vx2l,& & zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work) & p%precv(level)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -1553,13 +1420,14 @@ contains
end if end if
end if end if
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vtx,& if (level < nlev) then
& p%precv(level)%base_desc,info)
! !
! Apply the base preconditioner ! Apply the base preconditioner
! !
if (level < nlev) then
if (pre) then
if (trans == 'N') then if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
@ -1573,26 +1441,18 @@ contains
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info)
end if end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply') & a_err='Error during 2-PRE smoother_apply')
goto 9999 goto 9999
end if end if
endif
! !
! Compute the residual (at all levels but the coarsest one) ! Compute the residual and call recursively
! and call recursively
! !
if(level < nlev) then if (pre) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,& & zzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
@ -1605,6 +1465,7 @@ contains
& a_err='Error during residue') & a_err='Error during residue')
goto 9999 goto 9999
end if end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1617,10 +1478,16 @@ contains
! !
! Apply the prolongator ! Apply the prolongator
! !
if (pre) then
par = zone
else
par = zzero
endif
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,& & par,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -1630,6 +1497,10 @@ contains
! !
! Compute the residual ! Compute the residual
! !
if (post) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& & zone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans) & work=work,trans=trans)
@ -1660,19 +1531,88 @@ contains
& a_err='Error during 2-POST smoother_apply') & a_err='Error during 2-POST smoother_apply')
goto 9999 goto 9999
end if end if
endif endif
case default else if (level == nlev) then
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid smooth_pos',& sweeps = p%precv(level)%parms%sweeps
& i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999 goto 9999
end if
end select call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_inner_mult
recursive subroutine mld_z_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
type(psb_z_vect_type),intent(inout), optional :: u
type(psb_z_vect_type) :: res
type(psb_z_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle !V/W cycle
if (level > 1) then if (level > 1) then
@ -1703,7 +1643,7 @@ contains
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(zzero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1766,7 +1706,7 @@ contains
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, u=mlprec_wrk(level+1)%vy2l)
endif endif
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1826,8 +1766,69 @@ contains
endif endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_) call psb_erractionrestore(err_act)
return
9999 call psb_error_handler(err_act)
return
end subroutine mld_z_inner_vw_cycle
recursive subroutine mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
implicit none
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
type(psb_z_vect_type),intent(inout), optional :: u
type(psb_z_vect_type) :: res
type(psb_z_vect_type), pointer :: current
integer(psb_ipk_) :: sweeps_post, sweeps_pre
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
integer(psb_ipk_) :: debug_level, debug_unit
integer(psb_ipk_) :: nlev, ilev, sweeps
logical :: pre, post
character(len=20) :: name
name = 'inner_inner_add'
info = psb_success_
call psb_erractionsave(err_act)
debug_unit = psb_get_debug_unit()
debug_level = psb_get_debug_level()
nlev = size(p%precv)
if ((level < 1) .or. (level > nlev)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong call level to inner_add')
goto 9999
end if
ictxt = p%precv(level)%base_desc%get_context()
call psb_info(ictxt, me, np)
nc2l = p%precv(level)%base_desc%get_local_cols()
nr2l = p%precv(level)%base_desc%get_local_rows()
if(debug_level > 1) then
write(debug_unit,*) me,' inner_add at level ',level
end if
if ((level<1).or.(level>nlev)) then
info = psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='Invalid LEVEL>NLEV')
goto 9999
end if
!K cycle !K cycle
@ -1844,7 +1845,7 @@ contains
& zzero,mlprec_wrk(level)%vy2l,& & zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info) & p%precv(level)%base_desc,info)
else else
call mlprec_wrk(level)%vy2l%set(zzero) call mlprec_wrk(level)%vy2l%zero()
endif endif
res = mlprec_wrk(level)%vx2l res = mlprec_wrk(level)%vx2l
@ -1984,23 +1985,13 @@ contains
endif endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
& i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/))
goto 9999
end select
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
9999 call psb_error_handler(err_act) 9999 call psb_error_handler(err_act)
return return
end subroutine inner_ml_aply end subroutine mld_z_inner_k_cycle
recursive subroutine mld_zinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) recursive subroutine mld_zinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)

3
mlprec/mld_base_prec_type.F90
Unescape Escape View File

@ -227,6 +227,7 @@ module mld_base_prec_type
integer(psb_ipk_), parameter :: mld_kcycle_ml_ = 5 integer(psb_ipk_), parameter :: mld_kcycle_ml_ = 5
integer(psb_ipk_), parameter :: mld_kcyclesym_ml_ = 6 integer(psb_ipk_), parameter :: mld_kcyclesym_ml_ = 6
integer(psb_ipk_), parameter :: mld_new_ml_prec_ = 7 integer(psb_ipk_), parameter :: mld_new_ml_prec_ = 7
integer(psb_ipk_), parameter :: mld_mult_dev_ml_ = 7
integer(psb_ipk_), parameter :: mld_max_ml_type_ = 8 integer(psb_ipk_), parameter :: mld_max_ml_type_ = 8
! !
! Legal values for entry: mld_smoother_pos_ ! Legal values for entry: mld_smoother_pos_
@ -424,6 +425,8 @@ contains
val = mld_diag_scale_ val = mld_diag_scale_
case('ADD') case('ADD')
val = mld_add_ml_ val = mld_add_ml_
case('MULT_DEV')
val = mld_mult_dev_ml_
case('MULT') case('MULT')
val = mld_mult_ml_ val = mld_mult_ml_
case('VCYCLE') case('VCYCLE')

Write Preview
Loading…
Cancel
Save