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

mld2p4-2:

Browse Source 
mlprec/impl/mld_cmlprec_aply.f90
 mlprec/impl/mld_dmlprec_aply.f90
 mlprec/impl/mld_dprecbld.f90
 mlprec/impl/mld_dprecinit.F90
 mlprec/impl/mld_smlprec_aply.f90
 mlprec/impl/mld_zmlprec_aply.f90
 tests/fileread/cf_sample.f90
 tests/fileread/df_sample.f90
 tests/fileread/runs/dfs.inp
 tests/fileread/sf_sample.f90
 tests/fileread/zf_sample.f90
 tests/pdegen/ppde2d.f90
 tests/pdegen/ppde3d.f90
 tests/pdegen/runs/ppde.inp
 tests/pdegen/spde2d.f90
 tests/pdegen/spde3d.f90

New mlprec_aply.
New precinit interface & choice of levels.
This commit is contained in:
Salvatore Filippone
2016-05-24 09:12:25 +00:00
parent 90602d29fb
commit 285f77c8d6
16 changed files with 2037 additions and 578 deletions
Download Patch File Download Diff File Expand all files Collapse all files
mlprec/impl/mld_cmlprec_aply.f90
+496 -12
View File
@@ -1111,7 +1111,7 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
contains
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info, U)
implicit none
@@ -1122,6 +1122,10 @@ contains
character, intent(in) :: trans
complex(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
type(psb_c_vect_type),intent(inout), optional :: u
type(psb_c_vect_type) :: res
! Local variables
integer(psb_ipk_) :: ictxt,np,me
@@ -1263,7 +1267,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1287,7 +1291,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(cone,&
@@ -1299,7 +1303,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C')
@@ -1360,7 +1364,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
@@ -1438,7 +1442,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
@@ -1479,7 +1483,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
@@ -1491,7 +1495,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1582,7 +1586,7 @@ contains
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
@@ -1608,7 +1612,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
@@ -1616,13 +1620,13 @@ contains
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1667,6 +1671,321 @@ contains
end select
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level-1)%vty,&
& czero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
! call mlprec_wrk(level)%vy2l%set(u%get_vect())
call psb_geaxpby(cone,u,&
& czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(czero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
cone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vtx,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vtx,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
!K cycle
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call psb_geaxpby(cone,u,&
& czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(czero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
cone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-cone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
! Apply the restriction
call psb_map_X2Y(cone,mlprec_wrk(level)%vty,&
& czero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!Set the preconditioner
if ((level < nlev - 2)) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR')
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,&
& cone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-cone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& cone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vtx,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vtx,cone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
@@ -1683,5 +2002,170 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_cinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_c_inner_mod, mld_protect_name => mld_cmlprec_aply
implicit none
!Input/Oputput variables
type(mld_cprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv
complex(psb_spk_),target :: work(:)
!Other variables
type(psb_c_vect_type) :: v, w, rhs, v1, x
type(psb_c_vect_type), dimension(0:1) :: d
complex(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_spk_) :: l2_norm, delta, rtol=0.25
complex(psb_spk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call x%set(czero)
! rhs=vx2l and w=rhs
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,czero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,czero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_spk_)')
goto 9999
end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%set(czero)
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geaxpby(cone,mlprec_wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!FCG
if (innersolv == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR
else
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
endif
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, cone, w, p%precv(level)%base_desc, info)
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0
if (l2_norm <= rtol*delta) then
!Update solution x
call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(cone,w,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(cone,mlprec_wrk(level)%vy2l,czero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v1,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!tau1, tau2, tau3, tau4
!FCG
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
!CGR
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
endif
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),cone,x,p%precv(level)%base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info)
endif
!Free vectors
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_cinneritkcycle
end subroutine mld_cmlprec_aply_vect
mlprec/impl/mld_dmlprec_aply.f90
+406 -398
View File
@@ -581,7 +581,6 @@ contains
end if
else
! Here at coarse level
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,&
@@ -1112,7 +1111,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
contains
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info, u)
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info, U)
implicit none
@@ -1127,6 +1126,7 @@ contains
type(psb_d_vect_type) :: res
! Local variables
integer(psb_ipk_) :: ictxt,np,me
integer(psb_ipk_) :: i, nr2l,nc2l,err_act
@@ -1267,7 +1267,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1291,7 +1291,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(done,&
@@ -1303,7 +1303,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C')
@@ -1364,7 +1364,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
@@ -1442,7 +1442,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
@@ -1483,7 +1483,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
@@ -1495,7 +1495,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1586,7 +1586,7 @@ contains
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
@@ -1612,7 +1612,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
@@ -1620,13 +1620,13 @@ contains
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1672,311 +1672,319 @@ contains
end select
case(mld_vcycle_ml_, mld_wcycle_ml_)
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vty,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
!V/W cycle
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level-1)%vty,&
& dzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
! call mlprec_wrk(level)%vy2l%set(u%get_vect())
call psb_geaxpby(done,u,&
& dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(dzero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
done, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call mlprec_wrk(level)%vy2l%set(u%get_vect())
else
call mlprec_wrk(level)%vy2l%set(dzero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
done, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
endif
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
!K cycle
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call psb_geaxpby(done,u,&
& dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(dzero)
endif
res = mlprec_wrk(level)%vx2l
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
done, res, p%precv(level)%base_desc, info, work=work, trans=trans)
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l)
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!Set the preconditioner
if ((level < nlev - 2)) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR')
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!K cycle
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call mlprec_wrk(level)%vy2l%set(u%get_vect())
else
call mlprec_wrk(level)%vy2l%set(dzero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
done, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
! Apply the restriction
call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!Set the preconditioner
if ((level < nlev - 2)) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR')
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
endif
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,&
& done,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case default
info = psb_err_from_subroutine_ai_
@@ -1995,169 +2003,169 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_d_inner_mod, mld_protect_name => mld_dmlprec_aply
recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_d_inner_mod, mld_protect_name => mld_dmlprec_aply
implicit none
implicit none
!Input/Oputpu variables
type(mld_dprec_type), intent(inout) :: p
!Input/Oputput variables
type(mld_dprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv
real(psb_dpk_),target :: work(:)
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv
real(psb_dpk_),target :: work(:)
!Other variables
type(psb_d_vect_type) :: v, w, rhs, v1, x
type(psb_d_vect_type), dimension(0:1) :: d
real(psb_dpk_) :: delta_old, delta, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
!Other variables
type(psb_d_vect_type) :: v, w, rhs, v1, x
type(psb_d_vect_type), dimension(0:1) :: d
real(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, rtol=0.25
real(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
real(psb_dpk_) :: l2_norm, delta, rtol=0.25
real(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
!Assemble rhs, w, v, v1, x
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call x%set(dzero)
call x%set(dzero)
! rhs=vx2l and w=rhs
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,w,&
& p%precv(level)%base_desc,info)
! rhs=vx2l and w=rhs
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(done,mlprec_wrk(level)%vx2l,dzero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='real(psb_dpk_)')
& a_err='TYPE@(psb_dpk_)')
goto 9999
end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info)
delta = psb_gedot(w, w, p%precv(level)%base_desc, info)
!Apply the preconditioner
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%set(dzero)
call mlprec_wrk(level)%vy2l%set(dzero)
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!FCG
if (innersolv == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR
else
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
endif
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info)
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0
if (l2_norm <= rtol*delta) then
!Update solution x
call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(done,w,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v1,p%precv(level)%base_desc,info)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
& a_err='Error during residue')
goto 9999
end if
!tau1, tau2, tau3, tau4
!FCG
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
!CGR
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
endif
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),done,x,p%precv(level)%base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info)
endif
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info)
!Free vectors
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0
if (l2_norm <= rtol*delta) then
!Update solution x
call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(done,w,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(done,mlprec_wrk(level)%vy2l,dzero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v1,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!tau1, tau2, tau3, tau4
!FCG
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
!CGR
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
endif
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),done,x,p%precv(level)%base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info)
endif
!Free vectors
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end if
return
end subroutine mld_dinneritkcycle
end subroutine mld_dinneritkcycle
end subroutine mld_dmlprec_aply_vect
mlprec/impl/mld_dprecbld.f90
+2 -10
View File
@@ -173,23 +173,15 @@ subroutine mld_dprecbld(a,desc_a,p,info,amold,vmold,imold)
& a_err='One level preconditioner check.')
goto 9999
endif
call p%precv(1)%sm%build(a,desc_a,upd_,info,&
& amold=amold,vmold=vmold,imold=imold)
if (info == 0) then
if (allocated(p%precv(1)%sm2a)) then
call p%precv(1)%sm%build(a,desc_a,upd_,info,&
& amold=amold,vmold=vmold,imold=imold)
p%precv(1)%sm2 => p%precv(1)%sm2a
else
p%precv(1)%sm2 => p%precv(i)%sm
end if
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='One level preconditioner build.')
goto 9999
endif
!
! Number of levels > 1
!
mlprec/impl/mld_dprecinit.F90
+2 -2
View File
@@ -169,10 +169,10 @@ subroutine mld_dprecinit(p,ptype,info,nlev)
if (present(nlev)) then
nlev_ = max(1,nlev)
!p%n_prec_levs = nlev_
p%n_prec_levs = nlev_
else
nlev_ = 3
!p%n_prec_levs = -ione
p%n_prec_levs = -ione
end if
ilev_ = 1
allocate(p%precv(nlev_),stat=info)
mlprec/impl/mld_smlprec_aply.f90
+496 -12
View File
@@ -1111,7 +1111,7 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
contains
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info, U)
implicit none
@@ -1122,6 +1122,10 @@ contains
character, intent(in) :: trans
real(psb_spk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
type(psb_s_vect_type),intent(inout), optional :: u
type(psb_s_vect_type) :: res
! Local variables
integer(psb_ipk_) :: ictxt,np,me
@@ -1263,7 +1267,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1287,7 +1291,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(sone,&
@@ -1299,7 +1303,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C')
@@ -1360,7 +1364,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
@@ -1438,7 +1442,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
@@ -1479,7 +1483,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
@@ -1491,7 +1495,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1582,7 +1586,7 @@ contains
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
@@ -1608,7 +1612,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
@@ -1616,13 +1620,13 @@ contains
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1667,6 +1671,321 @@ contains
end select
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level-1)%vty,&
& szero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
! call mlprec_wrk(level)%vy2l%set(u%get_vect())
call psb_geaxpby(sone,u,&
& szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(szero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
sone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vtx,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vtx,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
!K cycle
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call psb_geaxpby(sone,u,&
& szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(szero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
sone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-sone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
! Apply the restriction
call psb_map_X2Y(sone,mlprec_wrk(level)%vty,&
& szero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!Set the preconditioner
if ((level < nlev - 2)) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR')
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,&
& sone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-sone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& sone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vtx,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vtx,sone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
@@ -1683,5 +2002,170 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_sinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_s_inner_mod, mld_protect_name => mld_smlprec_aply
implicit none
!Input/Oputput variables
type(mld_sprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv
real(psb_spk_),target :: work(:)
!Other variables
type(psb_s_vect_type) :: v, w, rhs, v1, x
type(psb_s_vect_type), dimension(0:1) :: d
real(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_spk_) :: l2_norm, delta, rtol=0.25
real(psb_spk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call x%set(szero)
! rhs=vx2l and w=rhs
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,szero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,szero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_spk_)')
goto 9999
end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%set(szero)
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geaxpby(sone,mlprec_wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!FCG
if (innersolv == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR
else
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
endif
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, sone, w, p%precv(level)%base_desc, info)
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0
if (l2_norm <= rtol*delta) then
!Update solution x
call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(sone,w,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(sone,mlprec_wrk(level)%vy2l,szero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v1,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!tau1, tau2, tau3, tau4
!FCG
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
!CGR
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
endif
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),sone,x,p%precv(level)%base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info)
endif
!Free vectors
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_sinneritkcycle
end subroutine mld_smlprec_aply_vect
mlprec/impl/mld_zmlprec_aply.f90
+496 -12
View File
@@ -1111,7 +1111,7 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info)
contains
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
recursive subroutine inner_ml_aply(level,p,mlprec_wrk,trans,work,info, U)
implicit none
@@ -1122,6 +1122,10 @@ contains
character, intent(in) :: trans
complex(psb_dpk_),target :: work(:)
integer(psb_ipk_), intent(out) :: info
type(psb_z_vect_type),intent(inout), optional :: u
type(psb_z_vect_type) :: res
! Local variables
integer(psb_ipk_) :: ictxt,np,me
@@ -1263,7 +1267,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1287,7 +1291,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
else
sweeps = p%precv(level)%parms%sweeps
call p%precv(level)%sm2%apply(zone,&
@@ -1299,7 +1303,7 @@ contains
& a_err='Error during POST smoother_apply')
goto 9999
end if
end if
case('T','C')
@@ -1360,7 +1364,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
@@ -1438,7 +1442,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
@@ -1479,7 +1483,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
@@ -1491,7 +1495,7 @@ contains
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1582,7 +1586,7 @@ contains
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
@@ -1608,7 +1612,7 @@ contains
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
@@ -1616,13 +1620,13 @@ contains
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
@@ -1667,6 +1671,321 @@ contains
end select
case(mld_vcycle_ml_, mld_wcycle_ml_)
!V/W cycle
if (level > 1) then
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level-1)%vty,&
& zzero,mlprec_wrk(level)%vx2l,&
& p%precv(level)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
end if
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
! call mlprec_wrk(level)%vy2l%set(u%get_vect())
call psb_geaxpby(zone,u,&
& zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(zzero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
zone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info)
if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then
call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vtx,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vtx,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case(mld_kcycle_ml_, mld_kcyclesym_ml_)
!K cycle
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vtx,&
& p%precv(level)%base_desc,info)
!
! Apply the base preconditioner
!
if (level < nlev) then
if (present(u)) then
call psb_geaxpby(zone,u,&
& zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc,info)
else
call mlprec_wrk(level)%vy2l%set(zzero)
endif
res = mlprec_wrk(level)%vx2l
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
zone, res, p%precv(level)%base_desc, info, work=work, trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
else
sweeps = p%precv(level)%parms%sweeps
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-PRE smoother_apply')
goto 9999
end if
!
! Compute the residual (at all levels but the coarsest one)
! and call recursively
!
if(level < nlev) then
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info)
if (info == psb_success_) call psb_spmm(-zone,p%precv(level)%base_a,&
& mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vty,&
& p%precv(level)%base_desc,info,work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
! Apply the restriction
call psb_map_X2Y(zone,mlprec_wrk(level)%vty,&
& zzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction')
goto 9999
end if
!Set the preconditioner
if ((level < nlev - 2)) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR')
endif
else
call inner_ml_aply(level + 1 ,p,mlprec_wrk,trans,work,info)
endif
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error in recursive call')
goto 9999
end if
!
! Apply the prolongator
!
call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,&
& zone,mlprec_wrk(level)%vy2l,&
& p%precv(level+1)%map,info,work=work)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during prolongation')
goto 9999
end if
!
! Compute the residual
!
call psb_spmm(-zone,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,&
& zone,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,&
& work=work,trans=trans)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!
! Apply the base preconditioner
!
if (trans == 'N') then
sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vtx,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
else
sweeps = p%precv(level)%parms%sweeps_pre
if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vtx,zone,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,&
& sweeps,work,info)
end if
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during 2-POST smoother_apply')
goto 9999
end if
endif
case default
info = psb_err_from_subroutine_ai_
call psb_errpush(info,name,a_err='invalid mltype',&
@@ -1683,5 +2002,170 @@ contains
end subroutine inner_ml_aply
recursive subroutine mld_zinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv)
use psb_base_mod
use mld_prec_mod
use mld_z_inner_mod, mld_protect_name => mld_zmlprec_aply
implicit none
!Input/Oputput variables
type(mld_zprec_type), intent(inout) :: p
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv
complex(psb_dpk_),target :: work(:)
!Other variables
type(psb_z_vect_type) :: v, w, rhs, v1, x
type(psb_z_vect_type), dimension(0:1) :: d
complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, delta, rtol=0.25
complex(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
!Assemble rhs, w, v, v1, x
call psb_geasb(rhs,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(w,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(v1,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call psb_geasb(x,&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
call x%set(zzero)
! rhs=vx2l and w=rhs
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,zzero,rhs,&
& p%precv(level)%base_desc,info)
call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,zzero,w,&
& p%precv(level)%base_desc,info)
if (psb_errstatus_fatal()) then
nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_dpk_)')
goto 9999
end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info)
!Apply the preconditioner
call mlprec_wrk(level)%vy2l%set(zzero)
idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geaxpby(zone,mlprec_wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!FCG
if (innersolv == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR
else
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
endif
alpha = delta_old/tau
!Update residual w
call psb_geaxpby(-alpha, v, zone, w, p%precv(level)%base_desc, info)
l2_norm = psb_gedot(w, w, p%precv(level)%base_desc, info)
iter = 0
if (l2_norm <= rtol*delta) then
!Update solution x
call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info)
else
iter = iter + 1
idx=mod(iter,2)
!Apply preconditioner
call psb_geaxpby(zone,w,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info)
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info)
call psb_geaxpby(zone,mlprec_wrk(level)%vy2l,zzero,d(idx),p%precv(level)%base_desc,info)
!Sparse matrix vector product
call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v1,p%precv(level)%base_desc,info)
if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during residue')
goto 9999
end if
!tau1, tau2, tau3, tau4
!FCG
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
!CGR
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, p%precv(level)%base_desc, info)
tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau
endif
!Update solution
alpha=alpha-(tau1*tau3)/(tau*tau4)
call psb_geaxpby(alpha,d(idx - 1),zone,x,p%precv(level)%base_desc,info)
alpha=tau3/tau4
call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info)
endif
!Free vectors
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
call psb_gefree(v, p%precv(level)%base_desc, info)
call psb_gefree(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info)
call psb_gefree(x, p%precv(level)%base_desc, info)
call psb_gefree(d(0), p%precv(level)%base_desc, info)
call psb_gefree(d(1), p%precv(level)%base_desc, info)
9999 continue
call psb_erractionrestore(err_act)
if (err_act.eq.psb_act_abort_) then
call psb_error()
return
end if
return
end subroutine mld_zinneritkcycle
end subroutine mld_zmlprec_aply_vect
tests/fileread/cf_sample.f90
+31 -30
View File
@@ -254,40 +254,41 @@ program cf_sample
!
if (psb_toupper(prec_choice%prec) == 'ML') then
nlv = prec_choice%nlev
call mld_precinit(prec,prec_choice%prec,info,nlev=nlv)
call mld_precset(prec,mld_smoother_type_, prec_choice%smther, info)
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,mld_aggr_kind_, prec_choice%aggrkind,info)
call mld_precset(prec,mld_aggr_alg_, prec_choice%aggr_alg,info)
call mld_precset(prec,mld_aggr_ord_, prec_choice%aggr_ord,info)
call mld_precset(prec,mld_ml_type_, prec_choice%mltype, info)
call mld_precset(prec,mld_smoother_pos_, prec_choice%smthpos, info)
call mld_precset(prec,mld_aggr_scale_, prec_choice%ascale, info)
call mld_precset(prec,mld_aggr_thresh_, prec_choice%athres, info)
call mld_precset(prec,mld_coarse_solve_, prec_choice%csolve, info)
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
call mld_precinit(prec,prec_choice%prec, info)
if (prec_choice%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prec_choice%nlev, info)
call mld_precset(prec,'smoother_type', prec_choice%smther, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
call mld_precset(prec,'aggr_kind', prec_choice%aggrkind,info)
call mld_precset(prec,'aggr_alg', prec_choice%aggr_alg,info)
call mld_precset(prec,'aggr_ord', prec_choice%aggr_ord,info)
call mld_precset(prec,'ml_type', prec_choice%mltype, info)
call mld_precset(prec,'smoother_pos', prec_choice%smthpos, info)
call mld_precset(prec,'aggr_scale', prec_choice%ascale, info)
call mld_precset(prec,'aggr_thresh', prec_choice%athres, info)
call mld_precset(prec,'coarse_solve', prec_choice%csolve, info)
call mld_precset(prec,'coarse_subsolve', prec_choice%csbsolve,info)
call mld_precset(prec,'coarse_mat', prec_choice%cmat, info)
call mld_precset(prec,'coarse_fillin', prec_choice%cfill, info)
call mld_precset(prec,'coarse_iluthrs', prec_choice%cthres, info)
call mld_precset(prec,'coarse_sweeps', prec_choice%cjswp, info)
else
nlv = 1
call mld_precinit(prec,prec_choice%prec,info)
if (psb_toupper(prec_choice%prec) /= 'NONE') then
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
end if
end if
! building the preconditioner
tests/fileread/df_sample.f90
+31 -30
View File
@@ -254,40 +254,41 @@ program df_sample
!
if (psb_toupper(prec_choice%prec) == 'ML') then
nlv = prec_choice%nlev
call mld_precinit(prec,prec_choice%prec,info,nlev=nlv)
call mld_precset(prec,mld_smoother_type_, prec_choice%smther, info)
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,mld_aggr_kind_, prec_choice%aggrkind,info)
call mld_precset(prec,mld_aggr_alg_, prec_choice%aggr_alg,info)
call mld_precset(prec,mld_aggr_ord_, prec_choice%aggr_ord,info)
call mld_precset(prec,mld_ml_type_, prec_choice%mltype, info)
call mld_precset(prec,mld_smoother_pos_, prec_choice%smthpos, info)
call mld_precset(prec,mld_aggr_scale_, prec_choice%ascale, info)
call mld_precset(prec,mld_aggr_thresh_, prec_choice%athres, info)
call mld_precset(prec,mld_coarse_solve_, prec_choice%csolve, info)
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
call mld_precinit(prec,prec_choice%prec, info)
if (prec_choice%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prec_choice%nlev, info)
call mld_precset(prec,'smoother_type', prec_choice%smther, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
call mld_precset(prec,'aggr_kind', prec_choice%aggrkind,info)
call mld_precset(prec,'aggr_alg', prec_choice%aggr_alg,info)
call mld_precset(prec,'aggr_ord', prec_choice%aggr_ord,info)
call mld_precset(prec,'ml_type', prec_choice%mltype, info)
call mld_precset(prec,'smoother_pos', prec_choice%smthpos, info)
call mld_precset(prec,'aggr_scale', prec_choice%ascale, info)
call mld_precset(prec,'aggr_thresh', prec_choice%athres, info)
call mld_precset(prec,'coarse_solve', prec_choice%csolve, info)
call mld_precset(prec,'coarse_subsolve', prec_choice%csbsolve,info)
call mld_precset(prec,'coarse_mat', prec_choice%cmat, info)
call mld_precset(prec,'coarse_fillin', prec_choice%cfill, info)
call mld_precset(prec,'coarse_iluthrs', prec_choice%cthres, info)
call mld_precset(prec,'coarse_sweeps', prec_choice%cjswp, info)
else
nlv = 1
call mld_precinit(prec,prec_choice%prec,info)
if (psb_toupper(prec_choice%prec) /= 'NONE') then
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
end if
end if
! building the preconditioner
tests/fileread/runs/dfs.inp
+1 -1
View File
@@ -18,7 +18,7 @@ ILU ! AS subdomain solver: DSCALE ILU MILU ILUT UMF SLU
0 ! Fill level P for ILU(P) and ILU(T,P)
1.d-4 ! Threshold T for ILU(T,P)
2 ! Number of Jacobi sweeps for base smoother
4 ! Number of levels in a multilevel preconditioner
-1 ! Number of levels in a multilevel preconditioner; if <0, lib default
BJAC ! Smoother type JACOBI BJAC AS ignored for non-ML
SMOOTHED ! Type of aggregation: SMOOTHED NONSMOOTHED
DEC ! Type of aggregation: DEC
tests/fileread/sf_sample.f90
+31 -30
View File
@@ -254,40 +254,41 @@ program sf_sample
!
if (psb_toupper(prec_choice%prec) == 'ML') then
nlv = prec_choice%nlev
call mld_precinit(prec,prec_choice%prec,info,nlev=nlv)
call mld_precset(prec,mld_smoother_type_, prec_choice%smther, info)
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,mld_aggr_kind_, prec_choice%aggrkind,info)
call mld_precset(prec,mld_aggr_alg_, prec_choice%aggr_alg,info)
call mld_precset(prec,mld_aggr_ord_, prec_choice%aggr_ord,info)
call mld_precset(prec,mld_ml_type_, prec_choice%mltype, info)
call mld_precset(prec,mld_smoother_pos_, prec_choice%smthpos, info)
call mld_precset(prec,mld_aggr_scale_, prec_choice%ascale, info)
call mld_precset(prec,mld_aggr_thresh_, prec_choice%athres, info)
call mld_precset(prec,mld_coarse_solve_, prec_choice%csolve, info)
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
call mld_precinit(prec,prec_choice%prec, info)
if (prec_choice%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prec_choice%nlev, info)
call mld_precset(prec,'smoother_type', prec_choice%smther, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
call mld_precset(prec,'aggr_kind', prec_choice%aggrkind,info)
call mld_precset(prec,'aggr_alg', prec_choice%aggr_alg,info)
call mld_precset(prec,'aggr_ord', prec_choice%aggr_ord,info)
call mld_precset(prec,'ml_type', prec_choice%mltype, info)
call mld_precset(prec,'smoother_pos', prec_choice%smthpos, info)
call mld_precset(prec,'aggr_scale', prec_choice%ascale, info)
call mld_precset(prec,'aggr_thresh', prec_choice%athres, info)
call mld_precset(prec,'coarse_solve', prec_choice%csolve, info)
call mld_precset(prec,'coarse_subsolve', prec_choice%csbsolve,info)
call mld_precset(prec,'coarse_mat', prec_choice%cmat, info)
call mld_precset(prec,'coarse_fillin', prec_choice%cfill, info)
call mld_precset(prec,'coarse_iluthrs', prec_choice%cthres, info)
call mld_precset(prec,'coarse_sweeps', prec_choice%cjswp, info)
else
nlv = 1
call mld_precinit(prec,prec_choice%prec,info)
if (psb_toupper(prec_choice%prec) /= 'NONE') then
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
end if
end if
! building the preconditioner
tests/fileread/zf_sample.f90
+31 -30
View File
@@ -254,40 +254,41 @@ program zf_sample
!
if (psb_toupper(prec_choice%prec) == 'ML') then
nlv = prec_choice%nlev
call mld_precinit(prec,prec_choice%prec,info,nlev=nlv)
call mld_precset(prec,mld_smoother_type_, prec_choice%smther, info)
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,mld_aggr_kind_, prec_choice%aggrkind,info)
call mld_precset(prec,mld_aggr_alg_, prec_choice%aggr_alg,info)
call mld_precset(prec,mld_aggr_ord_, prec_choice%aggr_ord,info)
call mld_precset(prec,mld_ml_type_, prec_choice%mltype, info)
call mld_precset(prec,mld_smoother_pos_, prec_choice%smthpos, info)
call mld_precset(prec,mld_aggr_scale_, prec_choice%ascale, info)
call mld_precset(prec,mld_aggr_thresh_, prec_choice%athres, info)
call mld_precset(prec,mld_coarse_solve_, prec_choice%csolve, info)
call mld_precset(prec,mld_coarse_subsolve_, prec_choice%csbsolve,info)
call mld_precset(prec,mld_coarse_mat_, prec_choice%cmat, info)
call mld_precset(prec,mld_coarse_fillin_, prec_choice%cfill, info)
call mld_precset(prec,mld_coarse_iluthrs_, prec_choice%cthres, info)
call mld_precset(prec,mld_coarse_sweeps_, prec_choice%cjswp, info)
call mld_precinit(prec,prec_choice%prec, info)
if (prec_choice%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prec_choice%nlev, info)
call mld_precset(prec,'smoother_type', prec_choice%smther, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
call mld_precset(prec,'aggr_kind', prec_choice%aggrkind,info)
call mld_precset(prec,'aggr_alg', prec_choice%aggr_alg,info)
call mld_precset(prec,'aggr_ord', prec_choice%aggr_ord,info)
call mld_precset(prec,'ml_type', prec_choice%mltype, info)
call mld_precset(prec,'smoother_pos', prec_choice%smthpos, info)
call mld_precset(prec,'aggr_scale', prec_choice%ascale, info)
call mld_precset(prec,'aggr_thresh', prec_choice%athres, info)
call mld_precset(prec,'coarse_solve', prec_choice%csolve, info)
call mld_precset(prec,'coarse_subsolve', prec_choice%csbsolve,info)
call mld_precset(prec,'coarse_mat', prec_choice%cmat, info)
call mld_precset(prec,'coarse_fillin', prec_choice%cfill, info)
call mld_precset(prec,'coarse_iluthrs', prec_choice%cthres, info)
call mld_precset(prec,'coarse_sweeps', prec_choice%cjswp, info)
else
nlv = 1
call mld_precinit(prec,prec_choice%prec,info)
if (psb_toupper(prec_choice%prec) /= 'NONE') then
call mld_precset(prec,mld_smoother_sweeps_, prec_choice%jsweeps, info)
call mld_precset(prec,mld_sub_ovr_, prec_choice%novr, info)
call mld_precset(prec,mld_sub_restr_, prec_choice%restr, info)
call mld_precset(prec,mld_sub_prol_, prec_choice%prol, info)
call mld_precset(prec,mld_sub_solve_, prec_choice%solve, info)
call mld_precset(prec,mld_sub_fillin_, prec_choice%fill, info)
call mld_precset(prec,mld_sub_iluthrs_, prec_choice%thr, info)
call mld_precset(prec,'smoother_sweeps', prec_choice%jsweeps, info)
call mld_precset(prec,'sub_ovr', prec_choice%novr, info)
call mld_precset(prec,'sub_restr', prec_choice%restr, info)
call mld_precset(prec,'sub_prol', prec_choice%prol, info)
call mld_precset(prec,'sub_solve', prec_choice%solve, info)
call mld_precset(prec,'sub_fillin', prec_choice%fill, info)
call mld_precset(prec,'sub_iluthrs', prec_choice%thr, info)
end if
end if
! building the preconditioner
tests/pdegen/ppde2d.f90
+3 -2
View File
@@ -229,8 +229,9 @@ program ppde2d
!
if (psb_toupper(prectype%prec) == 'ML') then
nlv = prectype%nlev
call mld_precinit(prec,prectype%prec, info, nlev=nlv)
call mld_precinit(prec,prectype%prec, info)
if (prectype%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prectype%nlev, info)
call mld_precset(prec,'smoother_type', prectype%smther, info)
call mld_precset(prec,'smoother_sweeps', prectype%jsweeps, info)
call mld_precset(prec,'sub_ovr', prectype%novr, info)
tests/pdegen/ppde3d.f90
+3 -3
View File
@@ -241,9 +241,9 @@ program ppde3d
! prepare the preconditioner.
!
if (psb_toupper(prectype%prec) == 'ML') then
nlv = prectype%nlev
call mld_precinit(prec,prectype%prec, info, nlev=nlv)
!!$ call mld_precset(prec,'n_prec_levs', prectype%nlev, info)
call mld_precinit(prec,prectype%prec, info)
if (prectype%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prectype%nlev, info)
call mld_precset(prec,'smoother_type', prectype%smther, info)
call mld_precset(prec,'smoother_sweeps', prectype%jsweeps, info)
call mld_precset(prec,'sub_ovr', prectype%novr, info)
tests/pdegen/runs/ppde.inp
+2 -2
View File
@@ -17,11 +17,11 @@ GS ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU GS M
1.d-4 ! Threshold T for ILU(T,P)
4 ! Smoother/Jacobi sweeps
BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML
3 ! Number of levels in a multilevel preconditioner
-1 ! Number of levels in a multilevel preconditioner; if <0, lib default.
SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY
SYMDEC ! Type of aggregation DEC SYMDEC
NATURAL ! Ordering of aggregation NATURAL DEGREE
KCYCLE ! Type of multilevel correction: ADD MULT
KCYCLE ! Type of multilevel correction: ADD MULT KCYCLE VCYCLE WCYCLE KCYCLESYM
TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD)
DIST ! Coarse level: matrix distribution DIST REPL
BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST MUMPS
tests/pdegen/spde2d.f90
+3 -2
View File
@@ -229,8 +229,9 @@ program spde2d
!
if (psb_toupper(prectype%prec) == 'ML') then
nlv = prectype%nlev
call mld_precinit(prec,prectype%prec, info, nlev=nlv)
call mld_precinit(prec,prectype%prec, info)
if (prectype%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prectype%nlev, info)
call mld_precset(prec,'smoother_type', prectype%smther, info)
call mld_precset(prec,'smoother_sweeps', prectype%jsweeps, info)
call mld_precset(prec,'sub_ovr', prectype%novr, info)
tests/pdegen/spde3d.f90
+3 -2
View File
@@ -242,8 +242,9 @@ program spde3d
!
if (psb_toupper(prectype%prec) == 'ML') then
nlv = prectype%nlev
call mld_precinit(prec,prectype%prec, info, nlev=nlv)
call mld_precinit(prec,prectype%prec, info)
if (prectype%nlev > 0) &
& call mld_precset(prec,'n_prec_levs', prectype%nlev, info)
call mld_precset(prec,'smoother_type', prectype%smther, info)
call mld_precset(prec,'smoother_sweeps', prectype%jsweeps, info)
call mld_precset(prec,'sub_ovr', prectype%novr, info)