diff --git a/mlprec/impl/mld_ccprecset.F90 b/mlprec/impl/mld_ccprecset.F90 index 6f627992..783b4582 100644 --- a/mlprec/impl/mld_ccprecset.F90 +++ b/mlprec/impl/mld_ccprecset.F90 @@ -143,6 +143,9 @@ subroutine mld_ccprecseti(p,what,val,info,ilev,pos) case('MAX_PREC_LEVS') p%max_prec_levs = max(val,1) return + case ('OUTER_SWEEPS') + p%outer_sweeps = max(val,1) + return end select diff --git a/mlprec/impl/mld_cmlprec_aply.f90 b/mlprec/impl/mld_cmlprec_aply.f90 index 79c6c3ff..03aaaf1a 100644 --- a/mlprec/impl/mld_cmlprec_aply.f90 +++ b/mlprec/impl/mld_cmlprec_aply.f90 @@ -314,10 +314,11 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) ! Local variables integer(psb_ipk_) :: ictxt, np, me - integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act character(len=20) :: name character :: trans_ - complex(psb_spk_) :: par + complex(psb_spk_) :: beta_ type mld_mlprec_wrk_type complex(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_c_vect_type) :: vtx, vty, vx2l, vy2l @@ -338,7 +339,7 @@ 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) + nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') @@ -366,12 +367,49 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do + ! + ! At first iteration we must use the input BETA + ! + beta_ = beta + level = 1 call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) - call mlprec_wrk(level)%vy2l%zero() + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,& + & p%precv(level)%base_desc,info) + ! all iterations after the first must use BETA = 1 + beta_ = cone + ! + ! Next iteration should use the current residual to compute a correction + ! + call psb_geaxpby(cone,x,czero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%base_desc,info) + call psb_spmm(-cone,p%precv(level)%base_a,y,& + & cone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) + end do + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -380,8 +418,9 @@ subroutine mld_cmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if - call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& + 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) diff --git a/mlprec/impl/mld_cprecset.F90 b/mlprec/impl/mld_cprecset.F90 index 4eefe7f0..38c1d791 100644 --- a/mlprec/impl/mld_cprecset.F90 +++ b/mlprec/impl/mld_cprecset.F90 @@ -142,6 +142,9 @@ subroutine mld_cprecseti(p,what,val,info,ilev,pos) case(mld_max_prec_levs_) p%max_prec_levs = max(val,1) return + case(mld_outer_sweeps_) + p%outer_sweeps = max(val,1) + return end select ! diff --git a/mlprec/impl/mld_dcprecset.F90 b/mlprec/impl/mld_dcprecset.F90 index feff56e4..bb61822a 100644 --- a/mlprec/impl/mld_dcprecset.F90 +++ b/mlprec/impl/mld_dcprecset.F90 @@ -149,6 +149,9 @@ subroutine mld_dcprecseti(p,what,val,info,ilev,pos) case('MAX_PREC_LEVS') p%max_prec_levs = max(val,1) return + case ('OUTER_SWEEPS') + p%outer_sweeps = max(val,1) + return end select diff --git a/mlprec/impl/mld_dmlprec_aply.f90 b/mlprec/impl/mld_dmlprec_aply.f90 index ef860cb4..00d432f8 100644 --- a/mlprec/impl/mld_dmlprec_aply.f90 +++ b/mlprec/impl/mld_dmlprec_aply.f90 @@ -314,10 +314,11 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) ! Local variables integer(psb_ipk_) :: ictxt, np, me - integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act character(len=20) :: name character :: trans_ - real(psb_dpk_) :: par + real(psb_dpk_) :: beta_ type mld_mlprec_wrk_type real(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_d_vect_type) :: vtx, vty, vx2l, vy2l @@ -338,7 +339,7 @@ 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) + nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') @@ -366,12 +367,49 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do + ! + ! At first iteration we must use the input BETA + ! + beta_ = beta + level = 1 call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) - call mlprec_wrk(level)%vy2l%zero() + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,& + & p%precv(level)%base_desc,info) + ! all iterations after the first must use BETA = 1 + beta_ = done + ! + ! Next iteration should use the current residual to compute a correction + ! + call psb_geaxpby(done,x,dzero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%base_desc,info) + call psb_spmm(-done,p%precv(level)%base_a,y,& + & done,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) + end do + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -380,8 +418,9 @@ subroutine mld_dmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if - call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& + 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) diff --git a/mlprec/impl/mld_dprecset.F90 b/mlprec/impl/mld_dprecset.F90 index 4991806c..37236b03 100644 --- a/mlprec/impl/mld_dprecset.F90 +++ b/mlprec/impl/mld_dprecset.F90 @@ -148,6 +148,9 @@ subroutine mld_dprecseti(p,what,val,info,ilev,pos) case(mld_max_prec_levs_) p%max_prec_levs = max(val,1) return + case(mld_outer_sweeps_) + p%outer_sweeps = max(val,1) + return end select ! diff --git a/mlprec/impl/mld_scprecset.F90 b/mlprec/impl/mld_scprecset.F90 index c39839fb..af641414 100644 --- a/mlprec/impl/mld_scprecset.F90 +++ b/mlprec/impl/mld_scprecset.F90 @@ -143,6 +143,9 @@ subroutine mld_scprecseti(p,what,val,info,ilev,pos) case('MAX_PREC_LEVS') p%max_prec_levs = max(val,1) return + case ('OUTER_SWEEPS') + p%outer_sweeps = max(val,1) + return end select diff --git a/mlprec/impl/mld_smlprec_aply.f90 b/mlprec/impl/mld_smlprec_aply.f90 index 4a24c246..1c57eec3 100644 --- a/mlprec/impl/mld_smlprec_aply.f90 +++ b/mlprec/impl/mld_smlprec_aply.f90 @@ -314,10 +314,11 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) ! Local variables integer(psb_ipk_) :: ictxt, np, me - integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act character(len=20) :: name character :: trans_ - real(psb_spk_) :: par + real(psb_spk_) :: beta_ type mld_mlprec_wrk_type real(psb_spk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_s_vect_type) :: vtx, vty, vx2l, vy2l @@ -338,7 +339,7 @@ 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) + nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') @@ -366,12 +367,49 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do + ! + ! At first iteration we must use the input BETA + ! + beta_ = beta + level = 1 call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) - call mlprec_wrk(level)%vy2l%zero() + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,& + & p%precv(level)%base_desc,info) + ! all iterations after the first must use BETA = 1 + beta_ = sone + ! + ! Next iteration should use the current residual to compute a correction + ! + call psb_geaxpby(sone,x,szero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%base_desc,info) + call psb_spmm(-sone,p%precv(level)%base_a,y,& + & sone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) + end do + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -380,8 +418,9 @@ subroutine mld_smlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if - call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& + 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) diff --git a/mlprec/impl/mld_sprecset.F90 b/mlprec/impl/mld_sprecset.F90 index d9400333..a8b5c289 100644 --- a/mlprec/impl/mld_sprecset.F90 +++ b/mlprec/impl/mld_sprecset.F90 @@ -142,6 +142,9 @@ subroutine mld_sprecseti(p,what,val,info,ilev,pos) case(mld_max_prec_levs_) p%max_prec_levs = max(val,1) return + case(mld_outer_sweeps_) + p%outer_sweeps = max(val,1) + return end select ! diff --git a/mlprec/impl/mld_zcprecset.F90 b/mlprec/impl/mld_zcprecset.F90 index b2aba3df..338616fb 100644 --- a/mlprec/impl/mld_zcprecset.F90 +++ b/mlprec/impl/mld_zcprecset.F90 @@ -149,6 +149,9 @@ subroutine mld_zcprecseti(p,what,val,info,ilev,pos) case('MAX_PREC_LEVS') p%max_prec_levs = max(val,1) return + case ('OUTER_SWEEPS') + p%outer_sweeps = max(val,1) + return end select diff --git a/mlprec/impl/mld_zmlprec_aply.f90 b/mlprec/impl/mld_zmlprec_aply.f90 index 1ab04619..689fdb8c 100644 --- a/mlprec/impl/mld_zmlprec_aply.f90 +++ b/mlprec/impl/mld_zmlprec_aply.f90 @@ -314,10 +314,11 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) ! Local variables integer(psb_ipk_) :: ictxt, np, me - integer(psb_ipk_) :: debug_level, debug_unit, nlev,nc2l,nr2l,level, err_act + integer(psb_ipk_) :: debug_level, debug_unit + integer(psb_ipk_) :: nlev,nc2l, level, isweep, err_act character(len=20) :: name character :: trans_ - complex(psb_dpk_) :: par + complex(psb_dpk_) :: beta_ type mld_mlprec_wrk_type complex(psb_dpk_), allocatable :: tx(:), ty(:), x2l(:), y2l(:) type(psb_z_vect_type) :: vtx, vty, vx2l, vy2l @@ -338,7 +339,7 @@ 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) + nlev = size(p%precv) allocate(mlprec_wrk(nlev),stat=info) if (info /= psb_success_) then call psb_errpush(psb_err_from_subroutine_,name,a_err='Allocate') @@ -366,12 +367,49 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if end do + ! + ! At first iteration we must use the input BETA + ! + beta_ = beta + level = 1 call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) - call mlprec_wrk(level)%vy2l%zero() + do isweep = 1, p%outer_sweeps - 1 + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! + call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) + + if (info /= psb_success_) then + call psb_errpush(psb_err_internal_error_,name,& + & a_err='Inner prec aply') + goto 9999 + end if + call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta_,y,& + & p%precv(level)%base_desc,info) + ! all iterations after the first must use BETA = 1 + beta_ = zone + ! + ! Next iteration should use the current residual to compute a correction + ! + call psb_geaxpby(zone,x,zzero,mlprec_wrk(level)%vx2l,& + & p%precv(level)%base_desc,info) + call psb_spmm(-zone,p%precv(level)%base_a,y,& + & zone,mlprec_wrk(level)%vx2l,p%precv(level)%base_desc,info) + end do + ! + ! If outer_sweeps == 1 we have just skipped the loop, and it's + ! equivalent to a single application. + ! + + ! + ! With the current implementation, y2l is zeroed internally at first smoother. + ! call mlprec_wrk(level)%vy2l%zero() + ! call inner_ml_aply(level,p,mlprec_wrk,trans_,work,info) if (info /= psb_success_) then @@ -380,8 +418,9 @@ subroutine mld_zmlprec_aply_vect(alpha,p,x,beta,y,desc_data,trans,work,info) goto 9999 end if - call psb_geaxpby(alpha,mlprec_wrk(level)%vy2l,beta,y,& + 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) diff --git a/mlprec/impl/mld_zprecset.F90 b/mlprec/impl/mld_zprecset.F90 index a4cd4fef..e57c4c38 100644 --- a/mlprec/impl/mld_zprecset.F90 +++ b/mlprec/impl/mld_zprecset.F90 @@ -148,6 +148,9 @@ subroutine mld_zprecseti(p,what,val,info,ilev,pos) case(mld_max_prec_levs_) p%max_prec_levs = max(val,1) return + case(mld_outer_sweeps_) + p%outer_sweeps = max(val,1) + return end select ! diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index 48cbe000..d76806e1 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -166,6 +166,7 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_n_prec_levs_ = 38 integer(psb_ipk_), parameter :: mld_max_prec_levs_ = 39 integer(psb_ipk_), parameter :: mld_min_aggr_ratio_ = 40 + integer(psb_ipk_), parameter :: mld_outer_sweeps_ = 41 integer(psb_ipk_), parameter :: mld_ifpsz_ = 42 ! @@ -495,6 +496,8 @@ contains val = mld_filter_mat_ case('NO_FILTER') val = mld_no_filter_mat_ + case('OUTER_SWEEPS') + val = mld_outer_sweeps_ case default val = -1 end select diff --git a/mlprec/mld_c_prec_type.f90 b/mlprec/mld_c_prec_type.f90 index 7bb565c8..3b9fc0f9 100644 --- a/mlprec/mld_c_prec_type.f90 +++ b/mlprec/mld_c_prec_type.f90 @@ -75,7 +75,10 @@ module mld_c_prec_type ! end type mld_Tprec_type ! ! Note that the levels are numbered in increasing order starting from - ! the finest one and the number of levels is given by size(precv(:)). + ! the finest one and the number of levels is given by size(precv(:)), + ! and that is the id of the coarsest level. + ! In the multigrid literature authors often number the levels in decreasing + ! order, with level 0 being the id of the coarsest level. ! ! @@ -84,15 +87,19 @@ module mld_c_prec_type ! ! Aggregation defaults: ! - ! 1. coarse_aggr_size = 0 Default target size will be computed (N_fine)**(1./3.) + ! 1. coarse_aggr_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) integer(psb_ipk_) :: coarse_aggr_size = izero ! 2. n_prec_levs = -1 Use aggregate size to stop integer(psb_ipk_) :: n_prec_levs = -ione - ! 3. Don't use more than 20 levels + ! 3. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_prec_levs = 20_psb_ipk_ ! 4. min_aggr_ratio = 1.5 real(psb_spk_) :: min_aggr_ratio = 1.5_psb_spk_ real(psb_spk_) :: op_complexity=szero + ! + ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. + ! + integer(psb_ipk_) :: outer_sweeps = 1 type(mld_c_onelev_type), allocatable :: precv(:) contains procedure, pass(prec) :: psb_c_apply2_vect => mld_c_apply2_vect @@ -492,6 +499,7 @@ contains ! if (nlev > 1) then write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) 'Outer sweeps:',p%outer_sweeps write(iout_,*) write(iout_,*) 'Base preconditioner (smoother) details' endif @@ -805,6 +813,10 @@ contains class is (mld_cprec_type) pout%ictxt = prec%ictxt pout%coarse_aggr_size = prec%coarse_aggr_size + pout%n_prec_levs = prec%n_prec_levs + pout%max_prec_levs = prec%max_prec_levs + pout%coarse_aggr_size = prec%coarse_aggr_size + pout%outer_sweeps = prec%outer_sweeps pout%op_complexity = prec%op_complexity if (allocated(prec%precv)) then ln = size(prec%precv) diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index f4f594fb..f843ea0b 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -75,7 +75,10 @@ module mld_d_prec_type ! end type mld_Tprec_type ! ! Note that the levels are numbered in increasing order starting from - ! the finest one and the number of levels is given by size(precv(:)). + ! the finest one and the number of levels is given by size(precv(:)), + ! and that is the id of the coarsest level. + ! In the multigrid literature authors often number the levels in decreasing + ! order, with level 0 being the id of the coarsest level. ! ! @@ -84,15 +87,19 @@ module mld_d_prec_type ! ! Aggregation defaults: ! - ! 1. coarse_aggr_size = 0 Default target size will be computed (N_fine)**(1./3.) + ! 1. coarse_aggr_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) integer(psb_ipk_) :: coarse_aggr_size = izero ! 2. n_prec_levs = -1 Use aggregate size to stop integer(psb_ipk_) :: n_prec_levs = -ione - ! 3. Don't use more than 20 levels + ! 3. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_prec_levs = 20_psb_ipk_ ! 4. min_aggr_ratio = 1.5 real(psb_dpk_) :: min_aggr_ratio = 1.5_psb_dpk_ real(psb_dpk_) :: op_complexity=dzero + ! + ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. + ! + integer(psb_ipk_) :: outer_sweeps = 1 type(mld_d_onelev_type), allocatable :: precv(:) contains procedure, pass(prec) :: psb_d_apply2_vect => mld_d_apply2_vect @@ -492,6 +499,7 @@ contains ! if (nlev > 1) then write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) 'Outer sweeps:',p%outer_sweeps write(iout_,*) write(iout_,*) 'Base preconditioner (smoother) details' endif @@ -805,6 +813,10 @@ contains class is (mld_dprec_type) pout%ictxt = prec%ictxt pout%coarse_aggr_size = prec%coarse_aggr_size + pout%n_prec_levs = prec%n_prec_levs + pout%max_prec_levs = prec%max_prec_levs + pout%coarse_aggr_size = prec%coarse_aggr_size + pout%outer_sweeps = prec%outer_sweeps pout%op_complexity = prec%op_complexity if (allocated(prec%precv)) then ln = size(prec%precv) diff --git a/mlprec/mld_s_prec_type.f90 b/mlprec/mld_s_prec_type.f90 index 623e1698..692ed165 100644 --- a/mlprec/mld_s_prec_type.f90 +++ b/mlprec/mld_s_prec_type.f90 @@ -75,7 +75,10 @@ module mld_s_prec_type ! end type mld_Tprec_type ! ! Note that the levels are numbered in increasing order starting from - ! the finest one and the number of levels is given by size(precv(:)). + ! the finest one and the number of levels is given by size(precv(:)), + ! and that is the id of the coarsest level. + ! In the multigrid literature authors often number the levels in decreasing + ! order, with level 0 being the id of the coarsest level. ! ! @@ -84,15 +87,19 @@ module mld_s_prec_type ! ! Aggregation defaults: ! - ! 1. coarse_aggr_size = 0 Default target size will be computed (N_fine)**(1./3.) + ! 1. coarse_aggr_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) integer(psb_ipk_) :: coarse_aggr_size = izero ! 2. n_prec_levs = -1 Use aggregate size to stop integer(psb_ipk_) :: n_prec_levs = -ione - ! 3. Don't use more than 20 levels + ! 3. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_prec_levs = 20_psb_ipk_ ! 4. min_aggr_ratio = 1.5 real(psb_spk_) :: min_aggr_ratio = 1.5_psb_spk_ real(psb_spk_) :: op_complexity=szero + ! + ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. + ! + integer(psb_ipk_) :: outer_sweeps = 1 type(mld_s_onelev_type), allocatable :: precv(:) contains procedure, pass(prec) :: psb_s_apply2_vect => mld_s_apply2_vect @@ -492,6 +499,7 @@ contains ! if (nlev > 1) then write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) 'Outer sweeps:',p%outer_sweeps write(iout_,*) write(iout_,*) 'Base preconditioner (smoother) details' endif @@ -805,6 +813,10 @@ contains class is (mld_sprec_type) pout%ictxt = prec%ictxt pout%coarse_aggr_size = prec%coarse_aggr_size + pout%n_prec_levs = prec%n_prec_levs + pout%max_prec_levs = prec%max_prec_levs + pout%coarse_aggr_size = prec%coarse_aggr_size + pout%outer_sweeps = prec%outer_sweeps pout%op_complexity = prec%op_complexity if (allocated(prec%precv)) then ln = size(prec%precv) diff --git a/mlprec/mld_z_prec_type.f90 b/mlprec/mld_z_prec_type.f90 index c75c6a1a..6afcc2bb 100644 --- a/mlprec/mld_z_prec_type.f90 +++ b/mlprec/mld_z_prec_type.f90 @@ -75,7 +75,10 @@ module mld_z_prec_type ! end type mld_Tprec_type ! ! Note that the levels are numbered in increasing order starting from - ! the finest one and the number of levels is given by size(precv(:)). + ! the finest one and the number of levels is given by size(precv(:)), + ! and that is the id of the coarsest level. + ! In the multigrid literature authors often number the levels in decreasing + ! order, with level 0 being the id of the coarsest level. ! ! @@ -84,15 +87,19 @@ module mld_z_prec_type ! ! Aggregation defaults: ! - ! 1. coarse_aggr_size = 0 Default target size will be computed (N_fine)**(1./3.) + ! 1. coarse_aggr_size = 0 Default target size will be computed as 40*(N_fine)**(1./3.) integer(psb_ipk_) :: coarse_aggr_size = izero ! 2. n_prec_levs = -1 Use aggregate size to stop integer(psb_ipk_) :: n_prec_levs = -ione - ! 3. Don't use more than 20 levels + ! 3. maximum number of levels. Defaults to 20 integer(psb_ipk_) :: max_prec_levs = 20_psb_ipk_ ! 4. min_aggr_ratio = 1.5 real(psb_dpk_) :: min_aggr_ratio = 1.5_psb_dpk_ real(psb_dpk_) :: op_complexity=dzero + ! + ! Number of outer sweeps. Sometimes 2 V-cycles may be better than 1 W-cycle. + ! + integer(psb_ipk_) :: outer_sweeps = 1 type(mld_z_onelev_type), allocatable :: precv(:) contains procedure, pass(prec) :: psb_z_apply2_vect => mld_z_apply2_vect @@ -492,6 +499,7 @@ contains ! if (nlev > 1) then write(iout_,*) 'Multilevel Schwarz' + write(iout_,*) 'Outer sweeps:',p%outer_sweeps write(iout_,*) write(iout_,*) 'Base preconditioner (smoother) details' endif @@ -805,6 +813,10 @@ contains class is (mld_zprec_type) pout%ictxt = prec%ictxt pout%coarse_aggr_size = prec%coarse_aggr_size + pout%n_prec_levs = prec%n_prec_levs + pout%max_prec_levs = prec%max_prec_levs + pout%coarse_aggr_size = prec%coarse_aggr_size + pout%outer_sweeps = prec%outer_sweeps pout%op_complexity = prec%op_complexity if (allocated(prec%precv)) then ln = size(prec%precv)