diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 4354f189..b9131014 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -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 diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index 11fe8643..b0d45972 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -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 + 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()) - else - call mlprec_wrk(level)%vy2l%set(dzero) - endif - res = mlprec_wrk(level)%vx2l + 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) - 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) + else + call mlprec_wrk(level)%vy2l%set(dzero) + endif + res = mlprec_wrk(level)%vx2l - 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) - 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_) 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)%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 2-PRE smoother_apply') + & a_err='Error during residue') goto 9999 end if - - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,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_) 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 + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + 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 + ! + ! 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 in recursive call') - goto 9999 - end if - + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + 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 + ! + ! 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 + 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 + endif - case(mld_kcycle_ml_, mld_kcyclesym_ml_) + case(mld_kcycle_ml_, mld_kcyclesym_ml_) - !K cycle + !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 + 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()) - else - call mlprec_wrk(level)%vy2l%set(dzero) - endif - res = mlprec_wrk(level)%vx2l + 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 - 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 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 (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 (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 2-PRE smoother_apply') + & a_err='Error during residue') goto 9999 end if - - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then + ! 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 - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) + !Set the preconditioner - 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 + 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 - !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 + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if - 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) - ! - ! 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 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 - 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(:) - - !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 - - 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 - - !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(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) + 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, rhs_norm, alpha, tau, tau1, tau2, tau3, tau4, beta + + 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 + + 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) + + ! 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) - - !Apply the preconditioner - - call mlprec_wrk(level)%vy2l%set(dzero) + delta = psb_gedot(w, w, p%precv(level)%base_desc, info) - idx=0 - call inner_ml_aply(level,p,mlprec_wrk,trans,work,info) + !Apply the preconditioner - !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 mlprec_wrk(level)%vy2l%set(dzero) + idx=0 + call inner_ml_aply(level,p,mlprec_wrk,trans,work,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 + !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) - !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) - 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) + + 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 - - !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) + 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) + 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 diff --git a/mlprec/impl/mld_dprecbld.f90 b/mlprec/impl/mld_dprecbld.f90 index 28259f3d..65238515 100644 --- a/mlprec/impl/mld_dprecbld.f90 +++ b/mlprec/impl/mld_dprecbld.f90 @@ -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 ! diff --git a/mlprec/impl/mld_dprecinit.F90 b/mlprec/impl/mld_dprecinit.F90 index 883957bb..7590f6ce 100644 --- a/mlprec/impl/mld_dprecinit.F90 +++ b/mlprec/impl/mld_dprecinit.F90 @@ -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) diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index f80db873..424e5861 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -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 diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index 6a4a3d02..56c83d16 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -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 diff --git a/tests/fileread/cf_sample.f90 b/tests/fileread/cf_sample.f90 index 8b409d7c..1643847f 100644 --- a/tests/fileread/cf_sample.f90 +++ b/tests/fileread/cf_sample.f90 @@ -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 diff --git a/tests/fileread/df_sample.f90 b/tests/fileread/df_sample.f90 index 713a64c0..4ab3b851 100644 --- a/tests/fileread/df_sample.f90 +++ b/tests/fileread/df_sample.f90 @@ -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 diff --git a/tests/fileread/runs/dfs.inp b/tests/fileread/runs/dfs.inp index 9aab8945..f3d03a9d 100644 --- a/tests/fileread/runs/dfs.inp +++ b/tests/fileread/runs/dfs.inp @@ -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 diff --git a/tests/fileread/sf_sample.f90 b/tests/fileread/sf_sample.f90 index e3a0e946..d9c97ae0 100644 --- a/tests/fileread/sf_sample.f90 +++ b/tests/fileread/sf_sample.f90 @@ -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 diff --git a/tests/fileread/zf_sample.f90 b/tests/fileread/zf_sample.f90 index 981e2d53..7b3112a4 100644 --- a/tests/fileread/zf_sample.f90 +++ b/tests/fileread/zf_sample.f90 @@ -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 diff --git a/tests/pdegen/ppde2d.f90 b/tests/pdegen/ppde2d.f90 index c96d4eae..a56b3018 100644 --- a/tests/pdegen/ppde2d.f90 +++ b/tests/pdegen/ppde2d.f90 @@ -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) diff --git a/tests/pdegen/ppde3d.f90 b/tests/pdegen/ppde3d.f90 index 0a302701..5f2c2b63 100644 --- a/tests/pdegen/ppde3d.f90 +++ b/tests/pdegen/ppde3d.f90 @@ -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) diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index d6056111..71b3c913 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -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 diff --git a/tests/pdegen/spde2d.f90 b/tests/pdegen/spde2d.f90 index ccff4dd4..3513aa2e 100644 --- a/tests/pdegen/spde2d.f90 +++ b/tests/pdegen/spde2d.f90 @@ -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) diff --git a/tests/pdegen/spde3d.f90 b/tests/pdegen/spde3d.f90 index 308cd03b..453095cc 100644 --- a/tests/pdegen/spde3d.f90 +++ b/tests/pdegen/spde3d.f90 @@ -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)