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

Fixed allocte_wrk & free_wrk for WV allocation.

Browse Source 
Modified interface of smoothers to use WV.
Initial tests.
Added WV to calls to MAP_X2Y & MAP_Y2X.
stopcriterion
Salvatore Filippone 7 years ago
parent 6f9a3c10d2
commit 823db4f943
40 changed files with 1212 additions and 1028 deletions
Show Stats Download Patch File Download Diff File

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

@ -255,7 +255,8 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(:))
! !
! At first iteration we must use the input BETA ! At first iteration we must use the input BETA
! !
@ -482,7 +483,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(cone,& call psb_geaxpby(cone,&
@ -494,12 +496,12 @@ contains
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& vy2l,czero,vty,& & vy2l,czero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(cone,& call p%precv(level)%sm2a%apply(cone,&
& vty,czero,vy2l,& & vty,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
@ -507,7 +509,7 @@ contains
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -519,7 +521,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(cone,vx2l,& call psb_map_X2Y(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,& & czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -538,7 +541,8 @@ contains
! !
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,& & cone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -605,7 +609,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level < nlev) then if (level < nlev) then
! !
! Apply the first smoother ! Apply the first smoother
@ -618,13 +623,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -651,7 +656,8 @@ contains
end if end if
call psb_map_X2Y(cone,vty,& call psb_map_X2Y(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,& & czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -661,7 +667,8 @@ contains
! Shortcut: just transfer x2l. ! Shortcut: just transfer x2l.
call psb_map_X2Y(cone,vx2l,& call psb_map_X2Y(cone,vx2l,&
& czero,p%precv(level+1)%wrk%vx2l,& & czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -676,7 +683,8 @@ contains
! !
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,& & cone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -693,7 +701,8 @@ contains
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(cone,vty,& if (info == psb_success_) call psb_map_X2Y(cone,vty,&
& czero,p%precv(level+1)%wrk%vx2l,& & czero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction') & a_err='Error during W-cycle restriction')
@ -704,7 +713,8 @@ contains
if (info == psb_success_) call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& if (info == psb_success_) call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,& & cone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -737,13 +747,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vty,cone,vy2l,& & vty,cone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vty,cone,vy2l,& & vty,cone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -760,7 +770,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info) & sweeps,work,wv,info)
else else
@ -836,7 +846,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -845,7 +856,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -854,13 +865,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vx2l,czero,vy2l,& & vx2l,czero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -890,7 +901,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(cone,vty,& call psb_map_X2Y(cone,vty,&
& czero,p%precv(level + 1)%wrk%vx2l,& & czero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -925,7 +937,8 @@ contains
! !
call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(cone,p%precv(level+1)%wrk%vy2l,&
& cone,vy2l,& & cone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -955,13 +968,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& vty,cone,vy2l,& & vty,cone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& vty,cone,vy2l,& & vty,cone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1014,7 +1027,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x

155
mlprec/impl/mld_cprecaply.f90
Unescape Escape View File

@ -323,6 +323,7 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
complex(psb_spk_), pointer :: work_(:) complex(psb_spk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_cprecaply' name='mld_cprecaply'
@ -358,6 +359,10 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
do_alloc_wrk = .not.allocated(prec%wrk)
if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
! !
! Number of levels > 1: apply the multilevel preconditioner ! Number of levels > 1: apply the multilevel preconditioner
@ -375,31 +380,29 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! !
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,&
! & wv => prec%precv(1)%wrk%wv)
! This is a kludge for handling the symmetrized GS case. if (allocated(prec%precv(1)%sm2a)) then
! Will need some rethinking. !
! ! This is a kludge for handling the symmetrized GS case.
twoside: block ! Will need some rethinking.
type(psb_c_vect_type) :: w1,w2 !
call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.)
call psb_geaxpby(cone,x,czero,w1,desc_data,info) call psb_geaxpby(cone,x,czero,w1,desc_data,info)
select case(trans_) select case(trans_)
case ('N') case ('N')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm%apply(cone,w1,czero,w2,desc_data,trans_,& call prec%precv(1)%sm%apply(cone,w1,czero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(cone,w2,czero,w1,desc_data,trans_,& call prec%precv(1)%sm2a%apply(cone,w2,czero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case('T','C') case('T','C')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm2a%apply(cone,w1,czero,w2,desc_data,trans_,& call prec%precv(1)%sm2a%apply(cone,w1,czero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm%apply(cone,w2,czero,w1,desc_data,trans_,& call prec%precv(1)%sm%apply(cone,w2,czero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case default case default
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
@ -407,14 +410,12 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end select end select
call psb_geaxpby(cone,w1,czero,y,desc_data,info) call psb_geaxpby(cone,w1,czero,y,desc_data,info)
call psb_gefree(w1,desc_data,info) else
call psb_gefree(w2,desc_data,info) call prec%precv(1)%sm%apply(cone,x,czero,y,desc_data,trans_,&
end block twoside & nswps,work_,wv,info)
end if
else end associate
call prec%precv(1)%sm%apply(cone,x,czero,y,desc_data,trans_,&
& nswps,work_,info)
end if
else else
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
@ -426,6 +427,7 @@ subroutine mld_cprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(y,desc_data,info,data=psb_comm_mov_) call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (do_alloc_wrk) call prec%free_wrk(info)
if (present(work)) then if (present(work)) then
else else
@ -459,10 +461,10 @@ subroutine mld_cprecaply1_vect(prec,x,desc_data,info,trans,work)
! Local variables ! Local variables
character :: trans_ character :: trans_
type(psb_c_vect_type) :: ww
complex(psb_spk_), pointer :: work_(:) complex(psb_spk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_cprecaply' name='mld_cprecaply'
@ -498,63 +500,68 @@ subroutine mld_cprecaply1_vect(prec,x,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
if (size(prec%precv) >1) then do_alloc_wrk = .not.allocated(prec%wrk)
! if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
! Number of levels > 1: apply the multilevel preconditioner
!
call mld_mlprec_aply(cone,prec,x,czero,ww,desc_data,trans_,work_,info)
if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_cmlprec_aply')
goto 9999
end if
else if (size(prec%precv) == 1) then associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv)
! call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
! Number of levels = 1: apply the base preconditioner if (size(prec%precv) >1) then
!
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
! !
! This is a kludge for handling the symmetrized GS case. ! Number of levels > 1: apply the multilevel preconditioner
! Will need some rethinking.
! !
select case(trans_) call mld_mlprec_aply(cone,prec,x,czero,ww,desc_data,trans_,work_,info)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm2a%apply(cone,ww,czero,x,desc_data,trans_,&
& ione, work_,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(cone,x,czero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm%apply(cone,ww,czero,x,desc_data,trans_,&
& ione, work_,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
else
call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,&
& nswps, work_,info)
if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info) if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info)
end if if(info /= psb_success_) then
else call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_cmlprec_aply')
goto 9999
end if
info = psb_err_from_subroutine_ai_ else if (size(prec%precv) == 1) then
call psb_errpush(info,name,a_err='Invalid size of precv',& !
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) ! Number of levels = 1: apply the base preconditioner
goto 9999 !
endif nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
!
! This is a kludge for handling the symmetrized GS case.
! Will need some rethinking.
!
select case(trans_)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(cone,ww,czero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(cone,x,czero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm%apply(cone,ww,czero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
if (info == 0) call psb_gefree(ww,desc_data,info) else
call prec%precv(1)%sm%apply(cone,x,czero,ww,desc_data,trans_,&
& nswps, work_,wv,info)
if (info == 0) call psb_geaxpby(cone,ww,czero,x,desc_data,info)
end if
else
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Invalid size of precv',&
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/))
goto 9999
endif
end associate
!!$ if (info == 0) call psb_gefree(ww,desc_data,info)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(x,desc_data,info,data=psb_comm_mov_) call psb_halo(x,desc_data,info,data=psb_comm_mov_)

4
mlprec/impl/mld_cprecset.F90
Unescape Escape View File

@ -432,7 +432,7 @@ subroutine mld_cprecsetsm(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsm'
info = psb_success_ info = psb_success_
@ -495,7 +495,7 @@ subroutine mld_cprecsetsv(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsv'
info = psb_success_ info = psb_success_

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

@ -255,7 +255,8 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(:))
! !
! At first iteration we must use the input BETA ! At first iteration we must use the input BETA
! !
@ -482,7 +483,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(done,& call psb_geaxpby(done,&
@ -494,12 +496,12 @@ contains
call p%precv(level)%sm%apply(done,& call p%precv(level)%sm%apply(done,&
& vy2l,dzero,vty,& & vy2l,dzero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(done,& call p%precv(level)%sm2a%apply(done,&
& vty,dzero,vy2l,& & vty,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
@ -507,7 +509,7 @@ contains
call p%precv(level)%sm%apply(done,& call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -519,7 +521,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(done,vx2l,& call psb_map_X2Y(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -538,7 +541,8 @@ contains
! !
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,& & done,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -605,7 +609,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level < nlev) then if (level < nlev) then
! !
! Apply the first smoother ! Apply the first smoother
@ -618,13 +623,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -651,7 +656,8 @@ contains
end if end if
call psb_map_X2Y(done,vty,& call psb_map_X2Y(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -661,7 +667,8 @@ contains
! Shortcut: just transfer x2l. ! Shortcut: just transfer x2l.
call psb_map_X2Y(done,vx2l,& call psb_map_X2Y(done,vx2l,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -676,7 +683,8 @@ contains
! !
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,& & done,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -693,7 +701,8 @@ contains
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(done,vty,& if (info == psb_success_) call psb_map_X2Y(done,vty,&
& dzero,p%precv(level+1)%wrk%vx2l,& & dzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction') & a_err='Error during W-cycle restriction')
@ -704,7 +713,8 @@ contains
if (info == psb_success_) call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& if (info == psb_success_) call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,& & done,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -737,13 +747,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vty,done,vy2l,& & vty,done,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vty,done,vy2l,& & vty,done,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -760,7 +770,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info) & sweeps,work,wv,info)
else else
@ -836,7 +846,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -845,7 +856,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -854,13 +865,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vx2l,dzero,vy2l,& & vx2l,dzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -890,7 +901,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(done,vty,& call psb_map_X2Y(done,vty,&
& dzero,p%precv(level + 1)%wrk%vx2l,& & dzero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -925,7 +937,8 @@ contains
! !
call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(done,p%precv(level+1)%wrk%vy2l,&
& done,vy2l,& & done,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -955,13 +968,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& vty,done,vy2l,& & vty,done,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& vty,done,vy2l,& & vty,done,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1014,7 +1027,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x

155
mlprec/impl/mld_dprecaply.f90
Unescape Escape View File

@ -323,6 +323,7 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
real(psb_dpk_), pointer :: work_(:) real(psb_dpk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_dprecaply' name='mld_dprecaply'
@ -358,6 +359,10 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
do_alloc_wrk = .not.allocated(prec%wrk)
if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
! !
! Number of levels > 1: apply the multilevel preconditioner ! Number of levels > 1: apply the multilevel preconditioner
@ -375,31 +380,29 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! !
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,&
! & wv => prec%precv(1)%wrk%wv)
! This is a kludge for handling the symmetrized GS case. if (allocated(prec%precv(1)%sm2a)) then
! Will need some rethinking. !
! ! This is a kludge for handling the symmetrized GS case.
twoside: block ! Will need some rethinking.
type(psb_d_vect_type) :: w1,w2 !
call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.)
call psb_geaxpby(done,x,dzero,w1,desc_data,info) call psb_geaxpby(done,x,dzero,w1,desc_data,info)
select case(trans_) select case(trans_)
case ('N') case ('N')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm%apply(done,w1,dzero,w2,desc_data,trans_,& call prec%precv(1)%sm%apply(done,w1,dzero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(done,w2,dzero,w1,desc_data,trans_,& call prec%precv(1)%sm2a%apply(done,w2,dzero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case('T','C') case('T','C')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm2a%apply(done,w1,dzero,w2,desc_data,trans_,& call prec%precv(1)%sm2a%apply(done,w1,dzero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm%apply(done,w2,dzero,w1,desc_data,trans_,& call prec%precv(1)%sm%apply(done,w2,dzero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case default case default
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
@ -407,14 +410,12 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end select end select
call psb_geaxpby(done,w1,dzero,y,desc_data,info) call psb_geaxpby(done,w1,dzero,y,desc_data,info)
call psb_gefree(w1,desc_data,info) else
call psb_gefree(w2,desc_data,info) call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,&
end block twoside & nswps,work_,wv,info)
end if
else end associate
call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,&
& nswps,work_,info)
end if
else else
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
@ -426,6 +427,7 @@ subroutine mld_dprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(y,desc_data,info,data=psb_comm_mov_) call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (do_alloc_wrk) call prec%free_wrk(info)
if (present(work)) then if (present(work)) then
else else
@ -459,10 +461,10 @@ subroutine mld_dprecaply1_vect(prec,x,desc_data,info,trans,work)
! Local variables ! Local variables
character :: trans_ character :: trans_
type(psb_d_vect_type) :: ww
real(psb_dpk_), pointer :: work_(:) real(psb_dpk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_dprecaply' name='mld_dprecaply'
@ -498,63 +500,68 @@ subroutine mld_dprecaply1_vect(prec,x,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
if (size(prec%precv) >1) then do_alloc_wrk = .not.allocated(prec%wrk)
! if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
! Number of levels > 1: apply the multilevel preconditioner
!
call mld_mlprec_aply(done,prec,x,dzero,ww,desc_data,trans_,work_,info)
if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply')
goto 9999
end if
else if (size(prec%precv) == 1) then associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv)
! call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
! Number of levels = 1: apply the base preconditioner if (size(prec%precv) >1) then
!
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
! !
! This is a kludge for handling the symmetrized GS case. ! Number of levels > 1: apply the multilevel preconditioner
! Will need some rethinking.
! !
select case(trans_) call mld_mlprec_aply(done,prec,x,dzero,ww,desc_data,trans_,work_,info)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm2a%apply(done,ww,dzero,x,desc_data,trans_,&
& ione, work_,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(done,x,dzero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm%apply(done,ww,dzero,x,desc_data,trans_,&
& ione, work_,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
else
call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,&
& nswps, work_,info)
if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info) if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info)
end if if(info /= psb_success_) then
else call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_dmlprec_aply')
goto 9999
end if
info = psb_err_from_subroutine_ai_ else if (size(prec%precv) == 1) then
call psb_errpush(info,name,a_err='Invalid size of precv',& !
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) ! Number of levels = 1: apply the base preconditioner
goto 9999 !
endif nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
!
! This is a kludge for handling the symmetrized GS case.
! Will need some rethinking.
!
select case(trans_)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(done,ww,dzero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(done,x,dzero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm%apply(done,ww,dzero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
if (info == 0) call psb_gefree(ww,desc_data,info) else
call prec%precv(1)%sm%apply(done,x,dzero,ww,desc_data,trans_,&
& nswps, work_,wv,info)
if (info == 0) call psb_geaxpby(done,ww,dzero,x,desc_data,info)
end if
else
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Invalid size of precv',&
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/))
goto 9999
endif
end associate
!!$ if (info == 0) call psb_gefree(ww,desc_data,info)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(x,desc_data,info,data=psb_comm_mov_) call psb_halo(x,desc_data,info,data=psb_comm_mov_)

4
mlprec/impl/mld_dprecset.F90
Unescape Escape View File

@ -465,7 +465,7 @@ subroutine mld_dprecsetsm(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsm'
info = psb_success_ info = psb_success_
@ -528,7 +528,7 @@ subroutine mld_dprecsetsv(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsv'
info = psb_success_ info = psb_success_

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

@ -255,7 +255,8 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(:))
! !
! At first iteration we must use the input BETA ! At first iteration we must use the input BETA
! !
@ -482,7 +483,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(sone,& call psb_geaxpby(sone,&
@ -494,12 +496,12 @@ contains
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& vy2l,szero,vty,& & vy2l,szero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(sone,& call p%precv(level)%sm2a%apply(sone,&
& vty,szero,vy2l,& & vty,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
@ -507,7 +509,7 @@ contains
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -519,7 +521,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(sone,vx2l,& call psb_map_X2Y(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,& & szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -538,7 +541,8 @@ contains
! !
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,& & sone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -605,7 +609,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level < nlev) then if (level < nlev) then
! !
! Apply the first smoother ! Apply the first smoother
@ -618,13 +623,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -651,7 +656,8 @@ contains
end if end if
call psb_map_X2Y(sone,vty,& call psb_map_X2Y(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,& & szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -661,7 +667,8 @@ contains
! Shortcut: just transfer x2l. ! Shortcut: just transfer x2l.
call psb_map_X2Y(sone,vx2l,& call psb_map_X2Y(sone,vx2l,&
& szero,p%precv(level+1)%wrk%vx2l,& & szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -676,7 +683,8 @@ contains
! !
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,& & sone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -693,7 +701,8 @@ contains
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(sone,vty,& if (info == psb_success_) call psb_map_X2Y(sone,vty,&
& szero,p%precv(level+1)%wrk%vx2l,& & szero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction') & a_err='Error during W-cycle restriction')
@ -704,7 +713,8 @@ contains
if (info == psb_success_) call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& if (info == psb_success_) call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,& & sone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -737,13 +747,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vty,sone,vy2l,& & vty,sone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vty,sone,vy2l,& & vty,sone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -760,7 +770,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info) & sweeps,work,wv,info)
else else
@ -836,7 +846,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -845,7 +856,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -854,13 +865,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vx2l,szero,vy2l,& & vx2l,szero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -890,7 +901,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(sone,vty,& call psb_map_X2Y(sone,vty,&
& szero,p%precv(level + 1)%wrk%vx2l,& & szero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -925,7 +937,8 @@ contains
! !
call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(sone,p%precv(level+1)%wrk%vy2l,&
& sone,vy2l,& & sone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -955,13 +968,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& vty,sone,vy2l,& & vty,sone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& vty,sone,vy2l,& & vty,sone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1014,7 +1027,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x

155
mlprec/impl/mld_sprecaply.f90
Unescape Escape View File

@ -323,6 +323,7 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
real(psb_spk_), pointer :: work_(:) real(psb_spk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_sprecaply' name='mld_sprecaply'
@ -358,6 +359,10 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
do_alloc_wrk = .not.allocated(prec%wrk)
if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
! !
! Number of levels > 1: apply the multilevel preconditioner ! Number of levels > 1: apply the multilevel preconditioner
@ -375,31 +380,29 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! !
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,&
! & wv => prec%precv(1)%wrk%wv)
! This is a kludge for handling the symmetrized GS case. if (allocated(prec%precv(1)%sm2a)) then
! Will need some rethinking. !
! ! This is a kludge for handling the symmetrized GS case.
twoside: block ! Will need some rethinking.
type(psb_s_vect_type) :: w1,w2 !
call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.)
call psb_geaxpby(sone,x,szero,w1,desc_data,info) call psb_geaxpby(sone,x,szero,w1,desc_data,info)
select case(trans_) select case(trans_)
case ('N') case ('N')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm%apply(sone,w1,szero,w2,desc_data,trans_,& call prec%precv(1)%sm%apply(sone,w1,szero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(sone,w2,szero,w1,desc_data,trans_,& call prec%precv(1)%sm2a%apply(sone,w2,szero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case('T','C') case('T','C')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm2a%apply(sone,w1,szero,w2,desc_data,trans_,& call prec%precv(1)%sm2a%apply(sone,w1,szero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm%apply(sone,w2,szero,w1,desc_data,trans_,& call prec%precv(1)%sm%apply(sone,w2,szero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case default case default
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
@ -407,14 +410,12 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end select end select
call psb_geaxpby(sone,w1,szero,y,desc_data,info) call psb_geaxpby(sone,w1,szero,y,desc_data,info)
call psb_gefree(w1,desc_data,info) else
call psb_gefree(w2,desc_data,info) call prec%precv(1)%sm%apply(sone,x,szero,y,desc_data,trans_,&
end block twoside & nswps,work_,wv,info)
end if
else end associate
call prec%precv(1)%sm%apply(sone,x,szero,y,desc_data,trans_,&
& nswps,work_,info)
end if
else else
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
@ -426,6 +427,7 @@ subroutine mld_sprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(y,desc_data,info,data=psb_comm_mov_) call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (do_alloc_wrk) call prec%free_wrk(info)
if (present(work)) then if (present(work)) then
else else
@ -459,10 +461,10 @@ subroutine mld_sprecaply1_vect(prec,x,desc_data,info,trans,work)
! Local variables ! Local variables
character :: trans_ character :: trans_
type(psb_s_vect_type) :: ww
real(psb_spk_), pointer :: work_(:) real(psb_spk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_sprecaply' name='mld_sprecaply'
@ -498,63 +500,68 @@ subroutine mld_sprecaply1_vect(prec,x,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
if (size(prec%precv) >1) then do_alloc_wrk = .not.allocated(prec%wrk)
! if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
! Number of levels > 1: apply the multilevel preconditioner
!
call mld_mlprec_aply(sone,prec,x,szero,ww,desc_data,trans_,work_,info)
if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_smlprec_aply')
goto 9999
end if
else if (size(prec%precv) == 1) then associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv)
! call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
! Number of levels = 1: apply the base preconditioner if (size(prec%precv) >1) then
!
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
! !
! This is a kludge for handling the symmetrized GS case. ! Number of levels > 1: apply the multilevel preconditioner
! Will need some rethinking.
! !
select case(trans_) call mld_mlprec_aply(sone,prec,x,szero,ww,desc_data,trans_,work_,info)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm2a%apply(sone,ww,szero,x,desc_data,trans_,&
& ione, work_,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(sone,x,szero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm%apply(sone,ww,szero,x,desc_data,trans_,&
& ione, work_,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
else
call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,&
& nswps, work_,info)
if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info) if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info)
end if if(info /= psb_success_) then
else call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_smlprec_aply')
goto 9999
end if
info = psb_err_from_subroutine_ai_ else if (size(prec%precv) == 1) then
call psb_errpush(info,name,a_err='Invalid size of precv',& !
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) ! Number of levels = 1: apply the base preconditioner
goto 9999 !
endif nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
!
! This is a kludge for handling the symmetrized GS case.
! Will need some rethinking.
!
select case(trans_)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(sone,ww,szero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(sone,x,szero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm%apply(sone,ww,szero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
if (info == 0) call psb_gefree(ww,desc_data,info) else
call prec%precv(1)%sm%apply(sone,x,szero,ww,desc_data,trans_,&
& nswps, work_,wv,info)
if (info == 0) call psb_geaxpby(sone,ww,szero,x,desc_data,info)
end if
else
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Invalid size of precv',&
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/))
goto 9999
endif
end associate
!!$ if (info == 0) call psb_gefree(ww,desc_data,info)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(x,desc_data,info,data=psb_comm_mov_) call psb_halo(x,desc_data,info,data=psb_comm_mov_)

4
mlprec/impl/mld_sprecset.F90
Unescape Escape View File

@ -432,7 +432,7 @@ subroutine mld_sprecsetsm(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsm'
info = psb_success_ info = psb_success_
@ -495,7 +495,7 @@ subroutine mld_sprecsetsv(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsv'
info = psb_success_ info = psb_success_

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

@ -255,7 +255,8 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv(:))
! !
! At first iteration we must use the input BETA ! At first iteration we must use the input BETA
! !
@ -482,7 +483,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (allocated(p%precv(level)%sm2a)) then if (allocated(p%precv(level)%sm2a)) then
call psb_geaxpby(zone,& call psb_geaxpby(zone,&
@ -494,12 +496,12 @@ contains
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& vy2l,zzero,vty,& & vy2l,zzero,vty,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
call p%precv(level)%sm2a%apply(zone,& call p%precv(level)%sm2a%apply(zone,&
& vty,zzero,vy2l,& & vty,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& ione,work,info,init='Z') & ione,work,wv,info,init='Z')
end do end do
else else
@ -507,7 +509,7 @@ contains
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -519,7 +521,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(zone,vx2l,& call psb_map_X2Y(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -538,7 +541,8 @@ contains
! !
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,& & zone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -605,7 +609,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level < nlev) then if (level < nlev) then
! !
! Apply the first smoother ! Apply the first smoother
@ -618,13 +623,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -651,7 +656,8 @@ contains
end if end if
call psb_map_X2Y(zone,vty,& call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -661,7 +667,8 @@ contains
! Shortcut: just transfer x2l. ! Shortcut: just transfer x2l.
call psb_map_X2Y(zone,vx2l,& call psb_map_X2Y(zone,vx2l,&
& zzero,p%precv(level+1)%wrk%vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -676,7 +683,8 @@ contains
! !
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,& & zone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation') & a_err='Error during prolongation')
@ -693,7 +701,8 @@ contains
& base_desc,info,work=work,trans=trans) & base_desc,info,work=work,trans=trans)
if (info == psb_success_) call psb_map_X2Y(zone,vty,& if (info == psb_success_) call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level+1)%wrk%vx2l,& & zzero,p%precv(level+1)%wrk%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during W-cycle restriction') & a_err='Error during W-cycle restriction')
@ -704,7 +713,8 @@ contains
if (info == psb_success_) call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& if (info == psb_success_) call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,& & zone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -737,13 +747,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vty,zone,vy2l,& & vty,zone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vty,zone,vy2l,& & vty,zone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -760,7 +770,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info) & sweeps,work,wv,info)
else else
@ -836,7 +846,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
if (level == nlev) then if (level == nlev) then
! !
! Apply smoother ! Apply smoother
@ -845,7 +856,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -854,13 +865,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vx2l,zzero,vy2l,& & vx2l,zzero,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -890,7 +901,8 @@ contains
! Apply the restriction ! Apply the restriction
call psb_map_X2Y(zone,vty,& call psb_map_X2Y(zone,vty,&
& zzero,p%precv(level + 1)%wrk%vx2l,& & zzero,p%precv(level + 1)%wrk%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -925,7 +937,8 @@ contains
! !
call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,& call psb_map_Y2X(zone,p%precv(level+1)%wrk%vy2l,&
& zone,vy2l,& & zone,vy2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work,&
& vtx=wv(1),vty=wv(2))
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -955,13 +968,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& vty,zone,vy2l,& & vty,zone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_pre sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& vty,zone,vy2l,& & vty,zone,vy2l,&
& base_desc, trans,& & base_desc, trans,&
& sweeps,work,info,init='Z') & sweeps,work,wv,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -1014,7 +1027,8 @@ contains
associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,& associate(vx2l => p%precv(level)%wrk%vx2l,vy2l => p%precv(level)%wrk%vy2l,&
& vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,& & vtx => p%precv(level)%wrk%vtx,vty => p%precv(level)%wrk%vty,&
& base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc) & base_a => p%precv(level)%base_a, base_desc=>p%precv(level)%base_desc,&
& wv => p%precv(level)%wrk%wv)
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x

155
mlprec/impl/mld_zprecaply.f90
Unescape Escape View File

@ -323,6 +323,7 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
complex(psb_dpk_), pointer :: work_(:) complex(psb_dpk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_zprecaply' name='mld_zprecaply'
@ -358,6 +359,10 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
do_alloc_wrk = .not.allocated(prec%wrk)
if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
if (size(prec%precv) >1) then if (size(prec%precv) >1) then
! !
! Number of levels > 1: apply the multilevel preconditioner ! Number of levels > 1: apply the multilevel preconditioner
@ -375,31 +380,29 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! !
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post) nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then associate(w1 => prec%precv(1)%wrk%vx2l, w2 => prec%precv(1)%wrk%vy2l,&
! & wv => prec%precv(1)%wrk%wv)
! This is a kludge for handling the symmetrized GS case. if (allocated(prec%precv(1)%sm2a)) then
! Will need some rethinking. !
! ! This is a kludge for handling the symmetrized GS case.
twoside: block ! Will need some rethinking.
type(psb_z_vect_type) :: w1,w2 !
call psb_geasb(w1,desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(w2,desc_data,info,mold=x%v,scratch=.true.)
call psb_geaxpby(zone,x,zzero,w1,desc_data,info) call psb_geaxpby(zone,x,zzero,w1,desc_data,info)
select case(trans_) select case(trans_)
case ('N') case ('N')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm%apply(zone,w1,zzero,w2,desc_data,trans_,& call prec%precv(1)%sm%apply(zone,w1,zzero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(zone,w2,zzero,w1,desc_data,trans_,& call prec%precv(1)%sm2a%apply(zone,w2,zzero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case('T','C') case('T','C')
do k=1, nswps do k=1, nswps
call prec%precv(1)%sm2a%apply(zone,w1,zzero,w2,desc_data,trans_,& call prec%precv(1)%sm2a%apply(zone,w1,zzero,w2,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
call prec%precv(1)%sm%apply(zone,w2,zzero,w1,desc_data,trans_,& call prec%precv(1)%sm%apply(zone,w2,zzero,w1,desc_data,trans_,&
& ione, work_,info) & ione, work_,wv,info)
end do end do
case default case default
info = psb_err_from_subroutine_ info = psb_err_from_subroutine_
@ -407,14 +410,12 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
goto 9999 goto 9999
end select end select
call psb_geaxpby(zone,w1,zzero,y,desc_data,info) call psb_geaxpby(zone,w1,zzero,y,desc_data,info)
call psb_gefree(w1,desc_data,info) else
call psb_gefree(w2,desc_data,info) call prec%precv(1)%sm%apply(zone,x,zzero,y,desc_data,trans_,&
end block twoside & nswps,work_,wv,info)
end if
else end associate
call prec%precv(1)%sm%apply(zone,x,zzero,y,desc_data,trans_,&
& nswps,work_,info)
end if
else else
info = psb_err_from_subroutine_ai_ info = psb_err_from_subroutine_ai_
@ -426,6 +427,7 @@ subroutine mld_zprecaply2_vect(prec,x,y,desc_data,info,trans,work)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(y,desc_data,info,data=psb_comm_mov_) call psb_halo(y,desc_data,info,data=psb_comm_mov_)
if (do_alloc_wrk) call prec%free_wrk(info)
if (present(work)) then if (present(work)) then
else else
@ -459,10 +461,10 @@ subroutine mld_zprecaply1_vect(prec,x,desc_data,info,trans,work)
! Local variables ! Local variables
character :: trans_ character :: trans_
type(psb_z_vect_type) :: ww
complex(psb_dpk_), pointer :: work_(:) complex(psb_dpk_), pointer :: work_(:)
integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: err_act,iwsz, k, nswps integer(psb_ipk_) :: err_act,iwsz, k, nswps
logical :: do_alloc_wrk
character(len=20) :: name character(len=20) :: name
name='mld_zprecaply' name='mld_zprecaply'
@ -498,63 +500,68 @@ subroutine mld_zprecaply1_vect(prec,x,desc_data,info,trans,work)
call psb_errpush(info,name) call psb_errpush(info,name)
goto 9999 goto 9999
end if end if
call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
if (size(prec%precv) >1) then do_alloc_wrk = .not.allocated(prec%wrk)
! if (do_alloc_wrk) call prec%allocate_wrk(info,vmold=x%v)
! Number of levels > 1: apply the multilevel preconditioner
!
call mld_mlprec_aply(zone,prec,x,zzero,ww,desc_data,trans_,work_,info)
if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info)
if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_zmlprec_aply')
goto 9999
end if
else if (size(prec%precv) == 1) then associate(ww => prec%precv(1)%wrk%vtx, wv => prec%precv(1)%wrk%wv)
! call psb_geasb(ww,desc_data,info,mold=x%v,scratch=.true.)
! Number of levels = 1: apply the base preconditioner if (size(prec%precv) >1) then
!
nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
! !
! This is a kludge for handling the symmetrized GS case. ! Number of levels > 1: apply the multilevel preconditioner
! Will need some rethinking.
! !
select case(trans_) call mld_mlprec_aply(zone,prec,x,zzero,ww,desc_data,trans_,work_,info)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm2a%apply(zone,ww,zzero,x,desc_data,trans_,&
& ione, work_,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(zone,x,zzero,ww,desc_data,trans_,&
& ione, work_,info)
call prec%precv(1)%sm%apply(zone,ww,zzero,x,desc_data,trans_,&
& ione, work_,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
else
call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,&
& nswps, work_,info)
if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info) if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info)
end if if(info /= psb_success_) then
else call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_zmlprec_aply')
goto 9999
end if
info = psb_err_from_subroutine_ai_ else if (size(prec%precv) == 1) then
call psb_errpush(info,name,a_err='Invalid size of precv',& !
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/)) ! Number of levels = 1: apply the base preconditioner
goto 9999 !
endif nswps = max(prec%precv(1)%parms%sweeps_pre,prec%precv(1)%parms%sweeps_post)
if (allocated(prec%precv(1)%sm2a)) then
!
! This is a kludge for handling the symmetrized GS case.
! Will need some rethinking.
!
select case(trans_)
case ('N')
do k=1, nswps
call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm2a%apply(zone,ww,zzero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case('T','C')
do k=1, nswps
call prec%precv(1)%sm2a%apply(zone,x,zzero,ww,desc_data,trans_,&
& ione, work_,wv,info)
call prec%precv(1)%sm%apply(zone,ww,zzero,x,desc_data,trans_,&
& ione, work_,wv,info)
end do
case default
info = psb_err_from_subroutine_
call psb_errpush(info,name,a_err='Invalid trans')
goto 9999
end select
if (info == 0) call psb_gefree(ww,desc_data,info) else
call prec%precv(1)%sm%apply(zone,x,zzero,ww,desc_data,trans_,&
& nswps, work_,wv,info)
if (info == 0) call psb_geaxpby(zone,ww,zzero,x,desc_data,info)
end if
else
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='Invalid size of precv',&
& i_Err=(/ione*size(prec%precv),izero,izero,izero,izero/))
goto 9999
endif
end associate
!!$ if (info == 0) call psb_gefree(ww,desc_data,info)
! If the original distribution has an overlap we should fix that. ! If the original distribution has an overlap we should fix that.
call psb_halo(x,desc_data,info,data=psb_comm_mov_) call psb_halo(x,desc_data,info,data=psb_comm_mov_)

4
mlprec/impl/mld_zprecset.F90
Unescape Escape View File

@ -465,7 +465,7 @@ subroutine mld_zprecsetsm(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsm'
info = psb_success_ info = psb_success_
@ -528,7 +528,7 @@ subroutine mld_zprecsetsv(p,val,info,ilev,ilmax,pos)
! Local variables ! Local variables
integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_ integer(psb_ipk_) :: ilev_, nlev_, ilmin_, ilmax_
character(len=*), parameter :: name='mld_precseti' character(len=*), parameter :: name='mld_precsetsv'
info = psb_success_ info = psb_success_

167
mlprec/impl/smoother/mld_c_as_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply_vect use mld_c_as_smoother, mld_protect_nam => mld_c_as_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_spk_), pointer :: aux(:) complex(psb_spk_), pointer :: aux(:)
@ -125,88 +125,95 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 3) then
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 goto 9999
endif end if
associate(tx => wv(1), ty => wv(2), ww => wv(3))
do i=1, sweeps-1
!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info) if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_) select case (init_)
case('Z')
if (info /= psb_success_) exit call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Z')
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y') case('Y')
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
if (info /= psb_success_) exit if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(cone,ww,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info) if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name,& & a_err='Error in sub_aply Jacobi Sweeps = 1')
& a_err='subsolve with Jacobi sweeps > 1') goto 9999
goto 9999 endif
end if
do i=1, sweeps-1
! !
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! ! block diagonal part and the remaining part of the local matrix
call psb_geaxpby(alpha,ty,beta,y,desc_data,info) ! and Y(j) is the approximate solution at sweep j.
!
if (info == 0) call psb_geaxpby(cone,tx,czero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-cone,sm%nd,ty,cone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(cone,ww,czero,ty,sm%desc_data,trans_,aux,info,init='Y')
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_
@ -220,9 +227,9 @@ subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (.not.(4*isz <= size(work))) then if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info) deallocate(aux,stat=info)
endif endif
if (info ==0) call ww%free(info) !!$ if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info) !!$ if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info) !!$ if (info ==0) call ty%free(info)
if (info /= 0) then if (info /= 0) then
info = psb_err_alloc_dealloc_ info = psb_err_alloc_dealloc_
call psb_errpush(info,name) call psb_errpush(info,name)

4
mlprec/impl/smoother/mld_c_base_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_c_base_smoother_mod, mld_protect_name => mld_c_base_smoother_apply_vect use mld_c_base_smoother_mod, mld_protect_name => mld_c_base_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='c_base_smoother_apply' character(len=20) :: name='c_base_smoother_apply'

133
mlprec/impl/smoother/mld_c_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply_vect use mld_c_jac_smoother, mld_protect_name => mld_c_jac_smoother_apply_vect
@ -49,10 +49,10 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_c_vect_type) :: tx, ty type(psb_c_vect_type) :: tx, ty
@ -122,77 +122,86 @@ subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 2) then
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) info = psb_err_internal_error_
write(0,*) 'Size (WV) : ',size(wv)
call psb_errpush(info,name,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one AXPBY and one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
select case (init_)
case('Z')
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one AXPBY and one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(cone,x,czero,tx,desc_data,info) select case (init_)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_) case('Z')
call sm%sv%apply(cone,x,czero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,y,czero,ty,desc_data,info)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_geaxpby(cone,initu,czero,ty,desc_data,info)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info /= psb_success_) exit do i=1, sweeps-1
!
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! block diagonal part and the remaining part of the local matrix
! and Y(j) is the approximate solution at sweep j.
!
call psb_geaxpby(cone,x,czero,tx,desc_data,info)
call psb_spmm(-cone,sm%nd,ty,cone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y') if (info /= psb_success_) exit
if (info /= psb_success_) exit call sm%sv%apply(cone,tx,czero,ty,desc_data,trans_,aux,info,init='Y')
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if (info /= psb_success_) exit
end do
if (info /= psb_success_) then if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
call tx%free(info) if (info /= psb_success_) then
if (info == psb_success_) call ty%free(info) info=psb_err_internal_error_
if (info /= psb_success_) then call psb_errpush(info,name,&
info=psb_err_internal_error_ & a_err='subsolve with Jacobi sweeps > 1')
call psb_errpush(info,name,& goto 9999
& a_err='final cleanup with Jacobi sweeps > 1') end if
goto 9999
end if
!!$ call tx%free(info)
!!$ if (info == psb_success_) call ty%free(info)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='final cleanup with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_

167
mlprec/impl/smoother/mld_d_as_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply_vect use mld_d_as_smoother, mld_protect_nam => mld_d_as_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i integer(psb_ipk_) :: n_row,n_col, nrow_d, i
real(psb_dpk_), pointer :: aux(:) real(psb_dpk_), pointer :: aux(:)
@ -125,88 +125,95 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 3) then
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 goto 9999
endif end if
associate(tx => wv(1), ty => wv(2), ww => wv(3))
do i=1, sweeps-1
!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info) if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_) select case (init_)
case('Z')
if (info /= psb_success_) exit call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Z')
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y') case('Y')
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
if (info /= psb_success_) exit if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(done,ww,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info) if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name,& & a_err='Error in sub_aply Jacobi Sweeps = 1')
& a_err='subsolve with Jacobi sweeps > 1') goto 9999
goto 9999 endif
end if
do i=1, sweeps-1
! !
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! ! block diagonal part and the remaining part of the local matrix
call psb_geaxpby(alpha,ty,beta,y,desc_data,info) ! and Y(j) is the approximate solution at sweep j.
!
if (info == 0) call psb_geaxpby(done,tx,dzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-done,sm%nd,ty,done,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(done,ww,dzero,ty,sm%desc_data,trans_,aux,info,init='Y')
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_
@ -220,9 +227,9 @@ subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (.not.(4*isz <= size(work))) then if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info) deallocate(aux,stat=info)
endif endif
if (info ==0) call ww%free(info) !!$ if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info) !!$ if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info) !!$ if (info ==0) call ty%free(info)
if (info /= 0) then if (info /= 0) then
info = psb_err_alloc_dealloc_ info = psb_err_alloc_dealloc_
call psb_errpush(info,name) call psb_errpush(info,name)

4
mlprec/impl/smoother/mld_d_base_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_d_base_smoother_mod, mld_protect_name => mld_d_base_smoother_apply_vect use mld_d_base_smoother_mod, mld_protect_name => mld_d_base_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='d_base_smoother_apply' character(len=20) :: name='d_base_smoother_apply'

133
mlprec/impl/smoother/mld_d_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply_vect use mld_d_jac_smoother, mld_protect_name => mld_d_jac_smoother_apply_vect
@ -49,10 +49,10 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_d_vect_type) :: tx, ty type(psb_d_vect_type) :: tx, ty
@ -122,77 +122,86 @@ subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 2) then
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) info = psb_err_internal_error_
write(0,*) 'Size (WV) : ',size(wv)
call psb_errpush(info,name,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one AXPBY and one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
select case (init_)
case('Z')
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one AXPBY and one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(done,x,dzero,tx,desc_data,info) select case (init_)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_) case('Z')
call sm%sv%apply(done,x,dzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,y,dzero,ty,desc_data,info)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_geaxpby(done,initu,dzero,ty,desc_data,info)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info /= psb_success_) exit do i=1, sweeps-1
!
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! block diagonal part and the remaining part of the local matrix
! and Y(j) is the approximate solution at sweep j.
!
call psb_geaxpby(done,x,dzero,tx,desc_data,info)
call psb_spmm(-done,sm%nd,ty,done,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y') if (info /= psb_success_) exit
if (info /= psb_success_) exit call sm%sv%apply(done,tx,dzero,ty,desc_data,trans_,aux,info,init='Y')
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if (info /= psb_success_) exit
end do
if (info /= psb_success_) then if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
call tx%free(info) if (info /= psb_success_) then
if (info == psb_success_) call ty%free(info) info=psb_err_internal_error_
if (info /= psb_success_) then call psb_errpush(info,name,&
info=psb_err_internal_error_ & a_err='subsolve with Jacobi sweeps > 1')
call psb_errpush(info,name,& goto 9999
& a_err='final cleanup with Jacobi sweeps > 1') end if
goto 9999
end if
!!$ call tx%free(info)
!!$ if (info == psb_success_) call ty%free(info)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='final cleanup with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_

167
mlprec/impl/smoother/mld_s_as_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply_vect use mld_s_as_smoother, mld_protect_nam => mld_s_as_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i integer(psb_ipk_) :: n_row,n_col, nrow_d, i
real(psb_spk_), pointer :: aux(:) real(psb_spk_), pointer :: aux(:)
@ -125,88 +125,95 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 3) then
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 goto 9999
endif end if
associate(tx => wv(1), ty => wv(2), ww => wv(3))
do i=1, sweeps-1
!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info) if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_) select case (init_)
case('Z')
if (info /= psb_success_) exit call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Z')
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y') case('Y')
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
if (info /= psb_success_) exit if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(sone,ww,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info) if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name,& & a_err='Error in sub_aply Jacobi Sweeps = 1')
& a_err='subsolve with Jacobi sweeps > 1') goto 9999
goto 9999 endif
end if
do i=1, sweeps-1
! !
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! ! block diagonal part and the remaining part of the local matrix
call psb_geaxpby(alpha,ty,beta,y,desc_data,info) ! and Y(j) is the approximate solution at sweep j.
!
if (info == 0) call psb_geaxpby(sone,tx,szero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-sone,sm%nd,ty,sone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(sone,ww,szero,ty,sm%desc_data,trans_,aux,info,init='Y')
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_
@ -220,9 +227,9 @@ subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (.not.(4*isz <= size(work))) then if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info) deallocate(aux,stat=info)
endif endif
if (info ==0) call ww%free(info) !!$ if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info) !!$ if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info) !!$ if (info ==0) call ty%free(info)
if (info /= 0) then if (info /= 0) then
info = psb_err_alloc_dealloc_ info = psb_err_alloc_dealloc_
call psb_errpush(info,name) call psb_errpush(info,name)

4
mlprec/impl/smoother/mld_s_base_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_s_base_smoother_mod, mld_protect_name => mld_s_base_smoother_apply_vect use mld_s_base_smoother_mod, mld_protect_name => mld_s_base_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='s_base_smoother_apply' character(len=20) :: name='s_base_smoother_apply'

133
mlprec/impl/smoother/mld_s_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply_vect use mld_s_jac_smoother, mld_protect_name => mld_s_jac_smoother_apply_vect
@ -49,10 +49,10 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_s_vect_type) :: tx, ty type(psb_s_vect_type) :: tx, ty
@ -122,77 +122,86 @@ subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 2) then
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) info = psb_err_internal_error_
write(0,*) 'Size (WV) : ',size(wv)
call psb_errpush(info,name,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one AXPBY and one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
select case (init_)
case('Z')
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one AXPBY and one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(sone,x,szero,tx,desc_data,info) select case (init_)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_) case('Z')
call sm%sv%apply(sone,x,szero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,y,szero,ty,desc_data,info)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_geaxpby(sone,initu,szero,ty,desc_data,info)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info /= psb_success_) exit do i=1, sweeps-1
!
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! block diagonal part and the remaining part of the local matrix
! and Y(j) is the approximate solution at sweep j.
!
call psb_geaxpby(sone,x,szero,tx,desc_data,info)
call psb_spmm(-sone,sm%nd,ty,sone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y') if (info /= psb_success_) exit
if (info /= psb_success_) exit call sm%sv%apply(sone,tx,szero,ty,desc_data,trans_,aux,info,init='Y')
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if (info /= psb_success_) exit
end do
if (info /= psb_success_) then if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
call tx%free(info) if (info /= psb_success_) then
if (info == psb_success_) call ty%free(info) info=psb_err_internal_error_
if (info /= psb_success_) then call psb_errpush(info,name,&
info=psb_err_internal_error_ & a_err='subsolve with Jacobi sweeps > 1')
call psb_errpush(info,name,& goto 9999
& a_err='final cleanup with Jacobi sweeps > 1') end if
goto 9999
end if
!!$ call tx%free(info)
!!$ if (info == psb_success_) call ty%free(info)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='final cleanup with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_

167
mlprec/impl/smoother/mld_z_as_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply_vect use mld_z_as_smoother, mld_protect_nam => mld_z_as_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
integer(psb_ipk_) :: n_row,n_col, nrow_d, i integer(psb_ipk_) :: n_row,n_col, nrow_d, i
complex(psb_dpk_), pointer :: aux(:) complex(psb_dpk_), pointer :: aux(:)
@ -125,88 +125,95 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 3) then
call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
select case (init_)
case('Z')
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply') & a_err='invalid wv size in smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in sub_aply Jacobi Sweeps = 1')
goto 9999 goto 9999
endif end if
associate(tx => wv(1), ty => wv(2), ww => wv(3))
do i=1, sweeps-1
!!$ call psb_geasb(tx,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,sm%desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ww,sm%desc_data,info,mold=x%v,scratch=.true.)
! Need to zero tx because of the apply_restr call.
call tx%zero()
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
if (info == 0) call sm%apply_restr(tx,trans_,aux,info)
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info) if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_) select case (init_)
case('Z')
if (info /= psb_success_) exit call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Z')
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y') case('Y')
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
if (info /= psb_success_) exit if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
if (info == 0) call sm%apply_restr(ty,trans_,aux,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
call sm%sv%apply(zone,ww,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info == 0) call sm%apply_prol(ty,trans_,aux,info) if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then if (info /= psb_success_) then
info=psb_err_internal_error_ call psb_errpush(psb_err_internal_error_,name,&
call psb_errpush(info,name,& & a_err='Error in sub_aply Jacobi Sweeps = 1')
& a_err='subsolve with Jacobi sweeps > 1') goto 9999
goto 9999 endif
end if
do i=1, sweeps-1
! !
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx) ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! ! block diagonal part and the remaining part of the local matrix
call psb_geaxpby(alpha,ty,beta,y,desc_data,info) ! and Y(j) is the approximate solution at sweep j.
!
if (info == 0) call psb_geaxpby(zone,tx,zzero,ww,sm%desc_data,info)
if (info == 0) call psb_spmm(-zone,sm%nd,ty,zone,ww,sm%desc_data,info,&
& work=aux,trans=trans_)
if (info /= psb_success_) exit
call sm%sv%apply(zone,ww,zzero,ty,sm%desc_data,trans_,aux,info,init='Y')
if (info /= psb_success_) exit
if (info == 0) call sm%apply_prol(ty,trans_,aux,info)
end do
if (info /= psb_success_) then
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
!
! Compute y = beta*y + alpha*ty (ty == K^(-1)*tx)
!
call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_
@ -220,9 +227,9 @@ subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
if (.not.(4*isz <= size(work))) then if (.not.(4*isz <= size(work))) then
deallocate(aux,stat=info) deallocate(aux,stat=info)
endif endif
if (info ==0) call ww%free(info) !!$ if (info ==0) call ww%free(info)
if (info ==0) call tx%free(info) !!$ if (info ==0) call tx%free(info)
if (info ==0) call ty%free(info) !!$ if (info ==0) call ty%free(info)
if (info /= 0) then if (info /= 0) then
info = psb_err_alloc_dealloc_ info = psb_err_alloc_dealloc_
call psb_errpush(info,name) call psb_errpush(info,name)

4
mlprec/impl/smoother/mld_z_base_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_z_base_smoother_mod, mld_protect_name => mld_z_base_smoother_apply_vect use mld_z_base_smoother_mod, mld_protect_name => mld_z_base_smoother_apply_vect
implicit none implicit none
@ -48,10 +48,10 @@ subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: err_act integer(psb_ipk_) :: err_act
character(len=20) :: name='z_base_smoother_apply' character(len=20) :: name='z_base_smoother_apply'

133
mlprec/impl/smoother/mld_z_jac_smoother_apply_vect.f90
Unescape Escape View File

@ -36,7 +36,7 @@
! !
! !
subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
use psb_base_mod use psb_base_mod
use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply_vect use mld_z_jac_smoother, mld_protect_name => mld_z_jac_smoother_apply_vect
@ -49,10 +49,10 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
! !
integer(psb_ipk_) :: n_row,n_col integer(psb_ipk_) :: n_row,n_col
type(psb_z_vect_type) :: tx, ty type(psb_z_vect_type) :: tx, ty
@ -122,77 +122,86 @@ subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
! to compute an approximate solution of a linear system. ! to compute an approximate solution of a linear system.
! !
! !
call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.) if (size(wv) < 2) then
call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.) info = psb_err_internal_error_
write(0,*) 'Size (WV) : ',size(wv)
call psb_errpush(info,name,&
& a_err='invalid wv size in smoother_apply')
goto 9999
end if
associate(tx => wv(1), ty => wv(2))
!!$ call psb_geasb(tx,desc_data,info,mold=x%v,scratch=.true.)
!!$ call psb_geasb(ty,desc_data,info,mold=x%v,scratch=.true.)
!
! Unroll the first iteration and fold it inside SELECT CASE
! this will save one AXPBY and one SPMM when INIT=Z, and will be
! significant when sweeps=1 (a common case)
!
select case (init_)
case('Z')
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
do i=1, sweeps-1
! !
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! Unroll the first iteration and fold it inside SELECT CASE
! block diagonal part and the remaining part of the local matrix ! this will save one AXPBY and one SPMM when INIT=Z, and will be
! and Y(j) is the approximate solution at sweep j. ! significant when sweeps=1 (a common case)
! !
call psb_geaxpby(zone,x,zzero,tx,desc_data,info) select case (init_)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_) case('Z')
call sm%sv%apply(zone,x,zzero,ty,desc_data,trans_,aux,info,init='Z')
case('Y')
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,y,zzero,ty,desc_data,info)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
case('U')
if (.not.present(initu)) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='missing initu to smoother_apply')
goto 9999
end if
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_geaxpby(zone,initu,zzero,ty,desc_data,info)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
case default
call psb_errpush(psb_err_internal_error_,name,&
& a_err='wrong init to smoother_apply')
goto 9999
end select
if (info /= psb_success_) exit do i=1, sweeps-1
!
! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the
! block diagonal part and the remaining part of the local matrix
! and Y(j) is the approximate solution at sweep j.
!
call psb_geaxpby(zone,x,zzero,tx,desc_data,info)
call psb_spmm(-zone,sm%nd,ty,zone,tx,desc_data,info,work=aux,trans=trans_)
call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y') if (info /= psb_success_) exit
if (info /= psb_success_) exit call sm%sv%apply(zone,tx,zzero,ty,desc_data,trans_,aux,info,init='Y')
end do
if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info) if (info /= psb_success_) exit
end do
if (info /= psb_success_) then if (info == psb_success_) call psb_geaxpby(alpha,ty,beta,y,desc_data,info)
info=psb_err_internal_error_
call psb_errpush(info,name,&
& a_err='subsolve with Jacobi sweeps > 1')
goto 9999
end if
call tx%free(info) if (info /= psb_success_) then
if (info == psb_success_) call ty%free(info) info=psb_err_internal_error_
if (info /= psb_success_) then call psb_errpush(info,name,&
info=psb_err_internal_error_ & a_err='subsolve with Jacobi sweeps > 1')
call psb_errpush(info,name,& goto 9999
& a_err='final cleanup with Jacobi sweeps > 1') end if
goto 9999
end if
!!$ call tx%free(info)
!!$ if (info == psb_success_) call ty%free(info)
!!$ if (info /= psb_success_) then
!!$ info=psb_err_internal_error_
!!$ call psb_errpush(info,name,&
!!$ & a_err='final cleanup with Jacobi sweeps > 1')
!!$ goto 9999
!!$ end if
end associate
else else
info = psb_err_iarg_neg_ info = psb_err_iarg_neg_

4
mlprec/mld_c_as_smoother.f90
Unescape Escape View File

@ -181,7 +181,7 @@ module mld_c_as_smoother
interface interface
subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_c_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, & import :: psb_cspmat_type, psb_c_vect_type, psb_c_base_vect_type, &
& psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, & & psb_spk_, mld_c_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_ & psb_desc_type, psb_ipk_
@ -194,10 +194,10 @@ module mld_c_as_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_c_as_smoother_apply_vect end subroutine mld_c_as_smoother_apply_vect
end interface end interface

4
mlprec/mld_c_base_smoother_mod.f90
Unescape Escape View File

@ -158,7 +158,7 @@ module mld_c_base_smoother_mod
interface interface
subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_c_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, & import :: psb_desc_type, psb_cspmat_type, psb_c_base_sparse_mat, &
& psb_c_vect_type, psb_c_base_vect_type, psb_spk_, & & psb_c_vect_type, psb_c_base_vect_type, psb_spk_, &
& mld_c_base_smoother_type, psb_ipk_ & mld_c_base_smoother_type, psb_ipk_
@ -170,10 +170,10 @@ module mld_c_base_smoother_mod
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_c_base_smoother_apply_vect end subroutine mld_c_base_smoother_apply_vect
end interface end interface

4
mlprec/mld_c_jac_smoother.f90
Unescape Escape View File

@ -85,7 +85,7 @@ module mld_c_jac_smoother
interface interface
subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_c_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, & import :: psb_desc_type, mld_c_jac_smoother_type, psb_c_vect_type, psb_spk_, &
& psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,& & psb_cspmat_type, psb_c_base_sparse_mat, psb_c_base_vect_type,&
& psb_ipk_ & psb_ipk_
@ -98,10 +98,10 @@ module mld_c_jac_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_spk_),target, intent(inout) :: work(:) complex(psb_spk_),target, intent(inout) :: work(:)
type(psb_c_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_c_vect_type),intent(inout), optional :: initu type(psb_c_vect_type),intent(inout), optional :: initu
type(psb_c_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_c_jac_smoother_apply_vect end subroutine mld_c_jac_smoother_apply_vect
end interface end interface

17
mlprec/mld_c_onelev_mod.f90
Unescape Escape View File

@ -583,10 +583,9 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_c_base_vect_type), intent(in), optional :: vmold class(psb_c_base_vect_type), intent(in), optional :: vmold
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv, i
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,& call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
@ -599,7 +598,12 @@ contains
call psb_geasb(lv%wrk%vty,& call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(lv%wrk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(lv%wrk%wv(i),&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end do
end subroutine c_base_onelev_allocate_wrk end subroutine c_base_onelev_allocate_wrk
@ -609,12 +613,17 @@ contains
class(mld_c_onelev_type), target, intent(inout) :: lv class(mld_c_onelev_type), target, intent(inout) :: lv
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv,i
info = psb_success_ info = psb_success_
call lv%wrk%vx2l%free(info) call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info) call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info) call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info) call lv%wrk%vty%free(info)
nwv = size(lv%wrk%wv)
do i=1,nwv
call lv%wrk%wv(i)%free(info)
end do
end subroutine c_base_onelev_free_wrk end subroutine c_base_onelev_free_wrk
end module mld_c_onelev_mod end module mld_c_onelev_mod

4
mlprec/mld_d_as_smoother.f90
Unescape Escape View File

@ -181,7 +181,7 @@ module mld_d_as_smoother
interface interface
subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_d_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, & import :: psb_dspmat_type, psb_d_vect_type, psb_d_base_vect_type, &
& psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, & & psb_dpk_, mld_d_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_ & psb_desc_type, psb_ipk_
@ -194,10 +194,10 @@ module mld_d_as_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_d_as_smoother_apply_vect end subroutine mld_d_as_smoother_apply_vect
end interface end interface

4
mlprec/mld_d_base_smoother_mod.f90
Unescape Escape View File

@ -158,7 +158,7 @@ module mld_d_base_smoother_mod
interface interface
subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_d_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, & import :: psb_desc_type, psb_dspmat_type, psb_d_base_sparse_mat, &
& psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, & & psb_d_vect_type, psb_d_base_vect_type, psb_dpk_, &
& mld_d_base_smoother_type, psb_ipk_ & mld_d_base_smoother_type, psb_ipk_
@ -170,10 +170,10 @@ module mld_d_base_smoother_mod
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_d_base_smoother_apply_vect end subroutine mld_d_base_smoother_apply_vect
end interface end interface

4
mlprec/mld_d_jac_smoother.f90
Unescape Escape View File

@ -85,7 +85,7 @@ module mld_d_jac_smoother
interface interface
subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_d_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, & import :: psb_desc_type, mld_d_jac_smoother_type, psb_d_vect_type, psb_dpk_, &
& psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,& & psb_dspmat_type, psb_d_base_sparse_mat, psb_d_base_vect_type,&
& psb_ipk_ & psb_ipk_
@ -98,10 +98,10 @@ module mld_d_jac_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_dpk_),target, intent(inout) :: work(:) real(psb_dpk_),target, intent(inout) :: work(:)
type(psb_d_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_d_vect_type),intent(inout), optional :: initu type(psb_d_vect_type),intent(inout), optional :: initu
type(psb_d_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_d_jac_smoother_apply_vect end subroutine mld_d_jac_smoother_apply_vect
end interface end interface

17
mlprec/mld_d_onelev_mod.f90
Unescape Escape View File

@ -583,10 +583,9 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_d_base_vect_type), intent(in), optional :: vmold class(psb_d_base_vect_type), intent(in), optional :: vmold
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv, i
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,& call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
@ -599,7 +598,12 @@ contains
call psb_geasb(lv%wrk%vty,& call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(lv%wrk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(lv%wrk%wv(i),&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end do
end subroutine d_base_onelev_allocate_wrk end subroutine d_base_onelev_allocate_wrk
@ -609,12 +613,17 @@ contains
class(mld_d_onelev_type), target, intent(inout) :: lv class(mld_d_onelev_type), target, intent(inout) :: lv
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv,i
info = psb_success_ info = psb_success_
call lv%wrk%vx2l%free(info) call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info) call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info) call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info) call lv%wrk%vty%free(info)
nwv = size(lv%wrk%wv)
do i=1,nwv
call lv%wrk%wv(i)%free(info)
end do
end subroutine d_base_onelev_free_wrk end subroutine d_base_onelev_free_wrk
end module mld_d_onelev_mod end module mld_d_onelev_mod

4
mlprec/mld_s_as_smoother.f90
Unescape Escape View File

@ -181,7 +181,7 @@ module mld_s_as_smoother
interface interface
subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_s_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, & import :: psb_sspmat_type, psb_s_vect_type, psb_s_base_vect_type, &
& psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, & & psb_spk_, mld_s_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_ & psb_desc_type, psb_ipk_
@ -194,10 +194,10 @@ module mld_s_as_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_s_as_smoother_apply_vect end subroutine mld_s_as_smoother_apply_vect
end interface end interface

4
mlprec/mld_s_base_smoother_mod.f90
Unescape Escape View File

@ -158,7 +158,7 @@ module mld_s_base_smoother_mod
interface interface
subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_s_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, & import :: psb_desc_type, psb_sspmat_type, psb_s_base_sparse_mat, &
& psb_s_vect_type, psb_s_base_vect_type, psb_spk_, & & psb_s_vect_type, psb_s_base_vect_type, psb_spk_, &
& mld_s_base_smoother_type, psb_ipk_ & mld_s_base_smoother_type, psb_ipk_
@ -170,10 +170,10 @@ module mld_s_base_smoother_mod
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_s_base_smoother_apply_vect end subroutine mld_s_base_smoother_apply_vect
end interface end interface

4
mlprec/mld_s_jac_smoother.f90
Unescape Escape View File

@ -85,7 +85,7 @@ module mld_s_jac_smoother
interface interface
subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_s_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, & import :: psb_desc_type, mld_s_jac_smoother_type, psb_s_vect_type, psb_spk_, &
& psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,& & psb_sspmat_type, psb_s_base_sparse_mat, psb_s_base_vect_type,&
& psb_ipk_ & psb_ipk_
@ -98,10 +98,10 @@ module mld_s_jac_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
real(psb_spk_),target, intent(inout) :: work(:) real(psb_spk_),target, intent(inout) :: work(:)
type(psb_s_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_s_vect_type),intent(inout), optional :: initu type(psb_s_vect_type),intent(inout), optional :: initu
type(psb_s_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_s_jac_smoother_apply_vect end subroutine mld_s_jac_smoother_apply_vect
end interface end interface

17
mlprec/mld_s_onelev_mod.f90
Unescape Escape View File

@ -583,10 +583,9 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_s_base_vect_type), intent(in), optional :: vmold class(psb_s_base_vect_type), intent(in), optional :: vmold
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv, i
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,& call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
@ -599,7 +598,12 @@ contains
call psb_geasb(lv%wrk%vty,& call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(lv%wrk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(lv%wrk%wv(i),&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end do
end subroutine s_base_onelev_allocate_wrk end subroutine s_base_onelev_allocate_wrk
@ -609,12 +613,17 @@ contains
class(mld_s_onelev_type), target, intent(inout) :: lv class(mld_s_onelev_type), target, intent(inout) :: lv
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv,i
info = psb_success_ info = psb_success_
call lv%wrk%vx2l%free(info) call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info) call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info) call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info) call lv%wrk%vty%free(info)
nwv = size(lv%wrk%wv)
do i=1,nwv
call lv%wrk%wv(i)%free(info)
end do
end subroutine s_base_onelev_free_wrk end subroutine s_base_onelev_free_wrk
end module mld_s_onelev_mod end module mld_s_onelev_mod

4
mlprec/mld_z_as_smoother.f90
Unescape Escape View File

@ -181,7 +181,7 @@ module mld_z_as_smoother
interface interface
subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_z_as_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, & import :: psb_zspmat_type, psb_z_vect_type, psb_z_base_vect_type, &
& psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, & & psb_dpk_, mld_z_as_smoother_type, psb_long_int_k_, &
& psb_desc_type, psb_ipk_ & psb_desc_type, psb_ipk_
@ -194,10 +194,10 @@ module mld_z_as_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_z_as_smoother_apply_vect end subroutine mld_z_as_smoother_apply_vect
end interface end interface

4
mlprec/mld_z_base_smoother_mod.f90
Unescape Escape View File

@ -158,7 +158,7 @@ module mld_z_base_smoother_mod
interface interface
subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,& subroutine mld_z_base_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,&
& trans,sweeps,work,info,init,initu,wv) & trans,sweeps,work,wv,info,init,initu)
import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, & import :: psb_desc_type, psb_zspmat_type, psb_z_base_sparse_mat, &
& psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, & & psb_z_vect_type, psb_z_base_vect_type, psb_dpk_, &
& mld_z_base_smoother_type, psb_ipk_ & mld_z_base_smoother_type, psb_ipk_
@ -170,10 +170,10 @@ module mld_z_base_smoother_mod
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_z_base_smoother_apply_vect end subroutine mld_z_base_smoother_apply_vect
end interface end interface

4
mlprec/mld_z_jac_smoother.f90
Unescape Escape View File

@ -85,7 +85,7 @@ module mld_z_jac_smoother
interface interface
subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,& subroutine mld_z_jac_smoother_apply_vect(alpha,sm,x,beta,y,desc_data,trans,&
& sweeps,work,info,init,initu,wv) & sweeps,work,wv,info,init,initu)
import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, & import :: psb_desc_type, mld_z_jac_smoother_type, psb_z_vect_type, psb_dpk_, &
& psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,& & psb_zspmat_type, psb_z_base_sparse_mat, psb_z_base_vect_type,&
& psb_ipk_ & psb_ipk_
@ -98,10 +98,10 @@ module mld_z_jac_smoother
character(len=1),intent(in) :: trans character(len=1),intent(in) :: trans
integer(psb_ipk_), intent(in) :: sweeps integer(psb_ipk_), intent(in) :: sweeps
complex(psb_dpk_),target, intent(inout) :: work(:) complex(psb_dpk_),target, intent(inout) :: work(:)
type(psb_z_vect_type),intent(inout) :: wv(:)
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
character, intent(in), optional :: init character, intent(in), optional :: init
type(psb_z_vect_type),intent(inout), optional :: initu type(psb_z_vect_type),intent(inout), optional :: initu
type(psb_z_vect_type),intent(inout), optional :: wv(:)
end subroutine mld_z_jac_smoother_apply_vect end subroutine mld_z_jac_smoother_apply_vect
end interface end interface

17
mlprec/mld_z_onelev_mod.f90
Unescape Escape View File

@ -583,10 +583,9 @@ contains
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
class(psb_z_base_vect_type), intent(in), optional :: vmold class(psb_z_base_vect_type), intent(in), optional :: vmold
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv, i
info = psb_success_ info = psb_success_
nwv = lv%get_wrksz() nwv = lv%get_wrksz()
write(0,*) 'Debug allocate_wrk: ',nwv
call psb_geasb(lv%wrk%vx2l,& call psb_geasb(lv%wrk%vx2l,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
@ -599,7 +598,12 @@ contains
call psb_geasb(lv%wrk%vty,& call psb_geasb(lv%wrk%vty,&
& lv%base_desc,info,& & lv%base_desc,info,&
& scratch=.true.,mold=vmold) & scratch=.true.,mold=vmold)
allocate(lv%wrk%wv(nwv),stat=info)
do i=1,nwv
call psb_geasb(lv%wrk%wv(i),&
& lv%base_desc,info,&
& scratch=.true.,mold=vmold)
end do
end subroutine z_base_onelev_allocate_wrk end subroutine z_base_onelev_allocate_wrk
@ -609,12 +613,17 @@ contains
class(mld_z_onelev_type), target, intent(inout) :: lv class(mld_z_onelev_type), target, intent(inout) :: lv
integer(psb_ipk_), intent(out) :: info integer(psb_ipk_), intent(out) :: info
! !
integer(psb_ipk_) :: nwv integer(psb_ipk_) :: nwv,i
info = psb_success_ info = psb_success_
call lv%wrk%vx2l%free(info) call lv%wrk%vx2l%free(info)
call lv%wrk%vy2l%free(info) call lv%wrk%vy2l%free(info)
call lv%wrk%vtx%free(info) call lv%wrk%vtx%free(info)
call lv%wrk%vty%free(info) call lv%wrk%vty%free(info)
nwv = size(lv%wrk%wv)
do i=1,nwv
call lv%wrk%wv(i)%free(info)
end do
end subroutine z_base_onelev_free_wrk end subroutine z_base_onelev_free_wrk
end module mld_z_onelev_mod end module mld_z_onelev_mod

Write Preview
Loading…
Cancel
Save