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

Update mlprec_aply to use WRK.

Browse Source
This commit is contained in:
Salvatore Filippone
2017-12-06 16:11:53 +00:00
parent f1f3240f27
commit 5e174d062e
8 changed files with 608 additions and 768 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/impl/mld_cmlprec_aply.f90
+150 -191
View File
@@ -225,11 +225,8 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name
character :: trans_
complex(psb_spk_) :: beta_
type mld_mlprec_wrk_type
complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
logical :: do_alloc_wrk
type(mld_cmlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_cmlprec_aply'
info = psb_success_
@@ -245,34 +242,15 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
nlev = size(p%precv)
do_alloc_wrk = .not.allocated(p%wrk)
if (do_alloc_wrk) call p%allocate_wrk(info,vmold=x%v)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
end do
!
! At first iteration we must use the input BETA
!
@@ -280,31 +258,35 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,x,czero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
end if
do isweep = 1, p%outer_sweeps - 1
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = cone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,&
call psb_geaxpby(cone,x,czero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,y,&
& cone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
& cone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
end do
!
@@ -314,40 +296,24 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info)
call mlprec_wrk(level)%vty%free(info)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
nc2l = p%precv(level)%base_desc%get_local_cols()
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_spk_)')
goto 9999
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
goto 9999
end if
if (do_alloc_wrk) call p%free_wrk(info)
call psb_erractionrestore(err_act)
return
@@ -379,14 +345,13 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_cprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -419,7 +384,7 @@ contains
call psb_info(ictxt, me, np)
if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if
select case(p%precv(level)%parms%ml_cycle)
@@ -434,15 +399,15 @@ contains
case(mld_add_ml_)
call mld_c_inner_add(p, mlprec_wrk, level, trans, work)
call mld_c_inner_add(p, level, trans, work)
case(mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_c_inner_mult(p, level, trans, work)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
call mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work)
call mld_c_inner_k_cycle(p, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -464,7 +429,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_c_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_c_inner_add(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -473,7 +438,6 @@ contains
!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(:)
@@ -517,18 +481,18 @@ contains
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vy2l,czero,mlprec_wrk(level)%vtx,&
& p%wrk(level)%vy2l,czero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(cone,&
& mlprec_wrk(level)%vtx,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vtx,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
end do
@@ -536,7 +500,7 @@ contains
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -548,8 +512,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -557,7 +521,7 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -567,8 +531,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -587,7 +551,7 @@ contains
end subroutine mld_c_inner_add
recursive subroutine mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_c_inner_mult(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -596,7 +560,6 @@ contains
!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(:)
@@ -633,7 +596,6 @@ contains
sweeps_pre = p%precv(level)%parms%sweeps_pre
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 < nlev) then
!
@@ -645,13 +607,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -662,25 +624,24 @@ contains
goto 9999
end if
endif
!
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& 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 psb_map_X2Y(cone,mlprec_wrk(level)%vty,&
& czero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -689,8 +650,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -699,13 +660,13 @@ contains
end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -715,14 +676,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(cone,mlprec_wrk(level)%vty,&
& czero,mlprec_wrk(level+1)%vx2l,&
if (info == psb_success_) call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -730,10 +691,10 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -746,12 +707,12 @@ contains
if (post) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
& p%wrk(level)%vy2l,&
& cone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -765,13 +726,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -788,7 +749,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
@@ -808,7 +769,7 @@ contains
end subroutine mld_c_inner_mult
recursive subroutine mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
recursive subroutine mld_c_inner_k_cycle(p, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
@@ -816,7 +777,6 @@ contains
!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(:)
@@ -870,7 +830,7 @@ contains
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
@@ -879,13 +839,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,czero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -901,12 +861,12 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,cone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -915,8 +875,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level)%vty,&
& czero,mlprec_wrk(level + 1)%vx2l,&
call psb_map_X2Y(cone,p%wrk(level)%vty,&
& czero,p%wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@@ -929,16 +889,16 @@ contains
if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
call mld_cinneritkcycle(p, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
call mld_cinneritkcycle(p, level + 1, trans, work, 'GCR')
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_cycle')
goto 9999
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level + 1 ,p,trans,work,info)
endif
if (info /= psb_success_) then
@@ -950,8 +910,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(cone,p%wrk(level+1)%vy2l,&
& cone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -963,11 +923,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,&
& czero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_spmm(-cone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& cone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -980,13 +940,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vty,cone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,cone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1014,7 +974,7 @@ contains
end subroutine mld_c_inner_k_cycle
recursive subroutine mld_cinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
recursive subroutine mld_cinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_c_inner_mod, mld_protect_name => mld_cmlprec_aply
@@ -1024,7 +984,6 @@ contains
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
@@ -1044,34 +1003,34 @@ contains
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,czero,rhs,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,czero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,czero,w,&
call psb_geaxpby(cone,p%wrk(level)%vx2l,czero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
@@ -1085,12 +1044,12 @@ contains
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call p%wrk(level)%vy2l%zero()
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,mlprec_wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(cone,p%wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1128,9 +1087,9 @@ contains
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(cone,w,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(cone,mlprec_wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(cone,w,czero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(cone,p%wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
@@ -1165,7 +1124,7 @@ contains
call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info)
endif
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(cone,x,czero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
@@ -1217,10 +1176,10 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level
character(len=20) :: name
character :: trans_
type mld_mlprec_wrk_type
type mld_mlwrk_type
complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
end type mld_mlwrk_type
type(mld_mlwrk_type), allocatable, target :: mlwrk(:)
name='mld_cmlprec_aply'
info = psb_success_
@@ -1238,7 +1197,7 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
allocate(mlwrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
@@ -1246,13 +1205,13 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%x2l,&
call psb_geasb(mlwrk(level)%x2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%y2l,&
call psb_geasb(mlwrk(level)%y2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%tx,&
call psb_geasb(mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%ty,&
call psb_geasb(mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
@@ -1263,10 +1222,10 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
end if
end do
mlprec_wrk(level)%x2l(:) = x(:)
mlprec_wrk(level)%y2l(:) = czero
mlwrk(level)%x2l(:) = x(:)
mlwrk(level)%y2l(:) = czero
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,mlwrk,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1274,7 +1233,7 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,&
call psb_geaxpby(alpha,mlwrk(level)%y2l,beta,y,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1315,14 +1274,14 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlwrk,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_cprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
type(mld_mlwrk_type), intent(inout), target :: mlwrk(:)
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -1370,15 +1329,15 @@ contains
case(mld_add_ml_)
call mld_c_inner_add(p, mlprec_wrk, level, trans, work)
call mld_c_inner_add(p, mlwrk, level, trans, work)
case(mld_mult_ml_, mld_vcycle_ml_, mld_wcycle_ml_)
call mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_c_inner_mult(p, mlwrk, level, trans, work)
! !$ case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !$
! !$ call mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work)
! !$ call mld_c_inner_k_cycle(p, mlwrk, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -1397,7 +1356,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_c_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_c_inner_add(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1406,7 +1365,7 @@ contains
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
@@ -1450,7 +1409,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%x2l,czero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,czero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
@@ -1461,17 +1420,17 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level)%x2l,&
& czero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(cone,mlwrk(level)%x2l,&
& czero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
mlprec_wrk(level+1)%y2l(:) = czero
mlwrk(level+1)%y2l(:) = czero
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -1481,8 +1440,8 @@ contains
!
! Apply the prolongator and add correction.
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%y2l,&
& cone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(cone,mlwrk(level+1)%y2l,&
& cone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1501,7 +1460,7 @@ contains
end subroutine mld_c_inner_add
recursive subroutine mld_c_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_c_inner_mult(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1510,7 +1469,7 @@ contains
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
@@ -1567,13 +1526,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%x2l,czero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,czero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%x2l,czero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,czero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
end if
@@ -1589,20 +1548,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(cone,mlprec_wrk(level)%x2l,&
& czero,mlprec_wrk(level)%ty,&
call psb_geaxpby(cone,mlwrk(level)%x2l,&
& czero,mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%y2l,cone,mlprec_wrk(level)%ty,&
& mlwrk(level)%y2l,cone,mlwrk(level)%ty,&
& 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 psb_map_X2Y(cone,mlprec_wrk(level)%ty,&
& czero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(cone,mlwrk(level)%ty,&
& czero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1611,8 +1570,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(cone,mlprec_wrk(level)%x2l,&
& czero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(cone,mlwrk(level)%x2l,&
& czero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1621,14 +1580,14 @@ contains
end if
endif
! First guess is zero
mlprec_wrk(level+1)%y2l(:) = czero
mlwrk(level+1)%y2l(:) = czero
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
! On second call will use output y2l as initial guess
if (info == psb_success_) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info == psb_success_) call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
endif
if (info /= psb_success_) then
@@ -1641,8 +1600,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%y2l,&
& cone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(cone,mlwrk(level+1)%y2l,&
& cone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1654,11 +1613,11 @@ contains
! Compute the residual
!
if (post) then
call psb_geaxpby(cone,mlprec_wrk(level)%x2l,&
& czero,mlprec_wrk(level)%tx,&
call psb_geaxpby(cone,mlwrk(level)%x2l,&
& czero,mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,&
& cone,mlprec_wrk(level)%tx,p%precv(level)%base_desc,info,&
call psb_spmm(-cone,p%precv(level)%base_a,mlwrk(level)%y2l,&
& cone,mlwrk(level)%tx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1671,13 +1630,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%tx,cone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,cone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%tx,cone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,cone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1694,7 +1653,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%x2l,czero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,czero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
mlprec/impl/mld_dmlprec_aply.f90
+150 -191
View File
@@ -225,11 +225,8 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name
character :: trans_
real(psb_dpk_) :: beta_
type mld_mlprec_wrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
logical :: do_alloc_wrk
type(mld_dmlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_dmlprec_aply'
info = psb_success_
@@ -245,34 +242,15 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
nlev = size(p%precv)
do_alloc_wrk = .not.allocated(p%wrk)
if (do_alloc_wrk) call p%allocate_wrk(info,vmold=x%v)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
end do
!
! At first iteration we must use the input BETA
!
@@ -280,31 +258,35 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,x,dzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
end if
do isweep = 1, p%outer_sweeps - 1
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = done
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,&
call psb_geaxpby(done,x,dzero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,y,&
& done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
& done,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
end do
!
@@ -314,40 +296,24 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info)
call mlprec_wrk(level)%vty%free(info)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
nc2l = p%precv(level)%base_desc%get_local_cols()
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
goto 9999
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
goto 9999
end if
if (do_alloc_wrk) call p%free_wrk(info)
call psb_erractionrestore(err_act)
return
@@ -379,14 +345,13 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_dprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -419,7 +384,7 @@ contains
call psb_info(ictxt, me, np)
if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if
select case(p%precv(level)%parms%ml_cycle)
@@ -434,15 +399,15 @@ contains
case(mld_add_ml_)
call mld_d_inner_add(p, mlprec_wrk, level, trans, work)
call mld_d_inner_add(p, level, trans, work)
case(mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_d_inner_mult(p, level, trans, work)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
call mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work)
call mld_d_inner_k_cycle(p, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -464,7 +429,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_d_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_d_inner_add(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -473,7 +438,6 @@ contains
!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(:)
@@ -517,18 +481,18 @@ contains
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vy2l,dzero,mlprec_wrk(level)%vtx,&
& p%wrk(level)%vy2l,dzero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(done,&
& mlprec_wrk(level)%vtx,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vtx,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
end do
@@ -536,7 +500,7 @@ contains
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -548,8 +512,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -557,7 +521,7 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -567,8 +531,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -587,7 +551,7 @@ contains
end subroutine mld_d_inner_add
recursive subroutine mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_d_inner_mult(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -596,7 +560,6 @@ contains
!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(:)
@@ -633,7 +596,6 @@ contains
sweeps_pre = p%precv(level)%parms%sweeps_pre
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 < nlev) then
!
@@ -645,13 +607,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -662,25 +624,24 @@ contains
goto 9999
end if
endif
!
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& 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 psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -689,8 +650,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -699,13 +660,13 @@ contains
end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -715,14 +676,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level+1)%vx2l,&
if (info == psb_success_) call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -730,10 +691,10 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -746,12 +707,12 @@ contains
if (post) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
& p%wrk(level)%vy2l,&
& done,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -765,13 +726,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vty,done,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vty,done,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -788,7 +749,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
@@ -808,7 +769,7 @@ contains
end subroutine mld_d_inner_mult
recursive subroutine mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
recursive subroutine mld_d_inner_k_cycle(p, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
@@ -816,7 +777,6 @@ contains
!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(:)
@@ -870,7 +830,7 @@ contains
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
@@ -879,13 +839,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,dzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -901,12 +861,12 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,done,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -915,8 +875,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level + 1)%vx2l,&
call psb_map_X2Y(done,p%wrk(level)%vty,&
& dzero,p%wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@@ -929,16 +889,16 @@ contains
if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
call mld_dinneritkcycle(p, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
call mld_dinneritkcycle(p, level + 1, trans, work, 'GCR')
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_cycle')
goto 9999
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level + 1 ,p,trans,work,info)
endif
if (info /= psb_success_) then
@@ -950,8 +910,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(done,p%wrk(level+1)%vy2l,&
& done,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -963,11 +923,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(done,p%wrk(level)%vx2l,&
& dzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_spmm(-done,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& done,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -980,13 +940,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vty,done,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vty,done,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,done,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1014,7 +974,7 @@ contains
end subroutine mld_d_inner_k_cycle
recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
recursive subroutine mld_dinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_d_inner_mod, mld_protect_name => mld_dmlprec_aply
@@ -1024,7 +984,6 @@ contains
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
@@ -1044,34 +1003,34 @@ contains
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,rhs,&
call psb_geaxpby(done,p%wrk(level)%vx2l,dzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,w,&
call psb_geaxpby(done,p%wrk(level)%vx2l,dzero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
@@ -1085,12 +1044,12 @@ contains
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call p%wrk(level)%vy2l%zero()
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(done,p%wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1128,9 +1087,9 @@ contains
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(done,w,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(done,w,dzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(done,p%wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
@@ -1165,7 +1124,7 @@ contains
call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info)
endif
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(done,x,dzero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
@@ -1217,10 +1176,10 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level
character(len=20) :: name
character :: trans_
type mld_mlprec_wrk_type
type mld_mlwrk_type
real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
end type mld_mlwrk_type
type(mld_mlwrk_type), allocatable, target :: mlwrk(:)
name='mld_dmlprec_aply'
info = psb_success_
@@ -1238,7 +1197,7 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
allocate(mlwrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
@@ -1246,13 +1205,13 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%x2l,&
call psb_geasb(mlwrk(level)%x2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%y2l,&
call psb_geasb(mlwrk(level)%y2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%tx,&
call psb_geasb(mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%ty,&
call psb_geasb(mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
@@ -1263,10 +1222,10 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
end if
end do
mlprec_wrk(level)%x2l(:) = x(:)
mlprec_wrk(level)%y2l(:) = dzero
mlwrk(level)%x2l(:) = x(:)
mlwrk(level)%y2l(:) = dzero
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,mlwrk,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1274,7 +1233,7 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,&
call psb_geaxpby(alpha,mlwrk(level)%y2l,beta,y,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1315,14 +1274,14 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlwrk,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_dprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
type(mld_mlwrk_type), intent(inout), target :: mlwrk(:)
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -1370,15 +1329,15 @@ contains
case(mld_add_ml_)
call mld_d_inner_add(p, mlprec_wrk, level, trans, work)
call mld_d_inner_add(p, mlwrk, level, trans, work)
case(mld_mult_ml_, mld_vcycle_ml_, mld_wcycle_ml_)
call mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_d_inner_mult(p, mlwrk, level, trans, work)
! !$ case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !$
! !$ call mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work)
! !$ call mld_d_inner_k_cycle(p, mlwrk, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -1397,7 +1356,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_d_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_d_inner_add(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1406,7 +1365,7 @@ contains
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
@@ -1450,7 +1409,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,dzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
@@ -1461,17 +1420,17 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%x2l,&
& dzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(done,mlwrk(level)%x2l,&
& dzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
mlprec_wrk(level+1)%y2l(:) = dzero
mlwrk(level+1)%y2l(:) = dzero
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -1481,8 +1440,8 @@ contains
!
! Apply the prolongator and add correction.
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%y2l,&
& done,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(done,mlwrk(level+1)%y2l,&
& done,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1501,7 +1460,7 @@ contains
end subroutine mld_d_inner_add
recursive subroutine mld_d_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_d_inner_mult(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1510,7 +1469,7 @@ contains
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_dpk_),target :: work(:)
@@ -1567,13 +1526,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,dzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,dzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
end if
@@ -1589,20 +1548,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(done,mlprec_wrk(level)%x2l,&
& dzero,mlprec_wrk(level)%ty,&
call psb_geaxpby(done,mlwrk(level)%x2l,&
& dzero,mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%y2l,done,mlprec_wrk(level)%ty,&
& mlwrk(level)%y2l,done,mlwrk(level)%ty,&
& 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 psb_map_X2Y(done,mlprec_wrk(level)%ty,&
& dzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(done,mlwrk(level)%ty,&
& dzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1611,8 +1570,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(done,mlprec_wrk(level)%x2l,&
& dzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(done,mlwrk(level)%x2l,&
& dzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1621,14 +1580,14 @@ contains
end if
endif
! First guess is zero
mlprec_wrk(level+1)%y2l(:) = dzero
mlwrk(level+1)%y2l(:) = dzero
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
! On second call will use output y2l as initial guess
if (info == psb_success_) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info == psb_success_) call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
endif
if (info /= psb_success_) then
@@ -1641,8 +1600,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%y2l,&
& done,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(done,mlwrk(level+1)%y2l,&
& done,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1654,11 +1613,11 @@ contains
! Compute the residual
!
if (post) then
call psb_geaxpby(done,mlprec_wrk(level)%x2l,&
& dzero,mlprec_wrk(level)%tx,&
call psb_geaxpby(done,mlwrk(level)%x2l,&
& dzero,mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%y2l,&
& done,mlprec_wrk(level)%tx,p%precv(level)%base_desc,info,&
call psb_spmm(-done,p%precv(level)%base_a,mlwrk(level)%y2l,&
& done,mlwrk(level)%tx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1671,13 +1630,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,done,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,done,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1694,7 +1653,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,dzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
mlprec/impl/mld_smlprec_aply.f90
+150 -191
View File
@@ -225,11 +225,8 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name
character :: trans_
real(psb_spk_) :: beta_
type mld_mlprec_wrk_type
real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
logical :: do_alloc_wrk
type(mld_smlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_smlprec_aply'
info = psb_success_
@@ -245,34 +242,15 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
nlev = size(p%precv)
do_alloc_wrk = .not.allocated(p%wrk)
if (do_alloc_wrk) call p%allocate_wrk(info,vmold=x%v)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
end do
!
! At first iteration we must use the input BETA
!
@@ -280,31 +258,35 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,x,szero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
end if
do isweep = 1, p%outer_sweeps - 1
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = sone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,&
call psb_geaxpby(sone,x,szero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,y,&
& sone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
& sone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
end do
!
@@ -314,40 +296,24 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info)
call mlprec_wrk(level)%vty%free(info)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
nc2l = p%precv(level)%base_desc%get_local_cols()
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_spk_)')
goto 9999
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
goto 9999
end if
if (do_alloc_wrk) call p%free_wrk(info)
call psb_erractionrestore(err_act)
return
@@ -379,14 +345,13 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_sprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -419,7 +384,7 @@ contains
call psb_info(ictxt, me, np)
if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if
select case(p%precv(level)%parms%ml_cycle)
@@ -434,15 +399,15 @@ contains
case(mld_add_ml_)
call mld_s_inner_add(p, mlprec_wrk, level, trans, work)
call mld_s_inner_add(p, level, trans, work)
case(mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_s_inner_mult(p, level, trans, work)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
call mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work)
call mld_s_inner_k_cycle(p, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -464,7 +429,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_s_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_s_inner_add(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -473,7 +438,6 @@ contains
!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(:)
@@ -517,18 +481,18 @@ contains
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vy2l,szero,mlprec_wrk(level)%vtx,&
& p%wrk(level)%vy2l,szero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(sone,&
& mlprec_wrk(level)%vtx,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vtx,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
end do
@@ -536,7 +500,7 @@ contains
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -548,8 +512,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -557,7 +521,7 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -567,8 +531,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -587,7 +551,7 @@ contains
end subroutine mld_s_inner_add
recursive subroutine mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_s_inner_mult(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -596,7 +560,6 @@ contains
!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(:)
@@ -633,7 +596,6 @@ contains
sweeps_pre = p%precv(level)%parms%sweeps_pre
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 < nlev) then
!
@@ -645,13 +607,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -662,25 +624,24 @@ contains
goto 9999
end if
endif
!
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& 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 psb_map_X2Y(sone,mlprec_wrk(level)%vty,&
& szero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -689,8 +650,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -699,13 +660,13 @@ contains
end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -715,14 +676,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(sone,mlprec_wrk(level)%vty,&
& szero,mlprec_wrk(level+1)%vx2l,&
if (info == psb_success_) call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -730,10 +691,10 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -746,12 +707,12 @@ contains
if (post) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
& p%wrk(level)%vy2l,&
& sone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -765,13 +726,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -788,7 +749,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
@@ -808,7 +769,7 @@ contains
end subroutine mld_s_inner_mult
recursive subroutine mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
recursive subroutine mld_s_inner_k_cycle(p, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
@@ -816,7 +777,6 @@ contains
!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(:)
@@ -870,7 +830,7 @@ contains
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
@@ -879,13 +839,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,szero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -901,12 +861,12 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,sone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -915,8 +875,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level)%vty,&
& szero,mlprec_wrk(level + 1)%vx2l,&
call psb_map_X2Y(sone,p%wrk(level)%vty,&
& szero,p%wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@@ -929,16 +889,16 @@ contains
if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
call mld_sinneritkcycle(p, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
call mld_sinneritkcycle(p, level + 1, trans, work, 'GCR')
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_cycle')
goto 9999
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level + 1 ,p,trans,work,info)
endif
if (info /= psb_success_) then
@@ -950,8 +910,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(sone,p%wrk(level+1)%vy2l,&
& sone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -963,11 +923,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,&
& szero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_spmm(-sone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& sone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -980,13 +940,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vty,sone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,sone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1014,7 +974,7 @@ contains
end subroutine mld_s_inner_k_cycle
recursive subroutine mld_sinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
recursive subroutine mld_sinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_s_inner_mod, mld_protect_name => mld_smlprec_aply
@@ -1024,7 +984,6 @@ contains
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
@@ -1044,34 +1003,34 @@ contains
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,szero,rhs,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,szero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,szero,w,&
call psb_geaxpby(sone,p%wrk(level)%vx2l,szero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
@@ -1085,12 +1044,12 @@ contains
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call p%wrk(level)%vy2l%zero()
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,mlprec_wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(sone,p%wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1128,9 +1087,9 @@ contains
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(sone,w,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(sone,mlprec_wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(sone,w,szero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(sone,p%wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
@@ -1165,7 +1124,7 @@ contains
call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info)
endif
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(sone,x,szero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
@@ -1217,10 +1176,10 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level
character(len=20) :: name
character :: trans_
type mld_mlprec_wrk_type
type mld_mlwrk_type
real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
end type mld_mlwrk_type
type(mld_mlwrk_type), allocatable, target :: mlwrk(:)
name='mld_smlprec_aply'
info = psb_success_
@@ -1238,7 +1197,7 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
allocate(mlwrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
@@ -1246,13 +1205,13 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%x2l,&
call psb_geasb(mlwrk(level)%x2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%y2l,&
call psb_geasb(mlwrk(level)%y2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%tx,&
call psb_geasb(mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%ty,&
call psb_geasb(mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
@@ -1263,10 +1222,10 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
end if
end do
mlprec_wrk(level)%x2l(:) = x(:)
mlprec_wrk(level)%y2l(:) = szero
mlwrk(level)%x2l(:) = x(:)
mlwrk(level)%y2l(:) = szero
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,mlwrk,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1274,7 +1233,7 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,&
call psb_geaxpby(alpha,mlwrk(level)%y2l,beta,y,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1315,14 +1274,14 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlwrk,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_sprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
type(mld_mlwrk_type), intent(inout), target :: mlwrk(:)
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -1370,15 +1329,15 @@ contains
case(mld_add_ml_)
call mld_s_inner_add(p, mlprec_wrk, level, trans, work)
call mld_s_inner_add(p, mlwrk, level, trans, work)
case(mld_mult_ml_, mld_vcycle_ml_, mld_wcycle_ml_)
call mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_s_inner_mult(p, mlwrk, level, trans, work)
! !$ case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !$
! !$ call mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work)
! !$ call mld_s_inner_k_cycle(p, mlwrk, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -1397,7 +1356,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_s_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_s_inner_add(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1406,7 +1365,7 @@ contains
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
@@ -1450,7 +1409,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%x2l,szero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,szero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
@@ -1461,17 +1420,17 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level)%x2l,&
& szero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(sone,mlwrk(level)%x2l,&
& szero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
mlprec_wrk(level+1)%y2l(:) = szero
mlwrk(level+1)%y2l(:) = szero
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -1481,8 +1440,8 @@ contains
!
! Apply the prolongator and add correction.
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%y2l,&
& sone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(sone,mlwrk(level+1)%y2l,&
& sone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1501,7 +1460,7 @@ contains
end subroutine mld_s_inner_add
recursive subroutine mld_s_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_s_inner_mult(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1510,7 +1469,7 @@ contains
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
@@ -1567,13 +1526,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%x2l,szero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,szero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%x2l,szero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,szero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
end if
@@ -1589,20 +1548,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(sone,mlprec_wrk(level)%x2l,&
& szero,mlprec_wrk(level)%ty,&
call psb_geaxpby(sone,mlwrk(level)%x2l,&
& szero,mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%y2l,sone,mlprec_wrk(level)%ty,&
& mlwrk(level)%y2l,sone,mlwrk(level)%ty,&
& 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 psb_map_X2Y(sone,mlprec_wrk(level)%ty,&
& szero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(sone,mlwrk(level)%ty,&
& szero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1611,8 +1570,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(sone,mlprec_wrk(level)%x2l,&
& szero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(sone,mlwrk(level)%x2l,&
& szero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1621,14 +1580,14 @@ contains
end if
endif
! First guess is zero
mlprec_wrk(level+1)%y2l(:) = szero
mlwrk(level+1)%y2l(:) = szero
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
! On second call will use output y2l as initial guess
if (info == psb_success_) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info == psb_success_) call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
endif
if (info /= psb_success_) then
@@ -1641,8 +1600,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%y2l,&
& sone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(sone,mlwrk(level+1)%y2l,&
& sone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1654,11 +1613,11 @@ contains
! Compute the residual
!
if (post) then
call psb_geaxpby(sone,mlprec_wrk(level)%x2l,&
& szero,mlprec_wrk(level)%tx,&
call psb_geaxpby(sone,mlwrk(level)%x2l,&
& szero,mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,&
& sone,mlprec_wrk(level)%tx,p%precv(level)%base_desc,info,&
call psb_spmm(-sone,p%precv(level)%base_a,mlwrk(level)%y2l,&
& sone,mlwrk(level)%tx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1671,13 +1630,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%tx,sone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,sone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%tx,sone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,sone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1694,7 +1653,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%x2l,szero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,szero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
mlprec/impl/mld_zmlprec_aply.f90
+150 -191
View File
@@ -225,11 +225,8 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
character(len=20) :: name
character :: trans_
complex(psb_dpk_) :: beta_
type mld_mlprec_wrk_type
complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
logical :: do_alloc_wrk
type(mld_zmlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
name='mld_zmlprec_aply'
info = psb_success_
@@ -245,34 +242,15 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
& ' Entry ', size(p%precv)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
nlev = size(p%precv)
do_alloc_wrk = .not.allocated(p%wrk)
if (do_alloc_wrk) call p%allocate_wrk(info,vmold=x%v)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
end if
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%vx2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
call psb_geasb(mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=x%v)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
end do
!
! At first iteration we must use the input BETA
!
@@ -280,31 +258,35 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='geaxbpy')
goto 9999
end if
do isweep = 1, p%outer_sweeps - 1
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
! all iterations after the first must use BETA = 1
beta_ = zone
!
! Next iteration should use the current residual to compute a correction
!
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,&
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vx2l,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,y,&
& zone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
& zone,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
end do
!
@@ -314,40 +296,24 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
!
! With the current implementation, y2l is zeroed internally at first smoother.
! call mlprec_wrk(level)%vy2l%zero()
! call p%wrk(level)%vy2l%zero()
!
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Inner prec aply')
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,&
call psb_geaxpby(alpha,p%wrk(level)%vy2l,beta_,y,&
& p%precv(level)%base_desc,info)
do level = 1, nlev
call mlprec_wrk(level)%vx2l%free(info)
call mlprec_wrk(level)%vy2l%free(info)
call mlprec_wrk(level)%vtx%free(info)
call mlprec_wrk(level)%vty%free(info)
if (psb_errstatus_fatal()) then
info=psb_err_alloc_request_
nc2l = p%precv(level)%base_desc%get_local_cols()
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='complex(psb_dpk_)')
goto 9999
end if
end do
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error final update')
goto 9999
end if
if (do_alloc_wrk) call p%free_wrk(info)
call psb_erractionrestore(err_act)
return
@@ -379,14 +345,13 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_zprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -419,7 +384,7 @@ contains
call psb_info(ictxt, me, np)
if(debug_level > 1) then
write(debug_unit,*) me,' Start inner_ml_aply at level ',level
write(debug_unit,*) me,' Start inner_ml_aply at level ',level, info
end if
select case(p%precv(level)%parms%ml_cycle)
@@ -434,15 +399,15 @@ contains
case(mld_add_ml_)
call mld_z_inner_add(p, mlprec_wrk, level, trans, work)
call mld_z_inner_add(p, level, trans, work)
case(mld_mult_ml_,mld_vcycle_ml_, mld_wcycle_ml_)
call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_z_inner_mult(p, level, trans, work)
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
call mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work)
call mld_z_inner_k_cycle(p, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -464,7 +429,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_z_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_z_inner_add(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -473,7 +438,6 @@ contains
!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(:)
@@ -517,18 +481,18 @@ contains
if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
sweeps = max(p%precv(level)%parms%sweeps_pre,p%precv(level)%parms%sweeps_post)
do k=1, sweeps
call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vy2l,zzero,mlprec_wrk(level)%vtx,&
& p%wrk(level)%vy2l,zzero,p%wrk(level)%vtx,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
call p%precv(level)%sm2a%apply(zone,&
& mlprec_wrk(level)%vtx,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vtx,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& ione,work,info,init='Z')
end do
@@ -536,7 +500,7 @@ contains
else
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -548,8 +512,8 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -557,7 +521,7 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -567,8 +531,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -587,7 +551,7 @@ contains
end subroutine mld_z_inner_add
recursive subroutine mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_z_inner_mult(p, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -596,7 +560,6 @@ contains
!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(:)
@@ -633,7 +596,6 @@ contains
sweeps_pre = p%precv(level)%parms%sweeps_pre
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 < nlev) then
!
@@ -645,13 +607,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -662,25 +624,24 @@ contains
goto 9999
end if
endif
!
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& 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 psb_map_X2Y(zone,mlprec_wrk(level)%vty,&
& zzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -689,8 +650,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level+1)%vx2l,&
call psb_map_X2Y(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -699,13 +660,13 @@ contains
end if
endif
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -715,14 +676,14 @@ contains
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(zone,mlprec_wrk(level)%vty,&
& zzero,mlprec_wrk(level+1)%vx2l,&
if (info == psb_success_) call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -730,10 +691,10 @@ contains
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,trans,work,info)
if (info == psb_success_) call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
if (info == psb_success_) call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -746,12 +707,12 @@ contains
if (post) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
& p%wrk(level)%vy2l,&
& zone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -765,13 +726,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -788,7 +749,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
@@ -808,7 +769,7 @@ contains
end subroutine mld_z_inner_mult
recursive subroutine mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work,u)
recursive subroutine mld_z_inner_k_cycle(p, level, trans, work,u)
use psb_base_mod
use mld_prec_mod
@@ -816,7 +777,6 @@ contains
!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(:)
@@ -870,7 +830,7 @@ contains
!
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
@@ -879,13 +839,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vx2l,zzero,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -901,12 +861,12 @@ contains
! Compute the residual and call recursively
!
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,&
& p%wrk(level)%vy2l,zone,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -915,8 +875,8 @@ contains
end if
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level)%vty,&
& zzero,mlprec_wrk(level + 1)%vx2l,&
call psb_map_X2Y(zone,p%wrk(level)%vty,&
& zzero,p%wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
@@ -929,16 +889,16 @@ contains
if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_cycle == mld_kcyclesym_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
call mld_zinneritkcycle(p, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_cycle == mld_kcycle_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
call mld_zinneritkcycle(p, level + 1, trans, work, 'GCR')
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_cycle')
goto 9999
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level + 1 ,p,trans,work,info)
endif
if (info /= psb_success_) then
@@ -950,8 +910,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
call psb_map_Y2X(zone,p%wrk(level+1)%vy2l,&
& zone,p%wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
@@ -963,11 +923,11 @@ contains
!
! Compute the residual
!
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,&
& zzero,p%wrk(level)%vty,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vty,p%precv(level)%base_desc,info,&
call psb_spmm(-zone,p%precv(level)%base_a,p%wrk(level)%vy2l,&
& zone,p%wrk(level)%vty,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -980,13 +940,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vty,zone,mlprec_wrk(level)%vy2l,&
& p%wrk(level)%vty,zone,p%wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1014,7 +974,7 @@ contains
end subroutine mld_z_inner_k_cycle
recursive subroutine mld_zinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
recursive subroutine mld_zinneritkcycle(p, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_z_inner_mod, mld_protect_name => mld_zmlprec_aply
@@ -1024,7 +984,6 @@ contains
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
@@ -1044,34 +1003,34 @@ contains
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
& scratch=.true.,mold=p%wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
& scratch=.true.,mold=p%wrk(level)%vy2l%v)
call x%zero()
! rhs=vx2l and w=rhs
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,zzero,rhs,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,zzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,zzero,w,&
call psb_geaxpby(zone,p%wrk(level)%vx2l,zzero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
@@ -1085,12 +1044,12 @@ contains
delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call p%wrk(level)%vy2l%zero()
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,mlprec_wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(zone,p%wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1128,9 +1087,9 @@ contains
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(zone,w,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(zone,mlprec_wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_geaxpby(zone,w,zzero,p%wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,trans,work,info)
call psb_geaxpby(zone,p%wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
@@ -1165,7 +1124,7 @@ contains
call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info)
endif
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_geaxpby(zone,x,zzero,p%wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
@@ -1217,10 +1176,10 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level
character(len=20) :: name
character :: trans_
type mld_mlprec_wrk_type
type mld_mlwrk_type
complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:)
end type mld_mlprec_wrk_type
type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:)
end type mld_mlwrk_type
type(mld_mlwrk_type), allocatable, target :: mlwrk(:)
name='mld_zmlprec_aply'
info = psb_success_
@@ -1238,7 +1197,7 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
trans_ = psb_toupper(trans)
nlev = size(p%precv)
allocate(mlprec_wrk(nlev),stat=info)
allocate(mlwrk(nlev),stat=info)
if (info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate')
goto 9999
@@ -1246,13 +1205,13 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
level = 1
do level = 1, nlev
call psb_geasb(mlprec_wrk(level)%x2l,&
call psb_geasb(mlwrk(level)%x2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%y2l,&
call psb_geasb(mlwrk(level)%y2l,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%tx,&
call psb_geasb(mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_geasb(mlprec_wrk(level)%ty,&
call psb_geasb(mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
@@ -1263,10 +1222,10 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
end if
end do
mlprec_wrk(level)%x2l(:) = x(:)
mlprec_wrk(level)%y2l(:) = zzero
mlwrk(level)%x2l(:) = x(:)
mlwrk(level)%y2l(:) = zzero
call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info)
call inner_ml_aply(level,p,mlwrk,trans_,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1274,7 +1233,7 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info)
goto 9999
end if
call psb_geaxpby(alpha,mlprec_wrk(level)%y2l,beta,y,&
call psb_geaxpby(alpha,mlwrk(level)%y2l,beta,y,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
@@ -1315,14 +1274,14 @@ contains
! between level and level+1 are stored at level+1.
!
!
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlwrk,trans,work,info)
implicit none
! Arguments
integer(psb_ipk_) :: level
type(mld_zprec_type), target, intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:)
type(mld_mlwrk_type), intent(inout), target :: mlwrk(:)
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
@@ -1370,15 +1329,15 @@ contains
case(mld_add_ml_)
call mld_z_inner_add(p, mlprec_wrk, level, trans, work)
call mld_z_inner_add(p, mlwrk, level, trans, work)
case(mld_mult_ml_, mld_vcycle_ml_, mld_wcycle_ml_)
call mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
call mld_z_inner_mult(p, mlwrk, level, trans, work)
! !$ case(mld_kcycle_ml_, mld_kcyclesym_ml_)
! !$
! !$ call mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work)
! !$ call mld_z_inner_k_cycle(p, mlwrk, level, trans, work)
case default
info = psb_err_from_subroutine_ai_
@@ -1397,7 +1356,7 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_z_inner_add(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_z_inner_add(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1406,7 +1365,7 @@ contains
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
@@ -1450,7 +1409,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,zzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
if (info /= psb_success_) then
@@ -1461,17 +1420,17 @@ contains
if (level < nlev) then
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level)%x2l,&
& zzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(zone,mlwrk(level)%x2l,&
& zzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
mlprec_wrk(level+1)%y2l(:) = zzero
mlwrk(level+1)%y2l(:) = zzero
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
@@ -1481,8 +1440,8 @@ contains
!
! Apply the prolongator and add correction.
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,&
& zone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(zone,mlwrk(level+1)%y2l,&
& zone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1501,7 +1460,7 @@ contains
end subroutine mld_z_inner_add
recursive subroutine mld_z_inner_mult(p, mlprec_wrk, level, trans, work)
recursive subroutine mld_z_inner_mult(p, mlwrk, level, trans, work)
use psb_base_mod
use mld_prec_mod
@@ -1510,7 +1469,7 @@ contains
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:)
type(mld_mlwrk_type), target, intent(inout) :: mlwrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
@@ -1567,13 +1526,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,zzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,zzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y')
end if
@@ -1589,20 +1548,20 @@ contains
! Compute the residual and call recursively
!
if (pre) then
call psb_geaxpby(zone,mlprec_wrk(level)%x2l,&
& zzero,mlprec_wrk(level)%ty,&
call psb_geaxpby(zone,mlwrk(level)%x2l,&
& zzero,mlwrk(level)%ty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%y2l,zone,mlprec_wrk(level)%ty,&
& mlwrk(level)%y2l,zone,mlwrk(level)%ty,&
& 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 psb_map_X2Y(zone,mlprec_wrk(level)%ty,&
& zzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(zone,mlwrk(level)%ty,&
& zzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1611,8 +1570,8 @@ contains
end if
else
! Shortcut: just transfer x2l.
call psb_map_X2Y(zone,mlprec_wrk(level)%x2l,&
& zzero,mlprec_wrk(level+1)%x2l,&
call psb_map_X2Y(zone,mlwrk(level)%x2l,&
& zzero,mlwrk(level+1)%x2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1621,14 +1580,14 @@ contains
end if
endif
! First guess is zero
mlprec_wrk(level+1)%y2l(:) = zzero
mlwrk(level+1)%y2l(:) = zzero
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
if (p%precv(level)%parms%ml_cycle == mld_wcycle_ml_) then
! On second call will use output y2l as initial guess
if (info == psb_success_) call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info == psb_success_) call inner_ml_aply(level+1,p,mlwrk,trans,work,info)
endif
if (info /= psb_success_) then
@@ -1641,8 +1600,8 @@ contains
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%y2l,&
& zone,mlprec_wrk(level)%y2l,&
call psb_map_Y2X(zone,mlwrk(level+1)%y2l,&
& zone,mlwrk(level)%y2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1654,11 +1613,11 @@ contains
! Compute the residual
!
if (post) then
call psb_geaxpby(zone,mlprec_wrk(level)%x2l,&
& zzero,mlprec_wrk(level)%tx,&
call psb_geaxpby(zone,mlwrk(level)%x2l,&
& zzero,mlwrk(level)%tx,&
& p%precv(level)%base_desc,info)
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%y2l,&
& zone,mlprec_wrk(level)%tx,p%precv(level)%base_desc,info,&
call psb_spmm(-zone,p%precv(level)%base_a,mlwrk(level)%y2l,&
& zone,mlwrk(level)%tx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
@@ -1671,13 +1630,13 @@ contains
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%tx,zone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,zone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%tx,zone,mlprec_wrk(level)%y2l,&
& mlwrk(level)%tx,zone,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Z')
end if
@@ -1694,7 +1653,7 @@ contains
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%x2l,zzero,mlprec_wrk(level)%y2l,&
& mlwrk(level)%x2l,zzero,mlwrk(level)%y2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
mlprec/mld_c_inner_mod.f90
+2 -1
View File
@@ -48,7 +48,8 @@ module mld_c_inner_mod
use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_i_base_vect_type, &
& psb_spk_, psb_c_base_sparse_mat, psb_c_base_vect_type, psb_ipk_, &
& psb_c_vect_type
use mld_c_prec_type, only : mld_cprec_type, mld_sml_parms, mld_c_onelev_type
use mld_c_prec_type, only : mld_cprec_type, mld_sml_parms, &
& mld_c_onelev_type, mld_cmlprec_wrk_type
interface mld_mlprec_bld
subroutine mld_cmlprec_bld(a,desc_a,prec,info, amold, vmold,imold)
mlprec/mld_d_inner_mod.f90
+2 -1
View File
@@ -48,7 +48,8 @@ module mld_d_inner_mod
use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_i_base_vect_type, &
& psb_dpk_, psb_d_base_sparse_mat, psb_d_base_vect_type, psb_ipk_, &
& psb_d_vect_type
use mld_d_prec_type, only : mld_dprec_type, mld_dml_parms, mld_d_onelev_type
use mld_d_prec_type, only : mld_dprec_type, mld_dml_parms, &
& mld_d_onelev_type, mld_dmlprec_wrk_type
interface mld_mlprec_bld
subroutine mld_dmlprec_bld(a,desc_a,prec,info, amold, vmold,imold)
mlprec/mld_s_inner_mod.f90
+2 -1
View File
@@ -48,7 +48,8 @@ module mld_s_inner_mod
use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_i_base_vect_type, &
& psb_spk_, psb_s_base_sparse_mat, psb_s_base_vect_type, psb_ipk_, &
& psb_s_vect_type
use mld_s_prec_type, only : mld_sprec_type, mld_sml_parms, mld_s_onelev_type
use mld_s_prec_type, only : mld_sprec_type, mld_sml_parms, &
& mld_s_onelev_type, mld_smlprec_wrk_type
interface mld_mlprec_bld
subroutine mld_smlprec_bld(a,desc_a,prec,info, amold, vmold,imold)
mlprec/mld_z_inner_mod.f90
+2 -1
View File
@@ -48,7 +48,8 @@ module mld_z_inner_mod
use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_i_base_vect_type, &
& psb_dpk_, psb_z_base_sparse_mat, psb_z_base_vect_type, psb_ipk_, &
& psb_z_vect_type
use mld_z_prec_type, only : mld_zprec_type, mld_dml_parms, mld_z_onelev_type
use mld_z_prec_type, only : mld_zprec_type, mld_dml_parms, &
& mld_z_onelev_type, mld_zmlprec_wrk_type
interface mld_mlprec_bld
subroutine mld_zmlprec_bld(a,desc_a,prec,info, amold, vmold,imold)