diff --git a/mlprec/mld_base_prec_type.f90 b/mlprec/mld_base_prec_type.f90 index a99f3e40..643d02ca 100644 --- a/mlprec/mld_base_prec_type.f90 +++ b/mlprec/mld_base_prec_type.f90 @@ -1,6 +1,6 @@ !!$ !!$ -!!$ MLD2P4 version 1.1 +!!$ MLD2P4 version 1.2 !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package !!$ based on PSBLAS (Parallel Sparse BLAS version 2.3.1) !!$ @@ -36,13 +36,11 @@ !!$ POSSIBILITY OF SUCH DAMAGE. !!$ !!$ -! File: mld_prec_type.f90 +! File: mld_base_prec_type.f90 ! -! Module: mld_prec_type +! Module: mld_base_prec_type ! -! This module defines: -! - the mld_prec_type data structure containing the preconditioner and related -! data structures; +! Constants and utilities in common to all type variants of MLD preconditioners. ! - integer constants defining the preconditioner; ! - character constants describing the preconditioner (used by the routines ! printing out a preconditioner description); @@ -99,7 +97,6 @@ module mld_base_prec_type integer, parameter :: mld_sub_ren_ = 5 integer, parameter :: mld_sub_ovr_ = 6 integer, parameter :: mld_sub_fillin_ = 8 - integer, parameter :: mld_smoother_sweeps_ = 9 !! 2 ints for 64 bit versions integer, parameter :: mld_slu_ptr_ = 10 integer, parameter :: mld_umf_symptr_ = 12 @@ -109,25 +106,28 @@ module mld_base_prec_type ! ! These are in onelev ! - integer, parameter :: mld_ml_type_ = 20 - integer, parameter :: mld_smoother_pos_ = 21 - integer, parameter :: mld_aggr_kind_ = 22 - integer, parameter :: mld_aggr_alg_ = 23 - integer, parameter :: mld_aggr_omega_alg_ = 24 - integer, parameter :: mld_aggr_eig_ = 25 - integer, parameter :: mld_aggr_filter_ = 26 - integer, parameter :: mld_coarse_mat_ = 27 - integer, parameter :: mld_coarse_solve_ = 28 - integer, parameter :: mld_coarse_sweeps_ = 29 - integer, parameter :: mld_coarse_fillin_ = 30 - integer, parameter :: mld_coarse_subsolve_ = 31 - integer, parameter :: mld_ifpsz_ = 32 + integer, parameter :: mld_ml_type_ = 20 + integer, parameter :: mld_smoother_sweeps_pre_ = 21 + integer, parameter :: mld_smoother_sweeps_post_ = 22 + integer, parameter :: mld_smoother_pos_ = 23 + integer, parameter :: mld_aggr_kind_ = 24 + integer, parameter :: mld_aggr_alg_ = 25 + integer, parameter :: mld_aggr_omega_alg_ = 26 + integer, parameter :: mld_aggr_eig_ = 27 + integer, parameter :: mld_aggr_filter_ = 28 + integer, parameter :: mld_coarse_mat_ = 29 + integer, parameter :: mld_coarse_solve_ = 30 + integer, parameter :: mld_coarse_sweeps_ = 31 + integer, parameter :: mld_coarse_fillin_ = 32 + integer, parameter :: mld_coarse_subsolve_ = 33 + integer, parameter :: mld_smoother_sweeps_ = 34 + integer, parameter :: mld_ifpsz_ = 36 ! ! Legal values for entry: mld_smoother_type_ ! - integer, parameter :: mld_min_prec_=0, mld_noprec_=0, mld_diag_=1, mld_bjac_=2,& - & mld_as_=3, mld_pjac_=4, mld_max_prec_=4 + integer, parameter :: mld_min_prec_=0, mld_noprec_=0, mld_jac_=1, mld_bjac_=2,& + & mld_as_=3, mld_max_prec_=3 ! VERY IMPORTANT: we are relying on the following to be true: ! mld_pjac_ == mld_diag_scale_ ! mld_bjac_ == mld_milu_n_ (or mld_ilu_n_ would be fine) @@ -136,8 +136,8 @@ module mld_base_prec_type ! ! Legal values for entry: mld_sub_solve_ ! - integer, parameter :: mld_f_none_=0,mld_ilu_n_=1,mld_milu_n_=2, mld_ilu_t_=3 - integer, parameter :: mld_diag_scale_=4, mld_slu_=5, mld_umf_=6, mld_sludist_=7 + integer, parameter :: mld_f_none_=0, mld_diag_scale_=1, mld_ilu_n_=2, mld_milu_n_=3 + integer, parameter :: mld_ilu_t_=4, mld_slu_=5, mld_umf_=6, mld_sludist_=7 integer, parameter :: mld_max_sub_solve_= 7 ! ! Legal values for entry: mld_sub_ren_ @@ -215,7 +215,7 @@ module mld_base_prec_type character(len=19), parameter, private :: & & eigen_estimates(0:0)=(/'infinity norm '/) character(len=19), parameter, private :: & - & smooth_names(1:3)=(/'pre-smoothing ','post-smoothing ',& + & smooth_pos_names(1:3)=(/'pre-smoothing ','post-smoothing ',& & 'pre/post-smoothing'/) character(len=15), parameter, private :: & & aggr_kinds(0:3)=(/'nonsmoothed ','smoothed ',& @@ -233,10 +233,10 @@ module mld_base_prec_type & ml_names(0:3)=(/'none ','additive ','multiplicative',& & 'new ML '/) character(len=15), parameter, private :: & - & fact_names(0:7)=(/'none ','ILU(n) ',& + & fact_names(0:7)=(/'none ','Point Jacobi ','ILU(n) ',& & 'MILU(n) ','ILU(t,n) ',& - & 'SuperLU ','UMFPACK LU ',& - & 'SuperLU_Dist ','DiagSc-PntJac '/) + & 'UMFPACK LU ',& + & 'SuperLU_Dist ','SuperLU '/) interface mld_check_def module procedure mld_icheck_def, mld_scheck_def, mld_dcheck_def @@ -319,12 +319,12 @@ contains val = mld_twoside_smooth_ case('NOPREC') val = mld_noprec_ - case('DIAG') - val = mld_diag_ +!!$ case('DIAG') +!!$ val = mld_diag_ case('BJAC') val = mld_bjac_ - case('PJAC', 'JACOBI') - val = mld_pjac_ + case('JAC','JACOBI') + val = mld_jac_ case('AS') val = mld_as_ case('RENUM_NONE') @@ -374,10 +374,8 @@ contains select case(iprcparm(mld_smoother_type_)) case(mld_noprec_) write(iout,*) ' No preconditioning' - case(mld_diag_) - write(iout,*) ' Diagonal scaling' - case(mld_pjac_) - write(iout,*) ' Point Jacobi ' + case(mld_jac_) + write(iout,*) ' Jacobi ' case(mld_bjac_) write(iout,*) ' Block Jacobi with ',& & fact_names(iprcparm(mld_sub_solve_)) @@ -428,6 +426,7 @@ contains integer, intent(out) :: info real(psb_spk_), intent(in), optional :: rprcparm(:) real(psb_dpk_), intent(in), optional :: dprcparm(:) + integer :: sweeps info = 0 if (count((/ present(rprcparm),present(dprcparm) /)) /= 1) then @@ -441,7 +440,26 @@ contains write(iout,*) ' Multilevel type: ',& & ml_names(iprcparm(mld_ml_type_)) write(iout,*) ' Smoother position: ',& - & smooth_names(iprcparm(mld_smoother_pos_)) + & smooth_pos_names(iprcparm(mld_smoother_pos_)) + if (iprcparm(mld_ml_type_) == mld_add_ml_) then + write(iout,*) ' Number of sweeps : ',& + & iprcparm(mld_smoother_sweeps_) + else + select case (iprcparm(mld_smoother_pos_)) + case (mld_pre_smooth_) + write(iout,*) ' Number of sweeps : ',& + & iprcparm(mld_smoother_sweeps_pre_) + case (mld_post_smooth_) + write(iout,*) ' Number of sweeps : ',& + & iprcparm(mld_smoother_sweeps_post_) + case (mld_twoside_smooth_) + write(iout,*) ' Number of sweeps : pre: ',& + & iprcparm(mld_smoother_sweeps_pre_) ,& + & ' post: ',& + & iprcparm(mld_smoother_sweeps_post_) + end select + end if + write(iout,*) ' Aggregation: ', & & aggr_names(iprcparm(mld_aggr_alg_)) write(iout,*) ' Aggregation type: ', & @@ -550,10 +568,18 @@ contains end if if (iprcparm(mld_coarse_mat_) == mld_distr_mat_ .and. & & iprcparm(mld_sub_solve_) /= mld_sludist_) then - write(iout,*) ' Coarsest matrix solver: block Jacobi with ', & - & fact_names(iprcparm2(mld_sub_solve_)) - write(iout,*) ' Number of Jacobi sweeps: ', & - & (iprcparm2(mld_smoother_sweeps_)) +!!$ write(iout,*) ' Coarsest matrix solver: ',& +!!$ & smoother_names(iprcparm2(mld_smoother_type_)) + select case (iprcparm2(mld_smoother_type_)) + case(mld_bjac_,mld_as_) + write(iout,*) ' subdomain solver: ',& + & fact_names(iprcparm2(mld_sub_solve_)) + write(iout,*) ' Number of smoother sweeps: ', & + & (iprcparm2(mld_smoother_sweeps_)) + case(mld_jac_) + write(iout,*) ' Number of smoother sweeps: ', & + & (iprcparm2(mld_smoother_sweeps_)) + end select else write(iout,*) ' Coarsest matrix solver: ', & & fact_names(iprcparm2(mld_sub_solve_)) @@ -837,8 +863,8 @@ contains select case(iprec) case(mld_noprec_) pr_to_str='NOPREC' - case(mld_diag_) - pr_to_str='DIAG' + case(mld_jac_) + pr_to_str='JAC' case(mld_bjac_) pr_to_str='BJAC' case(mld_as_) diff --git a/mlprec/mld_d_as_smoother.f03 b/mlprec/mld_d_as_smoother.f03 index 55eca1dc..1a943c5f 100644 --- a/mlprec/mld_d_as_smoother.f03 +++ b/mlprec/mld_d_as_smoother.f03 @@ -53,7 +53,7 @@ module mld_d_as_smoother ! type(psb_d_sparse_mat) :: nd type(psb_desc_type) :: desc_data - integer :: novr, sweeps, restr, prol + integer :: novr, restr, prol contains procedure, pass(sm) :: build => d_as_smoother_bld procedure, pass(sm) :: apply => d_as_smoother_apply @@ -79,7 +79,7 @@ module mld_d_as_smoother contains - subroutine d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,work,info) + subroutine d_as_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) use psb_sparse_mod type(psb_desc_type), intent(in) :: desc_data class(mld_d_as_smoother_type), intent(in) :: sm @@ -87,6 +87,7 @@ contains real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) integer, intent(out) :: info @@ -154,7 +155,7 @@ contains endif - if ((sm%novr == 0).and.(sm%sweeps == 1)) then + if ((sm%novr == 0).and.(sweeps == 1)) then ! ! Shortcut: in this case it's just the same ! as Block Jacobi. @@ -173,7 +174,7 @@ contains tx(nrow_d+1:isz) = dzero - if (sm%sweeps == 1) then + if (sweeps == 1) then select case(trans_) case('N') @@ -311,7 +312,7 @@ contains - else if (sm%sweeps > 1) then + else if (sweeps > 1) then ! ! @@ -320,7 +321,7 @@ contains ! ! ty = dzero - do i=1, sm%sweeps + do i=1, sweeps select case(trans_) case('N') ! @@ -473,7 +474,7 @@ contains info = 10 call psb_errpush(info,name,& - & i_err=(/2,sm%sweeps,0,0,0/)) + & i_err=(/2,sweeps,0,0,0/)) goto 9999 @@ -664,8 +665,8 @@ contains call psb_erractionsave(err_act) select case(what) - case(mld_smoother_sweeps_) - sm%sweeps = val +!!$ case(mld_smoother_sweeps_) +!!$ sm%sweeps = val case(mld_sub_ovr_) sm%novr = val case(mld_sub_restr_) @@ -834,7 +835,7 @@ contains endif write(iout_,*) ' Additive Schwarz with ',& - & sm%sweeps,' sweeps and ',sm%novr, ' overlap layers.' + & sm%novr, ' overlap layers.' write(iout_,*) ' Restrictor: ',restrict_names(sm%restr) write(iout_,*) ' Prolongator: ',prolong_names(sm%prol) write(iout_,*) ' Local solver:' diff --git a/mlprec/mld_d_jac_smoother.f03 b/mlprec/mld_d_jac_smoother.f03 index 67ef05be..e2fb92dc 100644 --- a/mlprec/mld_d_jac_smoother.f03 +++ b/mlprec/mld_d_jac_smoother.f03 @@ -52,7 +52,6 @@ module mld_d_jac_smoother ! class(mld_d_base_solver_type), allocatable :: sv ! type(psb_d_sparse_mat) :: nd - integer :: sweeps contains procedure, pass(sm) :: build => d_jac_smoother_bld procedure, pass(sm) :: apply => d_jac_smoother_apply @@ -74,7 +73,7 @@ module mld_d_jac_smoother contains - subroutine d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,work,info) + subroutine d_jac_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) use psb_sparse_mod type(psb_desc_type), intent(in) :: desc_data class(mld_d_jac_smoother_type), intent(in) :: sm @@ -82,6 +81,7 @@ contains real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) integer, intent(out) :: info @@ -136,7 +136,7 @@ contains end if endif - if (sm%sweeps == 1) then + if (sweeps == 1) then call sm%sv%apply(alpha,x,beta,y,desc_data,trans_,aux,info) @@ -145,11 +145,11 @@ contains goto 9999 endif - else if (sm%sweeps > 1) then + else if (sweeps > 1) then ! ! - ! Apply prec%iprcparm(mld_smoother_sweeps_) sweeps of a block-Jacobi solver + ! Apply multiple sweeps of a block-Jacobi solver ! to compute an approximate solution of a linear system. ! ! @@ -163,7 +163,7 @@ contains tx = dzero ty = dzero - do i=1, sm%sweeps + do i=1, sweeps ! ! Compute Y(j+1) = D^(-1)*(X-ND*Y(j)), where D and ND are the ! block diagonal part and the remaining part of the local matrix @@ -198,7 +198,7 @@ contains info = 10 call psb_errpush(info,name,& - & i_err=(/2,sm%sweeps,0,0,0/)) + & i_err=(/2,sweeps,0,0,0/)) goto 9999 endif @@ -309,8 +309,8 @@ contains call psb_erractionsave(err_act) select case(what) - case(mld_smoother_sweeps_) - sm%sweeps = val +!!$ case(mld_smoother_sweeps_) +!!$ sm%sweeps = val case default if (allocated(sm%sv)) then call sm%sv%set(what,val,info) @@ -472,8 +472,7 @@ contains iout_ = 6 endif - write(iout_,*) ' Block Jacobi smoother with ',& - & sm%sweeps,' sweeps.' + write(iout_,*) ' Block Jacobi smoother ' write(iout_,*) ' Local solver:' if (allocated(sm%sv)) then call sm%sv%descr(info,iout_) diff --git a/mlprec/mld_d_prec_type.f03 b/mlprec/mld_d_prec_type.f03 index dc1e5bed..994dcddf 100644 --- a/mlprec/mld_d_prec_type.f03 +++ b/mlprec/mld_d_prec_type.f03 @@ -672,7 +672,7 @@ contains end subroutine mld_dprec_free - subroutine d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,work,info) + subroutine d_base_smoother_apply(alpha,sm,x,beta,y,desc_data,trans,sweeps,work,info) use psb_sparse_mod type(psb_desc_type), intent(in) :: desc_data class(mld_d_base_smoother_type), intent(in) :: sm @@ -680,6 +680,7 @@ contains real(psb_dpk_),intent(inout) :: y(:) real(psb_dpk_),intent(in) :: alpha,beta character(len=1),intent(in) :: trans + integer, intent(in) :: sweeps real(psb_dpk_),target, intent(inout) :: work(:) integer, intent(out) :: info diff --git a/mlprec/mld_dmlprec_aply.f90 b/mlprec/mld_dmlprec_aply.f90 index 773ffbee..9c28d754 100644 --- a/mlprec/mld_dmlprec_aply.f90 +++ b/mlprec/mld_dmlprec_aply.f90 @@ -420,7 +420,8 @@ contains ! Local variables integer :: ictxt,np,me,i, nr2l,nc2l,err_act integer :: debug_level, debug_unit - integer :: nlev, ilev + integer :: nlev, ilev, sweeps + character(len=20) :: name name = 'inner_ml_aply' @@ -479,13 +480,14 @@ contains end if - + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) !!$ call mld_baseprec_aply(done,p%precv(level)%prec,& !!$ & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& !!$ & p%precv(level)%base_desc, trans,work,info) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) if (info /= 0) goto 9999 @@ -550,14 +552,17 @@ contains & work=work,trans=trans) if (info /= 0) goto 9999 + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) else - + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) end if @@ -583,10 +588,15 @@ contains ! ! Apply the base preconditioner ! + if (level < nlev) then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) - + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) if (info /= 0) goto 9999 @@ -637,10 +647,15 @@ contains ! ! Apply the base preconditioner ! + if (level < nlev) then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) - + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) if (info /= 0) goto 9999 @@ -699,15 +714,17 @@ contains & work=work,trans=trans) if (info /= 0) goto 9999 + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) else - + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) - + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) end if case default @@ -744,10 +761,19 @@ contains ! ! Apply the base preconditioner ! + if (level < nlev) then + if (trans == 'N') then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + end if + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_) + end if if (info == 0) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%x2l,dzero,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) - + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) ! ! Compute the residual (at all levels but the coarsest one) ! and call recursively @@ -782,11 +808,15 @@ contains ! ! Apply the base preconditioner ! + if (trans == 'N') then + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_post_) + else + sweeps = p%precv(level)%iprcparm(mld_smoother_sweeps_pre_) + end if if (info == 0) call p%precv(level)%sm%apply(done,& & mlprec_wrk(level)%tx,done,mlprec_wrk(level)%y2l,& - & p%precv(level)%base_desc, trans,work,info) - - + & p%precv(level)%base_desc, trans,& + & sweeps,work,info) if (info /= 0) then call psb_errpush(4001,name,a_err='Error: residual/baseprec_aply') goto 9999 @@ -822,7 +852,6 @@ contains return end subroutine inner_ml_aply - end subroutine mld_dmlprec_aply diff --git a/mlprec/mld_dmlprec_bld.f90 b/mlprec/mld_dmlprec_bld.f90 index 06c3c01c..674d555a 100644 --- a/mlprec/mld_dmlprec_bld.f90 +++ b/mlprec/mld_dmlprec_bld.f90 @@ -307,7 +307,11 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) call mld_check_def(p%precv(i)%prec%rprcparm(mld_sub_iluthrs_),& & 'Eps',dzero,is_legal_fact_thrs) end select - call mld_check_def(p%precv(i)%prec%iprcparm(mld_smoother_sweeps_),& + call mld_check_def(p%precv(i)%iprcparm(mld_smoother_sweeps_),& + & 'Jacobi sweeps',1,is_legal_jac_sweeps) + call mld_check_def(p%precv(i)%iprcparm(mld_smoother_sweeps_pre_),& + & 'Jacobi sweeps',1,is_legal_jac_sweeps) + call mld_check_def(p%precv(i)%iprcparm(mld_smoother_sweeps_post_),& & 'Jacobi sweeps',1,is_legal_jac_sweeps) ! @@ -322,7 +326,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) endif end if select case (p%precv(i)%prec%iprcparm(mld_smoother_type_)) - case(mld_diag_, mld_bjac_, mld_pjac_) + case(mld_bjac_, mld_jac_) allocate(mld_d_jac_smoother_type :: p%precv(i)%sm, stat=info) case(mld_as_) allocate(mld_d_as_smoother_type :: p%precv(i)%sm, stat=info) @@ -338,8 +342,6 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) call p%precv(i)%sm%set(mld_sub_restr_,p%precv(i)%prec%iprcparm(mld_sub_restr_),info) call p%precv(i)%sm%set(mld_sub_prol_,p%precv(i)%prec%iprcparm(mld_sub_prol_),info) call p%precv(i)%sm%set(mld_sub_ovr_,p%precv(i)%prec%iprcparm(mld_sub_ovr_),info) - call p%precv(i)%sm%set(mld_smoother_sweeps_,& - & p%precv(i)%prec%iprcparm(mld_smoother_sweeps_),info) select case (p%precv(i)%prec%iprcparm(mld_sub_solve_)) case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) @@ -416,6 +418,35 @@ contains ! val = prec%iprcparm(mld_coarse_solve_) select case (val) + case(mld_jac_) + + if (prec%prec%iprcparm(mld_sub_solve_) /= mld_diag_scale_) then + if (me == 0) write(debug_unit,*)& + & 'Warning: inconsistent coarse level specification.' + if (me == 0) write(debug_unit,*)& + & ' Resetting according to the value specified for mld_coarse_solve_.' + prec%prec%iprcparm(mld_sub_solve_) = mld_diag_scale_ + end if + prec%prec%iprcparm(mld_smoother_type_) = mld_jac_ + + case(mld_bjac_) + + if ((prec%prec%iprcparm(mld_sub_solve_) == mld_diag_scale_).or.& + & ( prec%prec%iprcparm(mld_smoother_type_) /= mld_bjac_)) then + if (me == 0) write(debug_unit,*)& + & 'Warning: inconsistent coarse level specification.' + if (me == 0) write(debug_unit,*)& + & ' Resetting according to the value specified for mld_coarse_solve_.' +!!$#if defined(HAVE_UMF_) +!!$ prec%prec%iprcparm(mld_sub_solve_) = mld_umf_ +!!$#elif defined(HAVE_SLU_) +!!$ prec%prec%iprcparm(mld_sub_solve_) = mld_slu_ +!!$#else + prec%prec%iprcparm(mld_sub_solve_) = mld_ilu_n_ +!!$#endif + end if + prec%prec%iprcparm(mld_smoother_type_) = mld_bjac_ + case(mld_umf_, mld_slu_) if ((prec%iprcparm(mld_coarse_mat_) /= mld_repl_mat_).or.& & (prec%prec%iprcparm(mld_sub_solve_) /= val)) then diff --git a/mlprec/mld_dprecaply.f90 b/mlprec/mld_dprecaply.f90 index e05941b7..43f04a29 100644 --- a/mlprec/mld_dprecaply.f90 +++ b/mlprec/mld_dprecaply.f90 @@ -141,7 +141,8 @@ subroutine mld_dprecaply(prec,x,y,desc_data,info,trans,work) ! Number of levels = 1: apply the base preconditioner ! !!$ call mld_baseprec_aply(done,prec%precv(1)%prec,x,dzero,y,desc_data,trans_, work_,info) - call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_, work_,info) + call prec%precv(1)%sm%apply(done,x,dzero,y,desc_data,trans_,& + & prec%precv(1)%iprcparm(mld_smoother_sweeps_), work_,info) else info = 4013 call psb_errpush(info,name,a_err='Invalid size of precv',& diff --git a/mlprec/mld_dprecbld.f90 b/mlprec/mld_dprecbld.f90 index bc39e2ad..e17b3491 100644 --- a/mlprec/mld_dprecbld.f90 +++ b/mlprec/mld_dprecbld.f90 @@ -184,7 +184,7 @@ subroutine mld_dprecbld(a,desc_a,p,info) call mld_check_def(p%precv(1)%prec%rprcparm(mld_sub_iluthrs_),& & 'Eps',dzero,is_legal_fact_thrs) end select - call mld_check_def(p%precv(1)%prec%iprcparm(mld_smoother_sweeps_),& + call mld_check_def(p%precv(1)%iprcparm(mld_smoother_sweeps_),& & 'Jacobi sweeps',1,is_legal_jac_sweeps) !!$ !!$ call mld_baseprec_bld(p%precv(1)%base_a,p%precv(1)%base_desc,& @@ -201,7 +201,7 @@ subroutine mld_dprecbld(a,desc_a,p,info) endif end if select case (p%precv(1)%prec%iprcparm(mld_smoother_type_)) - case(mld_diag_, mld_bjac_, mld_pjac_) + case(mld_jac_, mld_bjac_) allocate(mld_d_jac_smoother_type :: p%precv(1)%sm, stat=info) case(mld_as_) allocate(mld_d_as_smoother_type :: p%precv(1)%sm, stat=info) @@ -217,8 +217,6 @@ subroutine mld_dprecbld(a,desc_a,p,info) call p%precv(1)%sm%set(mld_sub_restr_,p%precv(1)%prec%iprcparm(mld_sub_restr_),info) call p%precv(1)%sm%set(mld_sub_prol_,p%precv(1)%prec%iprcparm(mld_sub_prol_),info) call p%precv(1)%sm%set(mld_sub_ovr_,p%precv(1)%prec%iprcparm(mld_sub_ovr_),info) - call p%precv(1)%sm%set(mld_smoother_sweeps_,& - & p%precv(1)%prec%iprcparm(mld_smoother_sweeps_),info) select case (p%precv(1)%prec%iprcparm(mld_sub_solve_)) case(mld_ilu_n_,mld_milu_n_,mld_ilu_t_) diff --git a/mlprec/mld_dprecinit.F90 b/mlprec/mld_dprecinit.F90 index 941dd8d2..2540ad39 100644 --- a/mlprec/mld_dprecinit.F90 +++ b/mlprec/mld_dprecinit.F90 @@ -51,6 +51,8 @@ ! ! 'DIAG' - diagonal preconditioner ! +! 'PJAC' - point Jacobi preconditioner +! ! 'BJAC' - block Jacobi preconditioner, with ILU(0) ! on the local blocks ! @@ -135,8 +137,11 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_ren_) = 0 p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 - case ('DIAG') + case ('JAC','DIAG','JACOBI') nlev_ = 1 ilev_ = 1 allocate(p%precv(nlev_),stat=info) @@ -149,13 +154,17 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%rprcparm(:) = dzero p%precv(ilev_)%prec%iprcparm(:) = 0 p%precv(ilev_)%prec%rprcparm(:) = dzero - p%precv(ilev_)%prec%iprcparm(mld_smoother_type_) = mld_diag_ - p%precv(ilev_)%prec%iprcparm(mld_sub_solve_) = mld_f_none_ + p%precv(ilev_)%prec%iprcparm(mld_smoother_type_) = mld_jac_ + p%precv(ilev_)%prec%iprcparm(mld_sub_solve_) = mld_diag_scale_ p%precv(ilev_)%prec%iprcparm(mld_sub_restr_) = psb_none_ p%precv(ilev_)%prec%iprcparm(mld_sub_prol_) = psb_none_ - p%precv(ilev_)%prec%iprcparm(mld_sub_ren_) = 0 + p%precv(ilev_)%prec%iprcparm(mld_sub_ren_) = 0 p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 0 + p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 case ('BJAC') nlev_ = 1 @@ -178,6 +187,9 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 0 p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 case ('AS') nlev_ = 1 @@ -200,6 +212,9 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 1 p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 case ('ML') @@ -224,7 +239,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_ p%precv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_ p%precv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_ - p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_post_smooth_ + p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_twoside_smooth_ p%precv(ilev_)%iprcparm(mld_aggr_omega_alg_) = mld_eig_est_ p%precv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_ p%precv(ilev_)%iprcparm(mld_aggr_filter_) = mld_no_filter_mat_ @@ -236,6 +251,9 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 1 p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 if (nlev_ == 1) return do ilev_ = 2, nlev_ -1 @@ -256,11 +274,14 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_solve_) = mld_ilu_n_ p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 1 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 1 p%precv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_ p%precv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_ p%precv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_ p%precv(ilev_)%iprcparm(mld_coarse_mat_) = mld_distr_mat_ - p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_post_smooth_ + p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_twoside_smooth_ p%precv(ilev_)%iprcparm(mld_aggr_omega_alg_) = mld_eig_est_ p%precv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_ p%precv(ilev_)%iprcparm(mld_aggr_filter_) = mld_no_filter_mat_ @@ -292,16 +313,18 @@ subroutine mld_dprecinit(p,ptype,info,nlev) p%precv(ilev_)%prec%iprcparm(mld_sub_ren_) = 0 p%precv(ilev_)%prec%iprcparm(mld_sub_ovr_) = 0 p%precv(ilev_)%prec%iprcparm(mld_sub_fillin_) = 0 - p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = 4 - p%precv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_ - p%precv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_ - p%precv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_ - p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_post_smooth_ - p%precv(ilev_)%iprcparm(mld_aggr_omega_alg_) = mld_eig_est_ - p%precv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_ - p%precv(ilev_)%iprcparm(mld_aggr_filter_) = mld_no_filter_mat_ - p%precv(ilev_)%rprcparm(mld_aggr_omega_val_) = dzero - p%precv(ilev_)%rprcparm(mld_aggr_thresh_) = dzero + p%precv(ilev_)%iprcparm(mld_ml_type_) = mld_mult_ml_ + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = 4 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = 4 + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = 4 + p%precv(ilev_)%iprcparm(mld_aggr_alg_) = mld_dec_aggr_ + p%precv(ilev_)%iprcparm(mld_aggr_kind_) = mld_smooth_prol_ + p%precv(ilev_)%iprcparm(mld_smoother_pos_) = mld_twoside_smooth_ + p%precv(ilev_)%iprcparm(mld_aggr_omega_alg_) = mld_eig_est_ + p%precv(ilev_)%iprcparm(mld_aggr_eig_) = mld_max_norm_ + p%precv(ilev_)%iprcparm(mld_aggr_filter_) = mld_no_filter_mat_ + p%precv(ilev_)%rprcparm(mld_aggr_omega_val_) = dzero + p%precv(ilev_)%rprcparm(mld_aggr_thresh_) = dzero case default write(0,*) name,': Warning: Unknown preconditioner type request "',ptype,'"' diff --git a/mlprec/mld_dprecset.f90 b/mlprec/mld_dprecset.f90 index 046f1cc4..912e8517 100644 --- a/mlprec/mld_dprecset.f90 +++ b/mlprec/mld_dprecset.f90 @@ -143,11 +143,17 @@ subroutine mld_dprecseti(p,what,val,info,ilev) case(mld_smoother_type_) p%precv(ilev_)%iprcparm(what) = val p%precv(ilev_)%prec%iprcparm(what) = val + case(mld_smoother_sweeps_) + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = val + p%precv(ilev_)%prec%iprcparm(what) = val case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& - & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,mld_smoother_sweeps_) + & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) p%precv(ilev_)%prec%iprcparm(what) = val case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_) + & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_) p%precv(ilev_)%iprcparm(what) = val case default write(0,*) name,': Error: invalid WHAT' @@ -159,12 +165,17 @@ subroutine mld_dprecseti(p,what,val,info,ilev) case(mld_smoother_type_) p%precv(ilev_)%iprcparm(what) = val p%precv(ilev_)%prec%iprcparm(what) = val + case(mld_smoother_sweeps_) + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = val + p%precv(ilev_)%prec%iprcparm(what) = val case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& - & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& - & mld_smoother_sweeps_) + & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) p%precv(ilev_)%prec%iprcparm(what) = val case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_) + & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_) p%precv(ilev_)%iprcparm(what) = val case(mld_coarse_mat_) if (ilev_ /= nlev_ .and. val /= mld_distr_mat_) then @@ -206,6 +217,9 @@ subroutine mld_dprecseti(p,what,val,info,ilev) return end if p%precv(ilev_)%prec%iprcparm(mld_smoother_sweeps_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = val case(mld_coarse_fillin_) if (ilev_ /= nlev_) then write(0,*) name,': Error: Inconsistent specification of WHAT vs. ILEV' @@ -227,8 +241,7 @@ subroutine mld_dprecseti(p,what,val,info,ilev) ! select case(what) case(mld_sub_solve_,mld_sub_restr_,mld_sub_prol_,& - & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& - & mld_smoother_sweeps_) + & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) do ilev_=1,max(1,nlev_-1) if (.not.allocated(p%precv(ilev_)%iprcparm)) then write(0,*) name,& @@ -239,6 +252,14 @@ subroutine mld_dprecseti(p,what,val,info,ilev) p%precv(ilev_)%iprcparm(what) = val p%precv(ilev_)%prec%iprcparm(what) = val end do + case(mld_smoother_sweeps_) + do ilev_=1,max(1,nlev_-1) + p%precv(ilev_)%iprcparm(what) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(ilev_)%iprcparm(mld_smoother_sweeps_post_) = val + p%precv(ilev_)%prec%iprcparm(what) = val + end do + case(mld_smoother_type_) do ilev_=1,nlev_ if (.not.allocated(p%precv(ilev_)%iprcparm)) then @@ -251,12 +272,13 @@ subroutine mld_dprecseti(p,what,val,info,ilev) p%precv(ilev_)%prec%iprcparm(what) = val end do case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,mld_aggr_filter_) + & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& + & mld_smoother_pos_,mld_aggr_omega_alg_,& + & mld_aggr_eig_,mld_aggr_filter_) do ilev_=1,nlev_ if (.not.allocated(p%precv(ilev_)%iprcparm)) then write(0,*) name,& - &': Error: uninitialized preconditioner component,',& - &' should call MLD_PRECINIT' + &': Error: uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 return endif @@ -291,9 +313,9 @@ subroutine mld_dprecseti(p,what,val,info,ilev) p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = val case(mld_sludist_) p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = val -!!$ case(mld_jac_) -!!$ p%precv(nlev_)%prec%iprcparm(mld_smoother_type_) = mld_jac_ -!!$ p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = mld_diag_scale_ + case(mld_jac_) + p%precv(nlev_)%prec%iprcparm(mld_smoother_type_) = mld_jac_ + p%precv(nlev_)%prec%iprcparm(mld_sub_solve_) = mld_diag_scale_ end select endif case(mld_coarse_subsolve_) @@ -314,7 +336,13 @@ subroutine mld_dprecseti(p,what,val,info,ilev) info = -1 return endif - if (nlev_ > 1) p%precv(nlev_)%prec%iprcparm(mld_smoother_sweeps_) = val + if (nlev_ > 1) then + p%precv(nlev_)%prec%iprcparm(mld_smoother_sweeps_) = val + p%precv(nlev_)%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(nlev_)%iprcparm(mld_smoother_sweeps_post_) = val + p%precv(nlev_)%prec%iprcparm(mld_smoother_sweeps_pre_) = val + p%precv(nlev_)%prec%iprcparm(mld_smoother_sweeps_post_) = val + end if case(mld_coarse_fillin_) if (.not.allocated(p%precv(nlev_)%iprcparm)) then write(0,*) name,&