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

mld2p4-2:

Browse Source 
mlprec/impl/mld_caggrmat_biz_asb.f90
 mlprec/impl/mld_caggrmat_minnrg_asb.f90
 mlprec/impl/mld_cmlprec_aply.f90
 mlprec/impl/mld_daggrmat_biz_asb.f90
 mlprec/impl/mld_daggrmat_minnrg_asb.f90
 mlprec/impl/mld_dmlprec_aply.f90
 mlprec/impl/mld_saggrmat_biz_asb.f90
 mlprec/impl/mld_saggrmat_minnrg_asb.f90
 mlprec/impl/mld_smlprec_aply.f90
 mlprec/impl/mld_zaggrmat_biz_asb.f90
 mlprec/impl/mld_zaggrmat_minnrg_asb.f90
 mlprec/impl/mld_zmlprec_aply.f90

Fixed use of op_prol vs tmp_prol in coarse matrix build. 
Fixed implementation of KCycle/KCycleSym. Seems to be working now.
stopcriterion
Salvatore Filippone 9 years ago
parent b33a57b3bb
commit 3c6fc5d14d
12 changed files with 228 additions and 204 deletions
Show Stats Download Patch File Download Diff File

14
mlprec/impl/mld_caggrmat_biz_asb.f90
Unescape Escape View File

@ -97,7 +97,7 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: am3, am4,tmp_prol type(psb_cspmat_type) :: am3, am4
type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_coo_sparse_mat) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_spk_), allocatable :: adiag(:) complex(psb_spk_), allocatable :: adiag(:)
@ -350,29 +350,29 @@ subroutine mld_caggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call op_prol%mv_from(acsr1) call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info) call psb_rwextd(ncol,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_
call tmp_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info) call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')

11
mlprec/impl/mld_caggrmat_minnrg_asb.f90
Unescape Escape View File

@ -120,7 +120,7 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name character(len=20) :: name
type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_cspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_cspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_cspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol type(psb_cspmat_type) :: dat, datp, datdatp, atmp3
type(psb_c_coo_sparse_mat) :: tmpcoo type(psb_c_coo_sparse_mat) :: tmpcoo
type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_c_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_c_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
@ -512,10 +512,9 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Now we have to gather the halo of op_prol, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call op_prol%clone(tmp_prol,info) call psb_sphalo(op_prol,desc_a,am4,info,&
if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
@ -548,13 +547,13 @@ subroutine mld_caggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'

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

@ -474,7 +474,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_ml_aply at level ',level write(debug_unit,*) me,' Start inner_ml_aply at level ',level
end if end if
select case(p%precv(level)%parms%ml_type) select case(p%precv(level)%parms%ml_type)
@ -538,6 +538,9 @@ contains
goto 9999 goto 9999
end select end select
if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -603,7 +606,7 @@ contains
call p%precv(level)%sm%apply(cone,& call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info,init='Z')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply') & a_err='Error during ADD smoother_apply')
@ -615,7 +618,6 @@ contains
call psb_map_X2Y(cone,mlprec_wrk(level)%vx2l,& call psb_map_X2Y(cone,mlprec_wrk(level)%vx2l,&
& czero,mlprec_wrk(level+1)%vx2l,& & czero,mlprec_wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
call mlprec_wrk(level+1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -917,7 +919,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,name,' at level ',level write(debug_unit,*) me,name,' start at level ',level
end if end if
if ((level<1).or.(level>nlev)) then if ((level<1).or.(level>nlev)) then
@ -937,7 +939,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -946,13 +948,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(cone,& if (info == psb_success_) call p%precv(level)%sm%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(cone,& if (info == psb_success_) call p%precv(level)%sm2%apply(cone,&
& mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -983,7 +985,7 @@ contains
call psb_map_X2Y(cone,mlprec_wrk(level)%vty,& call psb_map_X2Y(cone,mlprec_wrk(level)%vty,&
& czero,mlprec_wrk(level + 1)%vx2l,& & czero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work)
call mlprec_wrk(level + 1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -992,11 +994,11 @@ contains
!Set the preconditioner !Set the preconditioner
if (level < nlev ) then if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') call mld_cinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
else else
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_type') & a_err='Bad value for ml_type')
@ -1091,7 +1093,8 @@ contains
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
complex(psb_spk_),target :: work(:) complex(psb_spk_),target :: work(:)
!Other variables !Other variables
@ -1099,9 +1102,10 @@ contains
type(psb_c_vect_type), dimension(0:1) :: d type(psb_c_vect_type), dimension(0:1) :: d
complex(psb_spk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta complex(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_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
complex(psb_spk_), allocatable :: temp_v(:) complex(psb_spk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
@ -1120,6 +1124,14 @@ contains
call psb_geasb(x,& call psb_geasb(x,&
& p%precv(level)%base_desc,info,& & p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call x%zero() call x%zero()
@ -1133,30 +1145,20 @@ contains
nc2l = p%precv(level)%base_desc%get_local_cols() nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_spk_)') & a_err='complex(psb_spk_)')
goto 9999 goto 9999
end if end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info) delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner !Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call mlprec_wrk(level)%vy2l%set(czero)
idx=0 idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) 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_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) call psb_spmm(cone,p%precv(level)%base_a,d(idx),czero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1165,24 +1167,27 @@ contains
end if end if
!FCG !FCG
if (innersolv == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) 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) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR !GCR
else else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info) tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
alpha = delta_old/tau alpha = delta_old/tau
!Update residual w !Update residual w
call psb_geaxpby(-alpha, v, cone, w, p%precv(level)%base_desc, info) 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) l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
iter = 0 iter = 0
if (l2_norm <= rtol*delta0) then
if (l2_norm <= rtol*delta) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info) call psb_geaxpby(alpha, d(idx), cone, x, p%precv(level)%base_desc, info)
else else
@ -1204,18 +1209,20 @@ contains
end if end if
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
!FCG if (psb_toupper(trim(innersolv)) == 'FCG') then
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, 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) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
!CGR else if (psb_toupper(trim(innersolv)) == 'GCR') then
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, 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) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
!Update solution !Update solution
@ -1225,8 +1232,8 @@ contains
call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info) call psb_geaxpby(alpha,d(idx),cone,x,p%precv(level)%base_desc,info)
endif endif
!Free vectors
call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, 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(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info)

14
mlprec/impl/mld_daggrmat_biz_asb.f90
Unescape Escape View File

@ -97,7 +97,7 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: am3, am4,tmp_prol type(psb_dspmat_type) :: am3, am4
type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_coo_sparse_mat) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_dpk_), allocatable :: adiag(:) real(psb_dpk_), allocatable :: adiag(:)
@ -350,29 +350,29 @@ subroutine mld_daggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call op_prol%mv_from(acsr1) call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info) call psb_rwextd(ncol,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_
call tmp_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info) call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')

11
mlprec/impl/mld_daggrmat_minnrg_asb.f90
Unescape Escape View File

@ -120,7 +120,7 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name character(len=20) :: name
type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_dspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_dspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_dspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol type(psb_dspmat_type) :: dat, datp, datdatp, atmp3
type(psb_d_coo_sparse_mat) :: tmpcoo type(psb_d_coo_sparse_mat) :: tmpcoo
type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_d_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_d_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
@ -512,10 +512,9 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Now we have to gather the halo of op_prol, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call op_prol%clone(tmp_prol,info) call psb_sphalo(op_prol,desc_a,am4,info,&
if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
@ -548,13 +547,13 @@ subroutine mld_daggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'

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

@ -474,7 +474,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_ml_aply at level ',level write(debug_unit,*) me,' Start inner_ml_aply at level ',level
end if end if
select case(p%precv(level)%parms%ml_type) select case(p%precv(level)%parms%ml_type)
@ -538,6 +538,9 @@ contains
goto 9999 goto 9999
end select end select
if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -603,7 +606,7 @@ contains
call p%precv(level)%sm%apply(done,& call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info,init='Z')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply') & a_err='Error during ADD smoother_apply')
@ -615,7 +618,6 @@ contains
call psb_map_X2Y(done,mlprec_wrk(level)%vx2l,& call psb_map_X2Y(done,mlprec_wrk(level)%vx2l,&
& dzero,mlprec_wrk(level+1)%vx2l,& & dzero,mlprec_wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
call mlprec_wrk(level+1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -917,7 +919,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,name,' at level ',level write(debug_unit,*) me,name,' start at level ',level
end if end if
if ((level<1).or.(level>nlev)) then if ((level<1).or.(level>nlev)) then
@ -937,7 +939,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -946,13 +948,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(done,& if (info == psb_success_) call p%precv(level)%sm%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(done,& if (info == psb_success_) call p%precv(level)%sm2%apply(done,&
& mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -983,7 +985,7 @@ contains
call psb_map_X2Y(done,mlprec_wrk(level)%vty,& call psb_map_X2Y(done,mlprec_wrk(level)%vty,&
& dzero,mlprec_wrk(level + 1)%vx2l,& & dzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work)
call mlprec_wrk(level + 1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -992,11 +994,11 @@ contains
!Set the preconditioner !Set the preconditioner
if (level < nlev ) then if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
else else
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_type') & a_err='Bad value for ml_type')
@ -1091,7 +1093,8 @@ contains
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
real(psb_dpk_),target :: work(:) real(psb_dpk_),target :: work(:)
!Other variables !Other variables
@ -1099,9 +1102,10 @@ contains
type(psb_d_vect_type), dimension(0:1) :: d 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_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, delta, rtol=0.25 real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
real(psb_dpk_), allocatable :: temp_v(:) real(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
@ -1120,6 +1124,14 @@ contains
call psb_geasb(x,& call psb_geasb(x,&
& p%precv(level)%base_desc,info,& & p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call x%zero() call x%zero()
@ -1133,30 +1145,20 @@ contains
nc2l = p%precv(level)%base_desc%get_local_cols() nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_dpk_)') & a_err='real(psb_dpk_)')
goto 9999 goto 9999
end if end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info) delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner !Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call mlprec_wrk(level)%vy2l%set(dzero)
idx=0 idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) 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_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) call psb_spmm(done,p%precv(level)%base_a,d(idx),dzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1165,24 +1167,27 @@ contains
end if end if
!FCG !FCG
if (innersolv == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) 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) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR !GCR
else else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info) tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
alpha = delta_old/tau alpha = delta_old/tau
!Update residual w !Update residual w
call psb_geaxpby(-alpha, v, done, w, p%precv(level)%base_desc, info) 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) l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
iter = 0 iter = 0
if (l2_norm <= rtol*delta0) then
if (l2_norm <= rtol*delta) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info) call psb_geaxpby(alpha, d(idx), done, x, p%precv(level)%base_desc, info)
else else
@ -1204,18 +1209,20 @@ contains
end if end if
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
!FCG if (psb_toupper(trim(innersolv)) == 'FCG') then
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, 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) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
!CGR else if (psb_toupper(trim(innersolv)) == 'GCR') then
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, 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) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
!Update solution !Update solution
@ -1225,8 +1232,8 @@ contains
call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info) call psb_geaxpby(alpha,d(idx),done,x,p%precv(level)%base_desc,info)
endif endif
!Free vectors
call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, 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(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info)

14
mlprec/impl/mld_saggrmat_biz_asb.f90
Unescape Escape View File

@ -97,7 +97,7 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: am3, am4,tmp_prol type(psb_sspmat_type) :: am3, am4
type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_coo_sparse_mat) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
real(psb_spk_), allocatable :: adiag(:) real(psb_spk_), allocatable :: adiag(:)
@ -350,29 +350,29 @@ subroutine mld_saggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call op_prol%mv_from(acsr1) call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info) call psb_rwextd(ncol,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_
call tmp_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info) call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')

11
mlprec/impl/mld_saggrmat_minnrg_asb.f90
Unescape Escape View File

@ -120,7 +120,7 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name character(len=20) :: name
type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_sspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_sspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_sspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol type(psb_sspmat_type) :: dat, datp, datdatp, atmp3
type(psb_s_coo_sparse_mat) :: tmpcoo type(psb_s_coo_sparse_mat) :: tmpcoo
type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_s_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_s_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
@ -512,10 +512,9 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Now we have to gather the halo of op_prol, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call op_prol%clone(tmp_prol,info) call psb_sphalo(op_prol,desc_a,am4,info,&
if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
@ -548,13 +547,13 @@ subroutine mld_saggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'

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

@ -474,7 +474,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_ml_aply at level ',level write(debug_unit,*) me,' Start inner_ml_aply at level ',level
end if end if
select case(p%precv(level)%parms%ml_type) select case(p%precv(level)%parms%ml_type)
@ -538,6 +538,9 @@ contains
goto 9999 goto 9999
end select end select
if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -603,7 +606,7 @@ contains
call p%precv(level)%sm%apply(sone,& call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info,init='Z')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply') & a_err='Error during ADD smoother_apply')
@ -615,7 +618,6 @@ contains
call psb_map_X2Y(sone,mlprec_wrk(level)%vx2l,& call psb_map_X2Y(sone,mlprec_wrk(level)%vx2l,&
& szero,mlprec_wrk(level+1)%vx2l,& & szero,mlprec_wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
call mlprec_wrk(level+1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -917,7 +919,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,name,' at level ',level write(debug_unit,*) me,name,' start at level ',level
end if end if
if ((level<1).or.(level>nlev)) then if ((level<1).or.(level>nlev)) then
@ -937,7 +939,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -946,13 +948,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(sone,& if (info == psb_success_) call p%precv(level)%sm%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(sone,& if (info == psb_success_) call p%precv(level)%sm2%apply(sone,&
& mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -983,7 +985,7 @@ contains
call psb_map_X2Y(sone,mlprec_wrk(level)%vty,& call psb_map_X2Y(sone,mlprec_wrk(level)%vty,&
& szero,mlprec_wrk(level + 1)%vx2l,& & szero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work)
call mlprec_wrk(level + 1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -992,11 +994,11 @@ contains
!Set the preconditioner !Set the preconditioner
if (level < nlev ) then if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') call mld_sinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
else else
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_type') & a_err='Bad value for ml_type')
@ -1091,7 +1093,8 @@ contains
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
real(psb_spk_),target :: work(:) real(psb_spk_),target :: work(:)
!Other variables !Other variables
@ -1099,9 +1102,10 @@ contains
type(psb_s_vect_type), dimension(0:1) :: d 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_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_spk_) :: l2_norm, delta, rtol=0.25 real(psb_spk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
real(psb_spk_), allocatable :: temp_v(:) real(psb_spk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
@ -1120,6 +1124,14 @@ contains
call psb_geasb(x,& call psb_geasb(x,&
& p%precv(level)%base_desc,info,& & p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call x%zero() call x%zero()
@ -1133,30 +1145,20 @@ contains
nc2l = p%precv(level)%base_desc%get_local_cols() nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_spk_)') & a_err='real(psb_spk_)')
goto 9999 goto 9999
end if end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info) delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner !Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call mlprec_wrk(level)%vy2l%set(szero)
idx=0 idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) 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_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) call psb_spmm(sone,p%precv(level)%base_a,d(idx),szero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1165,24 +1167,27 @@ contains
end if end if
!FCG !FCG
if (innersolv == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) 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) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR !GCR
else else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info) tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
alpha = delta_old/tau alpha = delta_old/tau
!Update residual w !Update residual w
call psb_geaxpby(-alpha, v, sone, w, p%precv(level)%base_desc, info) 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) l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
iter = 0 iter = 0
if (l2_norm <= rtol*delta0) then
if (l2_norm <= rtol*delta) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info) call psb_geaxpby(alpha, d(idx), sone, x, p%precv(level)%base_desc, info)
else else
@ -1204,18 +1209,20 @@ contains
end if end if
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
!FCG if (psb_toupper(trim(innersolv)) == 'FCG') then
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, 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) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
!CGR else if (psb_toupper(trim(innersolv)) == 'GCR') then
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, 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) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
!Update solution !Update solution
@ -1225,8 +1232,8 @@ contains
call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info) call psb_geaxpby(alpha,d(idx),sone,x,p%precv(level)%base_desc,info)
endif endif
!Free vectors
call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, 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(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info)

14
mlprec/impl/mld_zaggrmat_biz_asb.f90
Unescape Escape View File

@ -97,7 +97,7 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
& naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw, err_act
integer(psb_ipk_) ::ictxt, np, me integer(psb_ipk_) ::ictxt, np, me
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: am3, am4,tmp_prol type(psb_zspmat_type) :: am3, am4
type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_coo_sparse_mat) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsrf, ptilde
complex(psb_dpk_), allocatable :: adiag(:) complex(psb_dpk_), allocatable :: adiag(:)
@ -350,29 +350,29 @@ subroutine mld_zaggrmat_biz_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_restr
call acsr1%set_dupl(psb_dupl_add_) call acsr1%set_dupl(psb_dupl_add_)
call op_prol%mv_from(acsr1) call op_prol%mv_from(acsr1)
call op_prol%clone(tmp_prol,info)
call psb_rwextd(ncol,tmp_prol,info) call psb_rwextd(ncol,op_prol,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol') call psb_errpush(psb_err_internal_error_,name,a_err='Halo of op_prol')
goto 9999 goto 9999
end if end if
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2') call psb_errpush(psb_err_from_subroutine_,name,a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_ & 'Done NUMBMM 2',parms%aggr_kind, mld_smooth_prol_
call tmp_prol%transp(op_restr) call op_prol%transp(op_restr)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call tmp_prol%free()
call psb_rwextd(ncol,am3,info) call psb_rwextd(ncol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3') call psb_errpush(psb_err_internal_error_,name,a_err='Extend am3')

11
mlprec/impl/mld_zaggrmat_minnrg_asb.f90
Unescape Escape View File

@ -120,7 +120,7 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
character(len=20) :: name character(len=20) :: name
type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp type(psb_zspmat_type) :: af, ptilde, rtilde, atran, atp, atdatp
type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da type(psb_zspmat_type) :: am3,am4, ap, adap,atmp,rada, ra, atmp2, dap, dadap, da
type(psb_zspmat_type) :: dat, datp, datdatp, atmp3, tmp_prol type(psb_zspmat_type) :: dat, datp, datdatp, atmp3
type(psb_z_coo_sparse_mat) :: tmpcoo type(psb_z_coo_sparse_mat) :: tmpcoo
type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf type(psb_z_csr_sparse_mat) :: acsr1, acsr2, acsr3, acsr, acsrf
type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc type(psb_z_csc_sparse_mat) :: csc_dap, csc_dadap, csc_datp, csc_datdatp, acsc
@ -512,10 +512,9 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
! Now we have to gather the halo of op_prol, and add it to itself ! Now we have to gather the halo of op_prol, and add it to itself
! to multiply it by A, ! to multiply it by A,
! !
call op_prol%clone(tmp_prol,info) call psb_sphalo(op_prol,desc_a,am4,info,&
if (info == psb_success_) call psb_sphalo(tmp_prol,desc_a,am4,info,&
& colcnv=.false.,rowscale=.true.) & colcnv=.false.,rowscale=.true.)
if (info == psb_success_) call psb_rwextd(ncol,tmp_prol,info,b=am4) if (info == psb_success_) call psb_rwextd(ncol,op_prol,info,b=am4)
if (info == psb_success_) call am4%free() if (info == psb_success_) call am4%free()
if(info /= psb_success_) then if(info /= psb_success_) then
@ -548,13 +547,13 @@ subroutine mld_zaggrmat_minnrg_asb(a,desc_a,ilaggr,nlaggr,parms,ac,op_prol,op_re
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'starting sphalo/ rwxtd' & 'starting sphalo/ rwxtd'
call psb_symbmm(a,tmp_prol,am3,info) call psb_symbmm(a,op_prol,am3,info)
if(info /= psb_success_) then if(info /= psb_success_) then
call psb_errpush(psb_err_from_subroutine_,name,& call psb_errpush(psb_err_from_subroutine_,name,&
& a_err='symbmm 2') & a_err='symbmm 2')
goto 9999 goto 9999
end if end if
call psb_numbmm(a,tmp_prol,am3) call psb_numbmm(a,op_prol,am3)
if (debug_level >= psb_debug_outer_) & if (debug_level >= psb_debug_outer_) &
& write(debug_unit,*) me,' ',trim(name),& & write(debug_unit,*) me,' ',trim(name),&
& 'Done NUMBMM 2' & 'Done NUMBMM 2'

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

@ -474,7 +474,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,' inner_ml_aply at level ',level write(debug_unit,*) me,' Start inner_ml_aply at level ',level
end if end if
select case(p%precv(level)%parms%ml_type) select case(p%precv(level)%parms%ml_type)
@ -538,6 +538,9 @@ contains
goto 9999 goto 9999
end select end select
if(debug_level > 1) then
write(debug_unit,*) me,' End inner_ml_aply at level ',level
end if
call psb_erractionrestore(err_act) call psb_erractionrestore(err_act)
return return
@ -603,7 +606,7 @@ contains
call p%precv(level)%sm%apply(zone,& call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info) & sweeps,work,info,init='Z')
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during ADD smoother_apply') & a_err='Error during ADD smoother_apply')
@ -615,7 +618,6 @@ contains
call psb_map_X2Y(zone,mlprec_wrk(level)%vx2l,& call psb_map_X2Y(zone,mlprec_wrk(level)%vx2l,&
& zzero,mlprec_wrk(level+1)%vx2l,& & zzero,mlprec_wrk(level+1)%vx2l,&
& p%precv(level+1)%map,info,work=work) & p%precv(level+1)%map,info,work=work)
call mlprec_wrk(level+1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -917,7 +919,7 @@ contains
call psb_info(ictxt, me, np) call psb_info(ictxt, me, np)
if(debug_level > 1) then if(debug_level > 1) then
write(debug_unit,*) me,name,' at level ',level write(debug_unit,*) me,name,' start at level ',level
end if end if
if ((level<1).or.(level>nlev)) then if ((level<1).or.(level>nlev)) then
@ -937,7 +939,7 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else if (level < nlev) then else if (level < nlev) then
@ -946,13 +948,13 @@ contains
if (info == psb_success_) call p%precv(level)%sm%apply(zone,& if (info == psb_success_) call p%precv(level)%sm%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
else else
sweeps = p%precv(level)%parms%sweeps_post sweeps = p%precv(level)%parms%sweeps_post
if (info == psb_success_) call p%precv(level)%sm2%apply(zone,& if (info == psb_success_) call p%precv(level)%sm2%apply(zone,&
& mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,&
& p%precv(level)%base_desc, trans,& & p%precv(level)%base_desc, trans,&
& sweeps,work,info,init='Y') & sweeps,work,info,init='Z')
end if end if
if (info /= psb_success_) then if (info /= psb_success_) then
@ -983,7 +985,7 @@ contains
call psb_map_X2Y(zone,mlprec_wrk(level)%vty,& call psb_map_X2Y(zone,mlprec_wrk(level)%vty,&
& zzero,mlprec_wrk(level + 1)%vx2l,& & zzero,mlprec_wrk(level + 1)%vx2l,&
& p%precv(level + 1)%map,info,work=work) & p%precv(level + 1)%map,info,work=work)
call mlprec_wrk(level + 1)%vy2l%zero()
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Error during restriction') & a_err='Error during restriction')
@ -992,11 +994,11 @@ contains
!Set the preconditioner !Set the preconditioner
if (level < nlev ) then if (level <= nlev - 2 ) then
if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG')
elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then
call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') call mld_zinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'GCR')
else else
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
& a_err='Bad value for ml_type') & a_err='Bad value for ml_type')
@ -1091,7 +1093,8 @@ contains
type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:)
integer(psb_ipk_), intent(in) :: level integer(psb_ipk_), intent(in) :: level
character, intent(in) :: trans, innersolv character, intent(in) :: trans
character(len=*), intent(in) :: innersolv
complex(psb_dpk_),target :: work(:) complex(psb_dpk_),target :: work(:)
!Other variables !Other variables
@ -1099,9 +1102,10 @@ contains
type(psb_z_vect_type), dimension(0:1) :: d type(psb_z_vect_type), dimension(0:1) :: d
complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta complex(psb_dpk_) :: delta_old, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta
real(psb_dpk_) :: l2_norm, delta, rtol=0.25 real(psb_dpk_) :: l2_norm, delta, rtol=0.25, delta0, tnrm
complex(psb_dpk_), allocatable :: temp_v(:) complex(psb_dpk_), allocatable :: temp_v(:)
integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx integer(psb_ipk_) :: info, nlev, i, iter, max_iter=2, idx
character(len=20) :: name = 'innerit_k_cycle'
!Assemble rhs, w, v, v1, x !Assemble rhs, w, v, v1, x
@ -1120,6 +1124,14 @@ contains
call psb_geasb(x,& call psb_geasb(x,&
& p%precv(level)%base_desc,info,& & p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vx2l%v) & scratch=.true.,mold=mlprec_wrk(level)%vx2l%v)
!Assemble d(0) and d(1)
call psb_geasb(d(0),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call psb_geasb(d(1),&
& p%precv(level)%base_desc,info,&
& scratch=.true.,mold=mlprec_wrk(level)%vy2l%v)
call x%zero() call x%zero()
@ -1133,30 +1145,20 @@ contains
nc2l = p%precv(level)%base_desc%get_local_cols() nc2l = p%precv(level)%base_desc%get_local_cols()
info=psb_err_alloc_request_ info=psb_err_alloc_request_
call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),& call psb_errpush(info,name,i_err=(/2*nc2l,izero,izero,izero,izero/),&
& a_err='TYPE@(psb_dpk_)') & a_err='complex(psb_dpk_)')
goto 9999 goto 9999
end if end if
delta = psb_gedot(w, w, p%precv(level)%base_desc, info) delta0 = psb_genrm2(w, p%precv(level)%base_desc, info)
!Apply the preconditioner !Apply the preconditioner
call mlprec_wrk(level)%vy2l%zero()
call mlprec_wrk(level)%vy2l%set(zzero)
idx=0 idx=0
call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) 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_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) call psb_spmm(zone,p%precv(level)%base_a,d(idx),zzero,v,p%precv(level)%base_desc,info)
if (info /= psb_success_) then if (info /= psb_success_) then
call psb_errpush(psb_err_internal_error_,name,& call psb_errpush(psb_err_internal_error_,name,&
@ -1165,24 +1167,27 @@ contains
end if end if
!FCG !FCG
if (innersolv == 'FCG') then if (psb_toupper(trim(innersolv)) == 'FCG') then
delta_old = psb_gedot(d(idx), w, p%precv(level)%base_desc, info) 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) tau = psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
!CGR !GCR
else else if (psb_toupper(trim(innersolv)) == 'GCR') then
delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info) delta_old = psb_gedot(v, w, p%precv(level)%base_desc, info)
tau = psb_gedot(v, v, p%precv(level)%base_desc, info) tau = psb_gedot(v, v, p%precv(level)%base_desc, info)
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
alpha = delta_old/tau alpha = delta_old/tau
!Update residual w !Update residual w
call psb_geaxpby(-alpha, v, zone, w, p%precv(level)%base_desc, info) 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) l2_norm = psb_genrm2(w, p%precv(level)%base_desc, info)
iter = 0 iter = 0
if (l2_norm <= rtol*delta0) then
if (l2_norm <= rtol*delta) then
!Update solution x !Update solution x
call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info) call psb_geaxpby(alpha, d(idx), zone, x, p%precv(level)%base_desc, info)
else else
@ -1204,18 +1209,20 @@ contains
end if end if
!tau1, tau2, tau3, tau4 !tau1, tau2, tau3, tau4
!FCG if (psb_toupper(trim(innersolv)) == 'FCG') then
if (innersolv == 'FCG') then
tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info) tau1= psb_gedot(d(idx), v, p%precv(level)%base_desc, info)
tau2= psb_gedot(d(idx), v1, 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) tau3= psb_gedot(d(idx), w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
!CGR else if (psb_toupper(trim(innersolv)) == 'GCR') then
else
tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info) tau1= psb_gedot(v1, v, p%precv(level)%base_desc, info)
tau2= psb_gedot(v1, v1, 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) tau3= psb_gedot(v1, w, p%precv(level)%base_desc, info)
tau4= tau2 - (tau1*tau1)/tau tau4= tau2 - (tau1*tau1)/tau
else
call psb_errpush(psb_err_internal_error_,name,&
& a_err='Invalid inner solver')
goto 9999
endif endif
!Update solution !Update solution
@ -1225,8 +1232,8 @@ contains
call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info) call psb_geaxpby(alpha,d(idx),zone,x,p%precv(level)%base_desc,info)
endif endif
!Free vectors
call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info) call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vy2l,p%precv(level)%base_desc,info)
!Free vectors
call psb_gefree(v, 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(v1, p%precv(level)%base_desc, info)
call psb_gefree(w, p%precv(level)%base_desc, info) call psb_gefree(w, p%precv(level)%base_desc, info)

Write Preview
Loading…
Cancel
Save