diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index b9131014..3941b1bf 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -39,7 +39,7 @@ ! File: mld_cmlprec_aply.f90 ! ! Subroutine: mld_cmlprec_aply -! Version: complex +! Version: real ! ! This routine computes ! @@ -84,7 +84,7 @@ ! The multilevel preconditioner data structure containing the ! local part of the preconditioner to be applied. ! Note that nlev = size(p%precv) = number of levels. -! p%precv(ilev)%prec - type(psb_dbaseprec_type) +! p%precv(ilev)%prec - type(psb_cbaseprec_type) ! The 'base preconditioner' for the current level ! p%precv(ilev)%ac - type(psb_cspmat_type) ! The local part of the matrix A(ilev). @@ -97,7 +97,7 @@ ! and prolongation operators described in the sequel. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! The integer parameters defining the multilevel -! strategy +! strategy ! p%precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable. ! The real parameters defining the multilevel strategy ! p%precv(ilev)%mlia - integer, dimension(:), allocatable. @@ -333,7 +333,7 @@ subroutine mld_cmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) type mld_mlprec_wrk_type complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_cmlprec_aply' info = psb_success_ @@ -1016,11 +1016,12 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act character(len=20) :: name character :: trans_ + complex(psb_spk_) :: par type mld_mlprec_wrk_type complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_cmlprec_aply' info = psb_success_ @@ -1036,7 +1037,6 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) & ' Entry ', size(p%precv) trans_ = psb_toupper(trans) - nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then @@ -1065,12 +1065,12 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do - level = 1 call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call mlprec_wrk(level)%vy2l%zero() + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -1082,6 +1082,7 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& & p%precv(level)%base_desc,info) do level = 1, nlev + call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vtx%free(info) @@ -1118,29 +1119,30 @@ contains ! Arguments integer(psb_ipk_) :: level type(mld_cprec_type), target, intent(inout) :: p - type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:) character, intent(in) :: trans complex(psb_spk_),target :: work(:) integer(psb_ipk_), intent(out) :: info type(psb_c_vect_type),intent(inout), optional :: u type(psb_c_vect_type) :: res - - + type(psb_c_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre ! Local variables integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: nlev, ilev, sweeps - + logical :: pre, post character(len=20) :: name + + name = 'inner_ml_aply' info = psb_success_ call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& @@ -1167,57 +1169,9 @@ contains goto 9999 case(mld_add_ml_) - ! - ! Additive multilevel - ! - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& - & 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 - - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error 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 + call mld_c_inner_add(p, mlprec_wrk, level, trans, work) - end if case(mld_mult_ml_) ! @@ -1230,492 +1184,250 @@ contains select case(p%precv(level)%parms%smoother_pos) case(mld_post_smooth_) + p%precv(level)%parms%sweeps_pre = 0 + call mld_c_inner_mult(p, mlprec_wrk, level, trans, work) + - select case (trans_) - case('N') - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& - & czero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + case(mld_pre_smooth_) + p%precv(level)%parms%sweeps_post = 0 + call mld_c_inner_mult(p, mlprec_wrk, level, trans, 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 + case(mld_twoside_smooth_) + call mld_c_inner_mult(p, mlprec_wrk, level, trans, work) - ! This is one step of post-smoothing + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid smooth_pos',& + & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) + goto 9999 - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + end select - ! - ! Apply the prolongator - ! - call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& - & czero,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)%vx2l,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 + case(mld_mult_dev_ml_) + call mld_c_inner_mult(p, mlprec_wrk, level, trans, work) - sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm2%apply(cone,& - & mlprec_wrk(level)%vx2l,cone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & mlprec_wrk(level)%vx2l,czero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + case(mld_vcycle_ml_, mld_wcycle_ml_) - end if + call mld_c_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u) - case('T','C') + case(mld_kcycle_ml_, mld_kcyclesym_ml_) - ! Post-smoothing transpose is pre-smoothing + call mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid mltype',& + & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) + goto 9999 + end select - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& - & czero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + call psb_erractionrestore(err_act) + return - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if +9999 call psb_error_handler(err_act) + return + end subroutine inner_ml_aply - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_post - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + recursive subroutine mld_c_inner_add(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod + implicit none - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-cone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vx2l,& - & 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 + !Input/Oputput variables + type(mld_cprec_type), intent(inout) :: p - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_spk_),target :: work(:) + type(psb_c_vect_type) :: res + type(psb_c_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 - end if + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if - case(mld_pre_smooth_) + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& + & czero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%map,info,work=work) - select case (trans_) - case('N') - ! One step of pre-smoothing + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + + sweeps = p%precv(level)%parms%sweeps + 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) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during ADD smoother_apply') + goto 9999 + end if - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& - & czero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + 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 + ! + ! 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 - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_pre - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - - - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-cone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,cone,mlprec_wrk(level)%vx2l,& - & 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 (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & 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 - - - end if - - - case('T','C') - - ! pre-smooth transpose is post-smoothing - - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(cone,mlprec_wrk(level-1)%vx2l,& - & 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 - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& - & czero,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)%vx2l,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 - - - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(cone,& - & mlprec_wrk(level)%vx2l,cone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - else - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') - goto 9999 - end if - end if - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_twoside_smooth_) - - 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 (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 + 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 + call psb_erractionrestore(err_act) + return - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then +9999 call psb_error_handler(err_act) + return - call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& - & czero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) + end subroutine mld_c_inner_add - 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 + recursive subroutine mld_c_inner_mult(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + implicit none + !Input/Oputput variables + type(mld_cprec_type), intent(inout) :: p - ! - ! 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) + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_spk_),target :: work(:) + type(psb_c_vect_type) :: res + type(psb_c_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 + name = 'inner_inner_mult' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_mult') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - endif + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_mult at level ',level + end if - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid smooth_pos',& - & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) - goto 9999 - end select + if ((level < nlev).or.(nlev == 1)) then + sweeps_post = p%precv(level)%parms%sweeps_post + sweeps_pre = p%precv(level)%parms%sweeps_pre + else + sweeps_post = p%precv(level-1)%parms%sweeps_post + sweeps_pre = p%precv(level-1)%parms%sweeps_pre + endif + pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) + post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N')) - 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 + if (level > 1) then + ! Apply the restriction + if (pre) then + current => mlprec_wrk(level-1)%vty + else + current => mlprec_wrk(level-1)%vx2l + endif + call psb_map_X2Y(cone,current,& + & 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 + + + if (level < nlev) then - 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 (pre) then if (trans == 'N') then sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(cone,& @@ -1729,31 +1441,22 @@ contains & 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 - + + 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 + endif ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively + ! Compute the residual and call recursively ! - if(level < nlev) then - + if (pre) 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) @@ -1762,36 +1465,42 @@ contains & 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 - ! + endif + + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + + + ! + ! Apply the prolongator + ! + + if (pre) then + par = cone + else + par = czero + endif + + call psb_map_Y2X(cone,mlprec_wrk(level+1)%vy2l,& + & par,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 + ! + if (post) then + call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& + & czero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) 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) @@ -1803,196 +1512,478 @@ contains ! ! Apply the base preconditioner ! - if (trans == 'N') then + 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 + 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 + + else if (level == nlev) then + + 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) + + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine mld_c_inner_mult + + recursive subroutine mld_c_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + + implicit none + + !Input/Oputput variables + type(mld_cprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_spk_),target :: work(:) + type(psb_c_vect_type),intent(inout), optional :: u + + + + type(psb_c_vect_type) :: res + type(psb_c_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name + + + + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) + + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + !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%zero() + 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_) 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_) 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, u=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 - case(mld_kcycle_ml_, mld_kcyclesym_ml_) + ! + ! 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) - !K cycle + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - call psb_geaxpby(cone,mlprec_wrk(level)%vx2l,& - & czero,mlprec_wrk(level)%vtx,& - & p%precv(level)%base_desc,info) + ! + ! 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 (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 + 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 + 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,& + & 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-PRE smoother_apply') + & a_err='Error during 2-POST smoother_apply') goto 9999 end if + endif - ! - ! 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) + call psb_erractionrestore(err_act) + return - 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 +9999 call psb_error_handler(err_act) + return - ! 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 + end subroutine mld_c_inner_vw_cycle - !Set the preconditioner + recursive subroutine mld_c_inner_k_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + implicit none - 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 + !Input/Oputput variables + type(mld_cprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_spk_),target :: work(:) + type(psb_c_vect_type),intent(inout), optional :: u - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(psb_c_vect_type) :: res + type(psb_c_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - ! - ! 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 + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - 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 + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + !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%zero() 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 - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) - goto 9999 + 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 - end select + 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 call psb_erractionrestore(err_act) return @@ -2000,7 +1991,7 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine inner_ml_aply + end subroutine mld_c_inner_k_cycle recursive subroutine mld_cinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index b0d45972..97932839 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -97,7 +97,7 @@ ! and prolongation operators described in the sequel. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! The integer parameters defining the multilevel -! strategy +! strategy ! p%precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable. ! The real parameters defining the multilevel strategy ! p%precv(ilev)%mlia - integer, dimension(:), allocatable. @@ -333,7 +333,7 @@ subroutine mld_dmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) type mld_mlprec_wrk_type real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_dmlprec_aply' info = psb_success_ @@ -1016,11 +1016,12 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act character(len=20) :: name character :: trans_ + real(psb_dpk_) :: par type mld_mlprec_wrk_type real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_dmlprec_aply' info = psb_success_ @@ -1036,7 +1037,6 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) & ' Entry ', size(p%precv) trans_ = psb_toupper(trans) - nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then @@ -1065,12 +1065,12 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do - level = 1 call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call mlprec_wrk(level)%vy2l%zero() + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -1082,6 +1082,7 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& & p%precv(level)%base_desc,info) do level = 1, nlev + call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vtx%free(info) @@ -1118,29 +1119,30 @@ contains ! Arguments integer(psb_ipk_) :: level type(mld_dprec_type), target, intent(inout) :: p - type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:) character, intent(in) :: trans real(psb_dpk_),target :: work(:) integer(psb_ipk_), intent(out) :: info type(psb_d_vect_type),intent(inout), optional :: u type(psb_d_vect_type) :: res - - + type(psb_d_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre ! Local variables integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: nlev, ilev, sweeps - + logical :: pre, post character(len=20) :: name + + name = 'inner_ml_aply' info = psb_success_ call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& @@ -1167,57 +1169,9 @@ contains goto 9999 case(mld_add_ml_) - ! - ! Additive multilevel - ! - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& - & 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 - - end if - - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& - & done,mlprec_wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if + call mld_d_inner_add(p, mlprec_wrk, level, trans, work) - end if case(mld_mult_ml_) ! @@ -1230,492 +1184,250 @@ contains select case(p%precv(level)%parms%smoother_pos) case(mld_post_smooth_) + p%precv(level)%parms%sweeps_pre = 0 + call mld_d_inner_mult(p, mlprec_wrk, level, trans, work) + - select case (trans_) - case('N') - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& - & dzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + case(mld_pre_smooth_) + p%precv(level)%parms%sweeps_post = 0 + call mld_d_inner_mult(p, mlprec_wrk, level, trans, 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 + case(mld_twoside_smooth_) + call mld_d_inner_mult(p, mlprec_wrk, level, trans, work) - ! This is one step of post-smoothing + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid smooth_pos',& + & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) + goto 9999 - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + end select - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& - & dzero,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)%vx2l,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 + case(mld_mult_dev_ml_) + call mld_d_inner_mult(p, mlprec_wrk, level, trans, work) - sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm2%apply(done,& - & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + case(mld_vcycle_ml_, mld_wcycle_ml_) - end if + call mld_d_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u) - case('T','C') + case(mld_kcycle_ml_, mld_kcyclesym_ml_) - ! Post-smoothing transpose is pre-smoothing + call mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid mltype',& + & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) + goto 9999 + end select - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& - & dzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + call psb_erractionrestore(err_act) + return - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if +9999 call psb_error_handler(err_act) + return + end subroutine inner_ml_aply - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_post - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + recursive subroutine mld_d_inner_add(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod + implicit none - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-done,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,& - & 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 + !Input/Oputput variables + type(mld_dprec_type), intent(inout) :: p - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_dpk_),target :: work(:) + type(psb_d_vect_type) :: res + type(psb_d_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 - end if + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if - case(mld_pre_smooth_) + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& + & dzero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%map,info,work=work) - select case (trans_) - case('N') - ! One step of pre-smoothing + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + + sweeps = p%precv(level)%parms%sweeps + 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) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during ADD smoother_apply') + goto 9999 + end if - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& - & dzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + 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 + ! + ! 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 - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_pre - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - - - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-done,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vx2l,& - & 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 (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & 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 - - - end if - - - case('T','C') - - ! pre-smooth transpose is post-smoothing - - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(done,mlprec_wrk(level-1)%vx2l,& - & 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 - - end if - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& - & dzero,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)%vx2l,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 - - - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vx2l,done,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - else - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') - goto 9999 - end if - end if - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_twoside_smooth_) - - 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 - end if - - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vtx,& - & p%precv(level)%base_desc,info) - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - if (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 + 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 + call psb_erractionrestore(err_act) + return - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then +9999 call psb_error_handler(err_act) + return - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) + end subroutine mld_d_inner_add - 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 + recursive subroutine mld_d_inner_mult(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + implicit none + !Input/Oputput variables + type(mld_dprec_type), intent(inout) :: p - ! - ! 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) + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_dpk_),target :: work(:) + type(psb_d_vect_type) :: res + type(psb_d_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 + name = 'inner_inner_mult' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_mult') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - endif + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_mult at level ',level + end if - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid smooth_pos',& - & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) - goto 9999 - end select + if ((level < nlev).or.(nlev == 1)) then + sweeps_post = p%precv(level)%parms%sweeps_post + sweeps_pre = p%precv(level)%parms%sweeps_pre + else + sweeps_post = p%precv(level-1)%parms%sweeps_post + sweeps_pre = p%precv(level-1)%parms%sweeps_pre + endif + pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) + post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N')) - 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) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if + if (level > 1) then + ! Apply the restriction + if (pre) then + current => mlprec_wrk(level-1)%vty + else + current => mlprec_wrk(level-1)%vx2l + endif + call psb_map_X2Y(done,current,& + & 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 + end if + + + 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()) - call psb_geaxpby(done,u,& - & dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc,info) - - else - call mlprec_wrk(level)%vy2l%set(dzero) - endif - res = mlprec_wrk(level)%vx2l - - call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& - done, res, p%precv(level)%base_desc, info, work=work, trans=trans) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if + if (pre) then if (trans == 'N') then sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(done,& @@ -1729,31 +1441,22 @@ contains & 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 - + + 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 + endif ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively + ! Compute the residual and call recursively ! - if(level < nlev) then - + if (pre) 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) @@ -1762,36 +1465,42 @@ contains & a_err='Error during residue') goto 9999 end if - - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - - if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, mlprec_wrk(level+1)%vy2l) - endif - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& - & done,mlprec_wrk(level)%vy2l,& - & p%precv(level+1)%map,info,work=work) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during prolongation') - goto 9999 - end if - - ! - ! Compute the residual - ! + endif + + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + + + ! + ! Apply the prolongator + ! + + if (pre) then + par = done + else + par = dzero + endif + + call psb_map_Y2X(done,mlprec_wrk(level+1)%vy2l,& + & par,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 + ! + if (post) then + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) 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) @@ -1803,196 +1512,478 @@ contains ! ! Apply the base preconditioner ! - if (trans == 'N') then + 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 + 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 + + else if (level == nlev) then + + 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) + + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine mld_d_inner_mult + + recursive subroutine mld_d_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + + implicit none + + !Input/Oputput variables + type(mld_dprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_dpk_),target :: work(:) + type(psb_d_vect_type),intent(inout), optional :: u + + + + type(psb_d_vect_type) :: res + type(psb_d_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name + + + + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) + + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + !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 + end if + + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) + ! + ! Apply the base preconditioner + ! + if (level < nlev) then + + if (present(u)) then + ! call mlprec_wrk(level)%vy2l%set(u%get_vect()) + call psb_geaxpby(done,u,& + & dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc,info) + + else + call mlprec_wrk(level)%vy2l%zero() + endif + res = mlprec_wrk(level)%vx2l + + call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& + done, res, p%precv(level)%base_desc, info, work=work, trans=trans) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if + + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + else + sweeps = p%precv(level)%parms%sweeps + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if + + + ! + ! Compute the residual (at all levels but the coarsest one) + ! and call recursively + ! + if(level < nlev) then + + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) - 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_) call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + + if (p%precv(level)%parms%ml_type == mld_wcycle_ml_) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info, u=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 - case(mld_kcycle_ml_, mld_kcyclesym_ml_) + ! + ! 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) - !K cycle + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& - & dzero,mlprec_wrk(level)%vtx,& - & p%precv(level)%base_desc,info) + ! + ! Compute the residual + ! + call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& + & done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& + & work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if ! ! Apply the base preconditioner ! - if (level < nlev) then - - if (present(u)) then - call psb_geaxpby(done,u,& - & dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc,info) - else - call mlprec_wrk(level)%vy2l%set(dzero) - endif - res = mlprec_wrk(level)%vx2l - - call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& - done, res, p%precv(level)%base_desc, info, work=work, trans=trans) - - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during residue') - goto 9999 - end if - - if (trans == 'N') then - sweeps = p%precv(level)%parms%sweeps_pre - if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - else - sweeps = p%precv(level)%parms%sweeps_post - if (info == psb_success_) call p%precv(level)%sm2%apply(done,& - & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - end if + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) else - sweeps = p%precv(level)%parms%sweeps + sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(done,& - & mlprec_wrk(level)%vx2l,dzero,mlprec_wrk(level)%vy2l,& + & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,& & p%precv(level)%base_desc, trans,& & sweeps,work,info) end if if (info /= psb_success_) then call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during 2-PRE smoother_apply') + & a_err='Error during 2-POST smoother_apply') goto 9999 end if + endif - ! - ! 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) + call psb_erractionrestore(err_act) + return - 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 +9999 call psb_error_handler(err_act) + return - ! 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 + end subroutine mld_d_inner_vw_cycle - !Set the preconditioner + recursive subroutine mld_d_inner_k_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + implicit none - 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 + !Input/Oputput variables + type(mld_dprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_dpk_),target :: work(:) + type(psb_d_vect_type),intent(inout), optional :: u - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(psb_d_vect_type) :: res + type(psb_d_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - ! - ! 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 + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - 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 + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + !K cycle + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) + ! + ! Apply the base preconditioner + ! + if (level < nlev) then + + if (present(u)) then + call psb_geaxpby(done,u,& + & dzero,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc,info) + else + call mlprec_wrk(level)%vy2l%zero() endif + res = mlprec_wrk(level)%vx2l + call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& + done, res, p%precv(level)%base_desc, info, work=work, trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) - goto 9999 + 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 - end select + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-PRE smoother_apply') + goto 9999 + end if + + + ! + ! Compute the residual (at all levels but the coarsest one) + ! and call recursively + ! + if(level < nlev) then + + call psb_geaxpby(done,mlprec_wrk(level)%vx2l,& + & dzero,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info) + + if (info == psb_success_) call psb_spmm(-done,p%precv(level)%base_a,& + & mlprec_wrk(level)%vy2l,done,mlprec_wrk(level)%vty,& + & p%precv(level)%base_desc,info,work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if + + ! Apply the restriction + call psb_map_X2Y(done,mlprec_wrk(level)%vty,& + & dzero,mlprec_wrk(level + 1)%vx2l,& + & p%precv(level + 1)%map,info,work=work) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + + !Set the preconditioner + + + if ((level < nlev - 2)) then + if (p%precv(level)%parms%ml_type == mld_kcyclesym_ml_) then + call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'FCG') + elseif (p%precv(level)%parms%ml_type == mld_kcycle_ml_) then + call mld_dinneritkcycle(p, mlprec_wrk, level + 1, trans, work, 'CGR') + endif + 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(done,mlprec_wrk(level+1)%vy2l,& + & done,mlprec_wrk(level)%vy2l,& + & p%precv(level+1)%map,info,work=work) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if + + ! + ! Compute the residual + ! + call psb_spmm(-done,p%precv(level)%base_a,mlprec_wrk(level)%vy2l,& + & done,mlprec_wrk(level)%vtx,p%precv(level)%base_desc,info,& + & work=work,trans=trans) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during residue') + goto 9999 + end if + ! + ! Apply the base preconditioner + ! + if (trans == 'N') then + sweeps = p%precv(level)%parms%sweeps_post + if (info == psb_success_) call p%precv(level)%sm2%apply(done,& + & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + else + sweeps = p%precv(level)%parms%sweeps_pre + if (info == psb_success_) call p%precv(level)%sm%apply(done,& + & mlprec_wrk(level)%vtx,done,mlprec_wrk(level)%vy2l,& + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) + end if + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during 2-POST smoother_apply') + goto 9999 + end if + + endif call psb_erractionrestore(err_act) return @@ -2000,7 +1991,7 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine inner_ml_aply + end subroutine mld_d_inner_k_cycle recursive subroutine mld_dinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 424e5861..26ab65a7 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -84,7 +84,7 @@ ! The multilevel preconditioner data structure containing the ! local part of the preconditioner to be applied. ! Note that nlev = size(p%precv) = number of levels. -! p%precv(ilev)%prec - type(psb_dbaseprec_type) +! p%precv(ilev)%prec - type(psb_sbaseprec_type) ! The 'base preconditioner' for the current level ! p%precv(ilev)%ac - type(psb_sspmat_type) ! The local part of the matrix A(ilev). @@ -97,7 +97,7 @@ ! and prolongation operators described in the sequel. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! The integer parameters defining the multilevel -! strategy +! strategy ! p%precv(ilev)%rprcparm - real(psb_spk_), dimension(:), allocatable. ! The real parameters defining the multilevel strategy ! p%precv(ilev)%mlia - integer, dimension(:), allocatable. @@ -333,7 +333,7 @@ subroutine mld_smlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) type mld_mlprec_wrk_type real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_smlprec_aply' info = psb_success_ @@ -1016,11 +1016,12 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act character(len=20) :: name character :: trans_ + real(psb_spk_) :: par type mld_mlprec_wrk_type real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_smlprec_aply' info = psb_success_ @@ -1036,7 +1037,6 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) & ' Entry ', size(p%precv) trans_ = psb_toupper(trans) - nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then @@ -1065,12 +1065,12 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do - level = 1 call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call mlprec_wrk(level)%vy2l%zero() + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -1082,6 +1082,7 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& & p%precv(level)%base_desc,info) do level = 1, nlev + call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vtx%free(info) @@ -1118,29 +1119,30 @@ contains ! Arguments integer(psb_ipk_) :: level type(mld_sprec_type), target, intent(inout) :: p - type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:) character, intent(in) :: trans real(psb_spk_),target :: work(:) integer(psb_ipk_), intent(out) :: info type(psb_s_vect_type),intent(inout), optional :: u type(psb_s_vect_type) :: res - - + type(psb_s_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre ! Local variables integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: nlev, ilev, sweeps - + logical :: pre, post character(len=20) :: name + + name = 'inner_ml_aply' info = psb_success_ call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& @@ -1167,57 +1169,9 @@ contains goto 9999 case(mld_add_ml_) - ! - ! Additive multilevel - ! - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& - & 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 - - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error 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 + call mld_s_inner_add(p, mlprec_wrk, level, trans, work) - end if case(mld_mult_ml_) ! @@ -1230,492 +1184,250 @@ contains select case(p%precv(level)%parms%smoother_pos) case(mld_post_smooth_) + p%precv(level)%parms%sweeps_pre = 0 + call mld_s_inner_mult(p, mlprec_wrk, level, trans, work) + - select case (trans_) - case('N') - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& - & szero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + case(mld_pre_smooth_) + p%precv(level)%parms%sweeps_post = 0 + call mld_s_inner_mult(p, mlprec_wrk, level, trans, 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 + case(mld_twoside_smooth_) + call mld_s_inner_mult(p, mlprec_wrk, level, trans, work) - ! This is one step of post-smoothing + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid smooth_pos',& + & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) + goto 9999 - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + end select - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& - & szero,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)%vx2l,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 + case(mld_mult_dev_ml_) + call mld_s_inner_mult(p, mlprec_wrk, level, trans, work) - sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm2%apply(sone,& - & mlprec_wrk(level)%vx2l,sone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & mlprec_wrk(level)%vx2l,szero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + case(mld_vcycle_ml_, mld_wcycle_ml_) - end if + call mld_s_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u) - case('T','C') + case(mld_kcycle_ml_, mld_kcyclesym_ml_) - ! Post-smoothing transpose is pre-smoothing + call mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid mltype',& + & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) + goto 9999 + end select - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& - & szero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + call psb_erractionrestore(err_act) + return - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if +9999 call psb_error_handler(err_act) + return + end subroutine inner_ml_aply - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_post - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + recursive subroutine mld_s_inner_add(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod + implicit none - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-sone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vx2l,& - & 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 + !Input/Oputput variables + type(mld_sprec_type), intent(inout) :: p - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_spk_),target :: work(:) + type(psb_s_vect_type) :: res + type(psb_s_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 - end if + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if - case(mld_pre_smooth_) + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& + & szero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%map,info,work=work) - select case (trans_) - case('N') - ! One step of pre-smoothing + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + + sweeps = p%precv(level)%parms%sweeps + 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) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during ADD smoother_apply') + goto 9999 + end if - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& - & szero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + 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 + ! + ! 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 - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_pre - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - - - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-sone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,sone,mlprec_wrk(level)%vx2l,& - & 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 (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & 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 - - - end if - - - case('T','C') - - ! pre-smooth transpose is post-smoothing - - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(sone,mlprec_wrk(level-1)%vx2l,& - & 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 - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& - & szero,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)%vx2l,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 - - - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(sone,& - & mlprec_wrk(level)%vx2l,sone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - else - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') - goto 9999 - end if - end if - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_twoside_smooth_) - - 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 (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 + 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 + call psb_erractionrestore(err_act) + return - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then +9999 call psb_error_handler(err_act) + return - call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& - & szero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) + end subroutine mld_s_inner_add - 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 + recursive subroutine mld_s_inner_mult(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + implicit none + !Input/Oputput variables + type(mld_sprec_type), intent(inout) :: p - ! - ! 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) + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_spk_),target :: work(:) + type(psb_s_vect_type) :: res + type(psb_s_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 + name = 'inner_inner_mult' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_mult') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - endif + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_mult at level ',level + end if - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid smooth_pos',& - & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) - goto 9999 - end select + if ((level < nlev).or.(nlev == 1)) then + sweeps_post = p%precv(level)%parms%sweeps_post + sweeps_pre = p%precv(level)%parms%sweeps_pre + else + sweeps_post = p%precv(level-1)%parms%sweeps_post + sweeps_pre = p%precv(level-1)%parms%sweeps_pre + endif + pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) + post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N')) - 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 + if (level > 1) then + ! Apply the restriction + if (pre) then + current => mlprec_wrk(level-1)%vty + else + current => mlprec_wrk(level-1)%vx2l + endif + call psb_map_X2Y(sone,current,& + & 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 + + + if (level < nlev) then - 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 (pre) then if (trans == 'N') then sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(sone,& @@ -1729,31 +1441,22 @@ contains & 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 - + + 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 + endif ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively + ! Compute the residual and call recursively ! - if(level < nlev) then - + if (pre) 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) @@ -1762,36 +1465,42 @@ contains & 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 - ! + endif + + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + + + ! + ! Apply the prolongator + ! + + if (pre) then + par = sone + else + par = szero + endif + + call psb_map_Y2X(sone,mlprec_wrk(level+1)%vy2l,& + & par,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 + ! + if (post) then + call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& + & szero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) 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) @@ -1803,196 +1512,478 @@ contains ! ! Apply the base preconditioner ! - if (trans == 'N') then + 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 + 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 + + else if (level == nlev) then + + 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) + + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine mld_s_inner_mult + + recursive subroutine mld_s_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + + implicit none + + !Input/Oputput variables + type(mld_sprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_spk_),target :: work(:) + type(psb_s_vect_type),intent(inout), optional :: u + + + + type(psb_s_vect_type) :: res + type(psb_s_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name + + + + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) + + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + !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%zero() + 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_) 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_) 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, u=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 - case(mld_kcycle_ml_, mld_kcyclesym_ml_) + ! + ! 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) - !K cycle + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - call psb_geaxpby(sone,mlprec_wrk(level)%vx2l,& - & szero,mlprec_wrk(level)%vtx,& - & p%precv(level)%base_desc,info) + ! + ! 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 (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 + 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 + 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,& + & 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-PRE smoother_apply') + & a_err='Error during 2-POST smoother_apply') goto 9999 end if + endif - ! - ! 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) + call psb_erractionrestore(err_act) + return - 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 +9999 call psb_error_handler(err_act) + return - ! 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 + end subroutine mld_s_inner_vw_cycle - !Set the preconditioner + recursive subroutine mld_s_inner_k_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + implicit none - 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 + !Input/Oputput variables + type(mld_sprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + real(psb_spk_),target :: work(:) + type(psb_s_vect_type),intent(inout), optional :: u - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(psb_s_vect_type) :: res + type(psb_s_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - ! - ! 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 + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - 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 + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + !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%zero() 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 - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) - goto 9999 + 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 - end select + 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 call psb_erractionrestore(err_act) return @@ -2000,7 +1991,7 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine inner_ml_aply + end subroutine mld_s_inner_k_cycle recursive subroutine mld_sinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index 56c83d16..6afe3609 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -39,7 +39,7 @@ ! File: mld_zmlprec_aply.f90 ! ! Subroutine: mld_zmlprec_aply -! Version: complex +! Version: real ! ! This routine computes ! @@ -84,7 +84,7 @@ ! The multilevel preconditioner data structure containing the ! local part of the preconditioner to be applied. ! Note that nlev = size(p%precv) = number of levels. -! p%precv(ilev)%prec - type(psb_dbaseprec_type) +! p%precv(ilev)%prec - type(psb_zbaseprec_type) ! The 'base preconditioner' for the current level ! p%precv(ilev)%ac - type(psb_zspmat_type) ! The local part of the matrix A(ilev). @@ -97,7 +97,7 @@ ! and prolongation operators described in the sequel. ! p%precv(ilev)%iprcparm - integer, dimension(:), allocatable. ! The integer parameters defining the multilevel -! strategy +! strategy ! p%precv(ilev)%rprcparm - real(psb_dpk_), dimension(:), allocatable. ! The real parameters defining the multilevel strategy ! p%precv(ilev)%mlia - integer, dimension(:), allocatable. @@ -333,7 +333,7 @@ subroutine mld_zmlprec_aply(alpha,p,x,beta,y,desc_data,trans,work,info) type mld_mlprec_wrk_type complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_zmlprec_aply' info = psb_success_ @@ -1016,11 +1016,12 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act character(len=20) :: name character :: trans_ + complex(psb_dpk_) :: par type mld_mlprec_wrk_type complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l end type mld_mlprec_wrk_type - type(mld_mlprec_wrk_type), allocatable :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), allocatable, target :: mlprec_wrk(:) name='mld_zmlprec_aply' info = psb_success_ @@ -1036,7 +1037,6 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) & ' Entry ', size(p%precv) trans_ = psb_toupper(trans) - nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then @@ -1065,12 +1065,12 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do - level = 1 call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) call mlprec_wrk(level)%vy2l%zero() + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -1082,6 +1082,7 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& & p%precv(level)%base_desc,info) do level = 1, nlev + call mlprec_wrk(level)%vx2l%free(info) call mlprec_wrk(level)%vy2l%free(info) call mlprec_wrk(level)%vtx%free(info) @@ -1118,29 +1119,30 @@ contains ! Arguments integer(psb_ipk_) :: level type(mld_zprec_type), target, intent(inout) :: p - type(mld_mlprec_wrk_type), intent(inout) :: mlprec_wrk(:) + type(mld_mlprec_wrk_type), intent(inout), target :: mlprec_wrk(:) character, intent(in) :: trans complex(psb_dpk_),target :: work(:) integer(psb_ipk_), intent(out) :: info type(psb_z_vect_type),intent(inout), optional :: u type(psb_z_vect_type) :: res - - + type(psb_z_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre ! Local variables integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: i, nr2l,nc2l,err_act integer(psb_ipk_) :: debug_level, debug_unit integer(psb_ipk_) :: nlev, ilev, sweeps - + logical :: pre, post character(len=20) :: name + + name = 'inner_ml_aply' info = psb_success_ call psb_erractionsave(err_act) debug_unit = psb_get_debug_unit() debug_level = psb_get_debug_level() - nlev = size(p%precv) if ((level < 1) .or. (level > nlev)) then call psb_errpush(psb_err_internal_error_,name,& @@ -1167,57 +1169,9 @@ contains goto 9999 case(mld_add_ml_) - ! - ! Additive multilevel - ! - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& - & 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 - - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during ADD smoother_apply') - goto 9999 - end if - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error 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 + call mld_z_inner_add(p, mlprec_wrk, level, trans, work) - end if case(mld_mult_ml_) ! @@ -1230,492 +1184,250 @@ contains select case(p%precv(level)%parms%smoother_pos) case(mld_post_smooth_) + p%precv(level)%parms%sweeps_pre = 0 + call mld_z_inner_mult(p, mlprec_wrk, level, trans, work) + - select case (trans_) - case('N') - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& - & zzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + case(mld_pre_smooth_) + p%precv(level)%parms%sweeps_post = 0 + call mld_z_inner_mult(p, mlprec_wrk, level, trans, 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 + case(mld_twoside_smooth_) + call mld_z_inner_mult(p, mlprec_wrk, level, trans, work) - ! This is one step of post-smoothing + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid smooth_pos',& + & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) + goto 9999 - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + end select - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& - & zzero,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)%vx2l,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 + case(mld_mult_dev_ml_) + call mld_z_inner_mult(p, mlprec_wrk, level, trans, work) - sweeps = p%precv(level)%parms%sweeps_post - call p%precv(level)%sm2%apply(zone,& - & mlprec_wrk(level)%vx2l,zone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & mlprec_wrk(level)%vx2l,zzero,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + case(mld_vcycle_ml_, mld_wcycle_ml_) - end if + call mld_z_inner_vw_cycle(p, mlprec_wrk, level, trans, work, u=u) - case('T','C') + case(mld_kcycle_ml_, mld_kcyclesym_ml_) - ! Post-smoothing transpose is pre-smoothing + call mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work, u=u) + + case default + info = psb_err_from_subroutine_ai_ + call psb_errpush(info,name,a_err='invalid mltype',& + & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) + goto 9999 + end select - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& - & zzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + call psb_erractionrestore(err_act) + return - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during restriction') - goto 9999 - end if +9999 call psb_error_handler(err_act) + return + end subroutine inner_ml_aply - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_post - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during POST smoother_apply') - goto 9999 - end if + recursive subroutine mld_z_inner_add(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod + implicit none - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-zone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vx2l,& - & 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 + !Input/Oputput variables + type(mld_zprec_type), intent(inout) :: p - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_dpk_),target :: work(:) + type(psb_z_vect_type) :: res + type(psb_z_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 - end if + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if - case(mld_pre_smooth_) + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + if (level > 1) then + ! Apply the restriction + call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& + & zzero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%map,info,work=work) - select case (trans_) - case('N') - ! One step of pre-smoothing + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during restriction') + goto 9999 + end if + end if + + sweeps = p%precv(level)%parms%sweeps + 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) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during ADD smoother_apply') + goto 9999 + end if - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& - & zzero,mlprec_wrk(level)%vx2l,& - & p%precv(level)%map,info,work=work) + if (level < nlev) then + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + 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 + ! + ! 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 - end if - ! - ! Apply the base preconditioner - ! - if (level < nlev) then - sweeps = p%precv(level)%parms%sweeps_pre - else - sweeps = p%precv(level)%parms%sweeps - end if - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - - - ! - ! Compute the residual (at all levels but the coarsest one) - ! - if (level < nlev) then - call psb_spmm(-zone,p%precv(level)%base_a,& - & mlprec_wrk(level)%vy2l,zone,mlprec_wrk(level)%vx2l,& - & 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 (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & 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,& - & 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 - - - end if - - - case('T','C') - - ! pre-smooth transpose is post-smoothing - - - if (level > 1) then - ! Apply the restriction - call psb_map_X2Y(zone,mlprec_wrk(level-1)%vx2l,& - & 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 - - if (level < nlev) then - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if - - ! - ! Apply the prolongator - ! - call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& - & zzero,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)%vx2l,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 - - - sweeps = p%precv(level)%parms%sweeps_pre - call p%precv(level)%sm%apply(zone,& - & mlprec_wrk(level)%vx2l,zone,mlprec_wrk(level)%vy2l,& - & p%precv(level)%base_desc, trans,& - & sweeps,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during PRE smoother_apply') - goto 9999 - end if - else - sweeps = p%precv(level)%parms%sweeps - 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) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error during smoother_apply') - goto 9999 - end if - end if - - case default - info = psb_err_internal_error_ - call psb_errpush(info,name,a_err='invalid trans') - goto 9999 - end select - - case(mld_twoside_smooth_) - - 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 (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 + 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 + call psb_erractionrestore(err_act) + return - ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively - ! - if(level < nlev) then +9999 call psb_error_handler(err_act) + return - call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& - & zzero,mlprec_wrk(level)%vty,& - & p%precv(level)%base_desc,info) + end subroutine mld_z_inner_add - 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 + recursive subroutine mld_z_inner_mult(p, mlprec_wrk, level, trans, work) + use psb_base_mod + use mld_prec_mod - call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + implicit none + !Input/Oputput variables + type(mld_zprec_type), intent(inout) :: p - ! - ! 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) + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_dpk_),target :: work(:) + type(psb_z_vect_type) :: res + type(psb_z_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - 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 + name = 'inner_inner_mult' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_mult') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - endif + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_mult at level ',level + end if - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid smooth_pos',& - & i_Err=(/p%precv(level)%parms%smoother_pos,izero,izero,izero,izero/)) - goto 9999 - end select + if ((level < nlev).or.(nlev == 1)) then + sweeps_post = p%precv(level)%parms%sweeps_post + sweeps_pre = p%precv(level)%parms%sweeps_pre + else + sweeps_post = p%precv(level-1)%parms%sweeps_post + sweeps_pre = p%precv(level-1)%parms%sweeps_pre + endif + pre = ((sweeps_pre>0).and.(trans=='N')).or.((sweeps_post>0).and.(trans/='N')) + post = ((sweeps_post>0).and.(trans=='N')).or.((sweeps_pre>0).and.(trans/='N')) - 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 + if (level > 1) then + ! Apply the restriction + if (pre) then + current => mlprec_wrk(level-1)%vty + else + current => mlprec_wrk(level-1)%vx2l + endif + call psb_map_X2Y(zone,current,& + & 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 + + + if (level < nlev) then - 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 (pre) then if (trans == 'N') then sweeps = p%precv(level)%parms%sweeps_pre if (info == psb_success_) call p%precv(level)%sm%apply(zone,& @@ -1729,31 +1441,22 @@ contains & 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 - + + 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 + endif ! - ! Compute the residual (at all levels but the coarsest one) - ! and call recursively + ! Compute the residual and call recursively ! - if(level < nlev) then - + if (pre) 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) @@ -1762,36 +1465,42 @@ contains & 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 - ! + endif + + call inner_ml_aply(level+1,p,mlprec_wrk,trans,work,info) + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error in recursive call') + goto 9999 + end if + + + ! + ! Apply the prolongator + ! + + if (pre) then + par = zone + else + par = zzero + endif + + call psb_map_Y2X(zone,mlprec_wrk(level+1)%vy2l,& + & par,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 + ! + if (post) then + call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& + & zzero,mlprec_wrk(level)%vtx,& + & p%precv(level)%base_desc,info) 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) @@ -1803,196 +1512,478 @@ contains ! ! Apply the base preconditioner ! - if (trans == 'N') then + 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 + 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 + + else if (level == nlev) then + + 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) + + else + + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 call psb_error_handler(err_act) + return + + end subroutine mld_z_inner_mult + + recursive subroutine mld_z_inner_vw_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + + implicit none + + !Input/Oputput variables + type(mld_zprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_dpk_),target :: work(:) + type(psb_z_vect_type),intent(inout), optional :: u + + + + type(psb_z_vect_type) :: res + type(psb_z_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name + + + + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) + + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + + !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%zero() + 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_) 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_) 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, u=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 - case(mld_kcycle_ml_, mld_kcyclesym_ml_) + ! + ! 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) - !K cycle + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Error during prolongation') + goto 9999 + end if - call psb_geaxpby(zone,mlprec_wrk(level)%vx2l,& - & zzero,mlprec_wrk(level)%vtx,& - & p%precv(level)%base_desc,info) + ! + ! 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 (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 + 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 + 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,& + & 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-PRE smoother_apply') + & a_err='Error during 2-POST smoother_apply') goto 9999 end if + endif - ! - ! 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) + call psb_erractionrestore(err_act) + return - 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 +9999 call psb_error_handler(err_act) + return - ! 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 + end subroutine mld_z_inner_vw_cycle - !Set the preconditioner + recursive subroutine mld_z_inner_k_cycle(p, mlprec_wrk, level, trans, work,u) + use psb_base_mod + use mld_prec_mod + implicit none - 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 + !Input/Oputput variables + type(mld_zprec_type), intent(inout) :: p + type(mld_mlprec_wrk_type), target, intent(inout) :: mlprec_wrk(:) + integer(psb_ipk_), intent(in) :: level + character, intent(in) :: trans + complex(psb_dpk_),target :: work(:) + type(psb_z_vect_type),intent(inout), optional :: u - if (info /= psb_success_) then - call psb_errpush(psb_err_internal_error_,name,& - & a_err='Error in recursive call') - goto 9999 - end if + type(psb_z_vect_type) :: res + type(psb_z_vect_type), pointer :: current + integer(psb_ipk_) :: sweeps_post, sweeps_pre + ! Local variables + integer(psb_ipk_) :: ictxt,np,me + integer(psb_ipk_) :: i, nr2l,nc2l,err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev, ilev, sweeps + logical :: pre, post + character(len=20) :: name - ! - ! 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 + name = 'inner_inner_add' + info = psb_success_ + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + nlev = size(p%precv) + if ((level < 1) .or. (level > nlev)) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='wrong call level to inner_add') + goto 9999 + end if + ictxt = p%precv(level)%base_desc%get_context() + call psb_info(ictxt, me, np) - 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 + nc2l = p%precv(level)%base_desc%get_local_cols() + nr2l = p%precv(level)%base_desc%get_local_rows() + if(debug_level > 1) then + write(debug_unit,*) me,' inner_add at level ',level + end if + + if ((level<1).or.(level>nlev)) then + info = psb_err_internal_error_ + call psb_errpush(info,name,& + & a_err='Invalid LEVEL>NLEV') + goto 9999 + end if + + !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%zero() 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 - case default - info = psb_err_from_subroutine_ai_ - call psb_errpush(info,name,a_err='invalid mltype',& - & i_Err=(/p%precv(level)%parms%ml_type,izero,izero,izero,izero/)) - goto 9999 + 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 - end select + 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 call psb_erractionrestore(err_act) return @@ -2000,7 +1991,7 @@ contains 9999 call psb_error_handler(err_act) return - end subroutine inner_ml_aply + end subroutine mld_z_inner_k_cycle recursive subroutine mld_zinneritkcycle(p, mlprec_wrk, level, trans, work, innersolv) diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 17c76645..4b3b66df 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -227,6 +227,7 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_kcycle_ml_ = 5 integer(psb_ipk_), parameter :: mld_kcyclesym_ml_ = 6 integer(psb_ipk_), parameter :: mld_new_ml_prec_ = 7 + integer(psb_ipk_), parameter :: mld_mult_dev_ml_ = 7 integer(psb_ipk_), parameter :: mld_max_ml_type_ = 8 ! ! Legal values for entry: mld_smoother_pos_ @@ -424,6 +425,8 @@ contains val = mld_diag_scale_ case('ADD') val = mld_add_ml_ + case('MULT_DEV') + val = mld_mult_dev_ml_ case('MULT') val = mld_mult_ml_ case('VCYCLE')