From 875443efe72e484b1ccd5cfd535a9bef5e8d2967 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 9 May 2016 16:00:12 +0000 Subject: [PATCH 1/2] mld2p4-2: docs/src/userinterface.tex mlprec/impl/level/mld_c_base_onelev_cseti.f90 mlprec/impl/level/mld_c_base_onelev_seti.f90 mlprec/impl/level/mld_d_base_onelev_cseti.f90 mlprec/impl/level/mld_d_base_onelev_seti.f90 mlprec/impl/level/mld_s_base_onelev_cseti.f90 mlprec/impl/level/mld_s_base_onelev_seti.f90 mlprec/impl/level/mld_z_base_onelev_cseti.f90 mlprec/impl/level/mld_z_base_onelev_seti.f90 mlprec/impl/mld_c_dec_map_bld.f90 mlprec/impl/mld_caggrmap_bld.f90 mlprec/impl/mld_ccoarse_bld.f90 mlprec/impl/mld_ccprecset.F90 mlprec/impl/mld_cprecinit.F90 mlprec/impl/mld_cprecset.F90 mlprec/impl/mld_d_dec_map_bld.f90 mlprec/impl/mld_daggrmap_bld.f90 mlprec/impl/mld_dcoarse_bld.f90 mlprec/impl/mld_dcprecset.F90 mlprec/impl/mld_dprecinit.F90 mlprec/impl/mld_dprecset.F90 mlprec/impl/mld_s_dec_map_bld.f90 mlprec/impl/mld_saggrmap_bld.f90 mlprec/impl/mld_scoarse_bld.f90 mlprec/impl/mld_scprecset.F90 mlprec/impl/mld_sprecinit.F90 mlprec/impl/mld_sprecset.F90 mlprec/impl/mld_z_dec_map_bld.f90 mlprec/impl/mld_zaggrmap_bld.f90 mlprec/impl/mld_zcoarse_bld.f90 mlprec/impl/mld_zcprecset.F90 mlprec/impl/mld_zprecinit.F90 mlprec/impl/mld_zprecset.F90 mlprec/mld_base_prec_type.F90 mlprec/mld_c_inner_mod.f90 mlprec/mld_c_onelev_mod.f90 mlprec/mld_d_base_solver_mod.f90 mlprec/mld_d_inner_mod.f90 mlprec/mld_d_onelev_mod.f90 mlprec/mld_d_prec_type.f90 mlprec/mld_d_umf_solver.F90 mlprec/mld_s_inner_mod.f90 mlprec/mld_s_onelev_mod.f90 mlprec/mld_z_inner_mod.f90 mlprec/mld_z_onelev_mod.f90 tests/pdegen/ppde2d.f90 tests/pdegen/ppde3d.f90 tests/pdegen/runs/ppde.inp tests/pdegen/spde2d.f90 tests/pdegen/spde3d.f90 Added option to apply ordering while aggregating. --- docs/src/userinterface.tex | 8 +++ mlprec/impl/level/mld_c_base_onelev_cseti.f90 | 3 + mlprec/impl/level/mld_c_base_onelev_seti.f90 | 3 + mlprec/impl/level/mld_d_base_onelev_cseti.f90 | 3 + mlprec/impl/level/mld_d_base_onelev_seti.f90 | 3 + mlprec/impl/level/mld_s_base_onelev_cseti.f90 | 3 + mlprec/impl/level/mld_s_base_onelev_seti.f90 | 3 + mlprec/impl/level/mld_z_base_onelev_cseti.f90 | 3 + mlprec/impl/level/mld_z_base_onelev_seti.f90 | 3 + mlprec/impl/mld_c_dec_map_bld.f90 | 51 ++++++++++------ mlprec/impl/mld_caggrmap_bld.f90 | 9 +-- mlprec/impl/mld_ccoarse_bld.f90 | 7 ++- mlprec/impl/mld_ccprecset.F90 | 21 ++++--- mlprec/impl/mld_cprecinit.F90 | 22 +++---- mlprec/impl/mld_cprecset.F90 | 39 ++++++++---- mlprec/impl/mld_d_dec_map_bld.f90 | 51 ++++++++++------ mlprec/impl/mld_daggrmap_bld.f90 | 9 +-- mlprec/impl/mld_dcoarse_bld.f90 | 4 +- mlprec/impl/mld_dcprecset.F90 | 60 +++++++++--------- mlprec/impl/mld_dprecinit.F90 | 16 ++--- mlprec/impl/mld_dprecset.F90 | 13 ++-- mlprec/impl/mld_s_dec_map_bld.f90 | 51 ++++++++++------ mlprec/impl/mld_saggrmap_bld.f90 | 9 +-- mlprec/impl/mld_scoarse_bld.f90 | 5 +- mlprec/impl/mld_scprecset.F90 | 18 +++--- mlprec/impl/mld_sprecinit.F90 | 22 +++---- mlprec/impl/mld_sprecset.F90 | 38 +++++++----- mlprec/impl/mld_z_dec_map_bld.f90 | 51 ++++++++++------ mlprec/impl/mld_zaggrmap_bld.f90 | 9 +-- mlprec/impl/mld_zcoarse_bld.f90 | 8 ++- mlprec/impl/mld_zcprecset.F90 | 61 ++++++++++--------- mlprec/impl/mld_zprecinit.F90 | 15 ++--- mlprec/impl/mld_zprecset.F90 | 36 +++++++---- mlprec/mld_base_prec_type.F90 | 54 +++++++++++----- mlprec/mld_c_inner_mod.f90 | 6 +- mlprec/mld_c_onelev_mod.f90 | 1 + mlprec/mld_d_base_solver_mod.f90 | 1 + mlprec/mld_d_inner_mod.f90 | 6 +- mlprec/mld_d_onelev_mod.f90 | 10 +-- mlprec/mld_d_prec_type.f90 | 6 +- mlprec/mld_d_umf_solver.F90 | 4 +- mlprec/mld_s_inner_mod.f90 | 6 +- mlprec/mld_s_onelev_mod.f90 | 1 + mlprec/mld_z_inner_mod.f90 | 6 +- mlprec/mld_z_onelev_mod.f90 | 1 + tests/pdegen/ppde2d.f90 | 4 ++ tests/pdegen/ppde3d.f90 | 5 +- tests/pdegen/runs/ppde.inp | 9 +-- tests/pdegen/spde2d.f90 | 4 ++ tests/pdegen/spde3d.f90 | 4 ++ 50 files changed, 489 insertions(+), 296 deletions(-) diff --git a/docs/src/userinterface.tex b/docs/src/userinterface.tex index 7771bfcc..0d03c64e 100644 --- a/docs/src/userinterface.tex +++ b/docs/src/userinterface.tex @@ -305,6 +305,14 @@ according to their needs. & \texttt{'DEC'} & Aggregation algorithm. Currently, only the decoupled aggregation is available. \\ \hline +\verb|mld_aggr_ord_| \break \verb|AGGR_ORD| & \verb|character(len=*)| + & \texttt{'NAT'} + & \texttt{'DEGREE'} + & Initial ordering of indices for aggregation + algorithm: natural ordering or sorted by + descending degree of the node in the + matrix graph. Since aggregation is + heuristics, results will be different. \\ \hline \verb|mld_aggr_kind_| \break \verb|AGGR_KIND| & \verb|character(len=*)| & \texttt{'SMOOTHED'} \hspace{2.5cm} \texttt{'NONSMOOTHED'} & \texttt{'SMOOTHED'} diff --git a/mlprec/impl/level/mld_c_base_onelev_cseti.f90 b/mlprec/impl/level/mld_c_base_onelev_cseti.f90 index e84c84e8..ecb4e6c9 100644 --- a/mlprec/impl/level/mld_c_base_onelev_cseti.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_cseti.f90 @@ -73,6 +73,9 @@ subroutine mld_c_base_onelev_cseti(lv,what,val,info) case ('AGGR_ALG') lv%parms%aggr_alg = val + case ('AGGR_ORD') + lv%parms%aggr_ord = val + case ('AGGR_KIND') lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_c_base_onelev_seti.f90 b/mlprec/impl/level/mld_c_base_onelev_seti.f90 index 56030073..859c1b5c 100644 --- a/mlprec/impl/level/mld_c_base_onelev_seti.f90 +++ b/mlprec/impl/level/mld_c_base_onelev_seti.f90 @@ -73,6 +73,9 @@ subroutine mld_c_base_onelev_seti(lv,what,val,info) case (mld_aggr_alg_) lv%parms%aggr_alg = val + case (mld_aggr_ord_) + lv%parms%aggr_ord = val + case (mld_aggr_kind_) lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_d_base_onelev_cseti.f90 b/mlprec/impl/level/mld_d_base_onelev_cseti.f90 index 39e32ae5..a2e3255b 100644 --- a/mlprec/impl/level/mld_d_base_onelev_cseti.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_cseti.f90 @@ -73,6 +73,9 @@ subroutine mld_d_base_onelev_cseti(lv,what,val,info) case ('AGGR_ALG') lv%parms%aggr_alg = val + case ('AGGR_ORD') + lv%parms%aggr_ord = val + case ('AGGR_KIND') lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_d_base_onelev_seti.f90 b/mlprec/impl/level/mld_d_base_onelev_seti.f90 index 3ca46d90..4e27cbca 100644 --- a/mlprec/impl/level/mld_d_base_onelev_seti.f90 +++ b/mlprec/impl/level/mld_d_base_onelev_seti.f90 @@ -73,6 +73,9 @@ subroutine mld_d_base_onelev_seti(lv,what,val,info) case (mld_aggr_alg_) lv%parms%aggr_alg = val + case (mld_aggr_ord_) + lv%parms%aggr_ord = val + case (mld_aggr_kind_) lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_s_base_onelev_cseti.f90 b/mlprec/impl/level/mld_s_base_onelev_cseti.f90 index 6f4bb091..b35e92ba 100644 --- a/mlprec/impl/level/mld_s_base_onelev_cseti.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_cseti.f90 @@ -73,6 +73,9 @@ subroutine mld_s_base_onelev_cseti(lv,what,val,info) case ('AGGR_ALG') lv%parms%aggr_alg = val + case ('AGGR_ORD') + lv%parms%aggr_ord = val + case ('AGGR_KIND') lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_s_base_onelev_seti.f90 b/mlprec/impl/level/mld_s_base_onelev_seti.f90 index 39acbdb0..ed687707 100644 --- a/mlprec/impl/level/mld_s_base_onelev_seti.f90 +++ b/mlprec/impl/level/mld_s_base_onelev_seti.f90 @@ -73,6 +73,9 @@ subroutine mld_s_base_onelev_seti(lv,what,val,info) case (mld_aggr_alg_) lv%parms%aggr_alg = val + case (mld_aggr_ord_) + lv%parms%aggr_ord = val + case (mld_aggr_kind_) lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_z_base_onelev_cseti.f90 b/mlprec/impl/level/mld_z_base_onelev_cseti.f90 index 3281f5ad..6352b4b9 100644 --- a/mlprec/impl/level/mld_z_base_onelev_cseti.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_cseti.f90 @@ -73,6 +73,9 @@ subroutine mld_z_base_onelev_cseti(lv,what,val,info) case ('AGGR_ALG') lv%parms%aggr_alg = val + case ('AGGR_ORD') + lv%parms%aggr_ord = val + case ('AGGR_KIND') lv%parms%aggr_kind = val diff --git a/mlprec/impl/level/mld_z_base_onelev_seti.f90 b/mlprec/impl/level/mld_z_base_onelev_seti.f90 index 61014f13..2d0a36a1 100644 --- a/mlprec/impl/level/mld_z_base_onelev_seti.f90 +++ b/mlprec/impl/level/mld_z_base_onelev_seti.f90 @@ -73,6 +73,9 @@ subroutine mld_z_base_onelev_seti(lv,what,val,info) case (mld_aggr_alg_) lv%parms%aggr_alg = val + case (mld_aggr_ord_) + lv%parms%aggr_ord = val + case (mld_aggr_kind_) lv%parms%aggr_kind = val diff --git a/mlprec/impl/mld_c_dec_map_bld.f90 b/mlprec/impl/mld_c_dec_map_bld.f90 index 687f7161..7d0c5f94 100644 --- a/mlprec/impl/mld_c_dec_map_bld.f90 +++ b/mlprec/impl/mld_c_dec_map_bld.f90 @@ -37,7 +37,7 @@ !!$ !!$ -subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_c_dec_map_bld @@ -45,6 +45,7 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -52,9 +53,11 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), intent(out) :: info ! Local variables - integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:) + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:) complex(psb_spk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + type(psb_c_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl logical :: recovery integer(psb_ipk_) :: debug_level, debug_unit,err_act @@ -75,7 +78,7 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) ncol = desc_a%get_local_cols() nr = a%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info) if(info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& @@ -90,11 +93,20 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + call a%cp_to(acsr) + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = acsr%irp(i+1) - acsr%irp(i) + end do + call acsr%free() + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i @@ -106,7 +118,8 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) nlp = 0 do icnt = 0 - do i=1, nr + do ii=1, nr + i = idxs(ii) if (ilaggr(i) == -(nr+1)) then ! ! 1. Untouched nodes are marked >0 together @@ -124,14 +137,16 @@ subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) end if do k=1, nz - j = icol(k) - ilg = ilaggr(j) - if ((ilg<0).and.(1<=j).and.(j<=nr).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif + j = icol(k) + if ((1<=j).and.(j<=nr)) then + ilg = ilaggr(j) + if ((ilg<0).and.(i /= j)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ilaggr(j) = naggr + else + ilaggr(j) = -naggr + endif + end if end if enddo diff --git a/mlprec/impl/mld_caggrmap_bld.f90 b/mlprec/impl/mld_caggrmap_bld.f90 index 81e8be11..e2ed909b 100644 --- a/mlprec/impl/mld_caggrmap_bld.f90 +++ b/mlprec/impl/mld_caggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod use mld_c_inner_mod, mld_protect_name => mld_caggrmap_bld @@ -87,6 +87,7 @@ subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_spk_), intent(in) :: theta type(psb_cspmat_type), intent(in) :: a @@ -118,7 +119,7 @@ subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) select case (aggr_type) case (mld_dec_aggr_) - call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + call mld_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) case (mld_sym_dec_aggr_) nr = a%get_nrows() @@ -126,14 +127,14 @@ subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) & rscale=.false.,cscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) - if (info == psb_success_) call atrans%transp(atmp) + if (info == psb_success_) call atmp%transp(atrans) if (info == psb_success_) call atrans%cscnv(info,type='COO') if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') - if (info == psb_success_) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call mld_dec_map_bld(iorder,theta,atmp,desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() case default diff --git a/mlprec/impl/mld_ccoarse_bld.f90 b/mlprec/impl/mld_ccoarse_bld.f90 index 9f184722..ab9fc56d 100644 --- a/mlprec/impl/mld_ccoarse_bld.f90 +++ b/mlprec/impl/mld_ccoarse_bld.f90 @@ -39,7 +39,7 @@ ! File: mld_ccoarse_bld.f90 ! ! Subroutine: mld_ccoarse_bld -! Version: complex +! Version: real ! ! This routine builds the matrix associated to the current level of the ! multilevel preconditioner from the matrix associated to the previous level, @@ -95,6 +95,8 @@ subroutine mld_ccoarse_bld(a,desc_a,p,info) & mld_mult_ml_,is_legal_ml_type) call mld_check_def(p%parms%aggr_alg,'Aggregation',& & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(p%parms%aggr_kind,'Smoother',& & mld_smooth_prol_,is_legal_ml_aggr_kind) call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& @@ -117,7 +119,7 @@ subroutine mld_ccoarse_bld(a,desc_a,p,info) ! aggregation algorithm. This also defines a tentative prolongator from ! the coarse to the fine level. ! - call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_thresh,& + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& & a,desc_a,ilaggr,nlaggr,info) if (info /= psb_success_) then @@ -136,6 +138,7 @@ subroutine mld_ccoarse_bld(a,desc_a,p,info) call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') goto 9999 end if + ! ! Fix the base_a and base_desc pointers for handling of residuals. ! This is correct because this routine is only called at levels >=2. diff --git a/mlprec/impl/mld_ccprecset.F90 b/mlprec/impl/mld_ccprecset.F90 index 29235f51..e4858649 100644 --- a/mlprec/impl/mld_ccprecset.F90 +++ b/mlprec/impl/mld_ccprecset.F90 @@ -93,6 +93,7 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) use mld_c_mumps_solver #endif + implicit none ! Arguments @@ -148,9 +149,10 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN') call p%precv(ilev_)%set(what,val,info) @@ -166,9 +168,10 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN',& & 'COARSE_MAT') @@ -214,7 +217,7 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) case(mld_mumps_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) - call p%precv(nlev_)%set('COARSE_MAT',mld_distr_mat_,info) + call p%precv(nlev_)%set('COARSE_MAT',mld_distr_mat_,info) case(mld_jac_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),mld_diag_scale_,info) @@ -280,7 +283,7 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case('ML_TYPE','AGGR_ALG','AGGR_KIND',& + case('ML_TYPE','AGGR_ALG','AGGR_ORD','AGGR_KIND',& & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& & 'SMOOTHER_POS','AGGR_OMEGA_ALG',& & 'AGGR_EIG','AGGR_FILTER') @@ -302,7 +305,6 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) #if defined(HAVE_SLU_) call onelev_set_solver(p%precv(nlev_),mld_slu_,info) - call onelev_set_solver(p%precv(nlev_),mld_slu_,info) #elif defined(HAVE_MUMPS_) call onelev_set_solver(p%precv(nlev_),mld_mumps_,info) #else @@ -321,7 +323,6 @@ subroutine mld_ccprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set('COARSE_MAT',mld_distr_mat_,info) - case(mld_jac_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),mld_diag_scale_,info) diff --git a/mlprec/impl/mld_cprecinit.F90 b/mlprec/impl/mld_cprecinit.F90 index a5f6ffd2..b1347ea7 100644 --- a/mlprec/impl/mld_cprecinit.F90 +++ b/mlprec/impl/mld_cprecinit.F90 @@ -2,9 +2,9 @@ !!$ !!$ MLD2P4 version 2.0 !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) !!$ -!!$ (C) Copyright 2008,2009,2010,2012,2013 +!!$ (C) Copyright 2008, 2010, 2012, 2015 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -61,7 +61,7 @@ ! ! 'ML' - Multilevel hybrid preconditioner (additive on the ! same level and multiplicative through the levels), -! with 2 levels and post-smoothing only. RAS with +! with 2 levels, pre and post-smoothing, RAS with ! overlap 1 and ILU(0) on the local blocks is ! applied as post-smoother at each level, but the ! coarsest one; four sweeps of the block-Jacobi solver, @@ -97,9 +97,6 @@ subroutine mld_cprecinit(p,ptype,info,nlev) use mld_c_id_solver use mld_c_diag_solver use mld_c_ilu_solver -#if defined(HAVE_UMF_) && 0 - use mld_c_umf_solver -#endif #if defined(HAVE_SLU_) use mld_c_slu_solver #endif @@ -115,7 +112,7 @@ subroutine mld_cprecinit(p,ptype,info,nlev) ! Local variables integer(psb_ipk_) :: nlev_, ilev_ - real(psb_spk_) :: thr + real(psb_spk_) :: thr, scale character(len=*), parameter :: name='mld_precinit' info = psb_success_ @@ -191,10 +188,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev) ilev_ = nlev_ allocate(mld_c_jac_smoother_type :: p%precv(ilev_)%sm, stat=info) if (info /= psb_success_) return -#if defined(HAVE_UMF_) && 0 - allocate(mld_c_umf_solver_type :: p%precv(ilev_)%sm%sv, stat=info) -#elif defined(HAVE_SLU_) - allocate(mld_c_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) +#if defined(HAVE_SLU_) + allocate(mld_c_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #else allocate(mld_c_ilu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #endif @@ -205,10 +200,11 @@ subroutine mld_cprecinit(p,ptype,info,nlev) call p%precv(ilev_)%set(mld_sub_prol_,psb_none_,info) call p%precv(ilev_)%set(mld_sub_ovr_,izero,info) - thr = 0.16d0 + thr = 0.05 + scale = 1.0 do ilev_=1,nlev_ call p%precv(ilev_)%set(mld_aggr_thresh_,thr,info) - thr = thr/2 + call p%precv(ilev_)%set(mld_aggr_scale_,scale,info) end do case default diff --git a/mlprec/impl/mld_cprecset.F90 b/mlprec/impl/mld_cprecset.F90 index 7e16077d..6728ca12 100644 --- a/mlprec/impl/mld_cprecset.F90 +++ b/mlprec/impl/mld_cprecset.F90 @@ -92,6 +92,7 @@ subroutine mld_cprecseti(p,what,val,info,ilev) #if defined(HAVE_MUMPS_) use mld_c_mumps_solver #endif + implicit none ! Arguments @@ -148,8 +149,8 @@ subroutine mld_cprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) @@ -166,8 +167,8 @@ subroutine mld_cprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& @@ -207,7 +208,7 @@ subroutine mld_cprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set(mld_coarse_mat_,mld_repl_mat_,info) - case(mld_sludist_, mld_mumps_) + case(mld_sludist_,mld_mumps_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set(mld_coarse_mat_,mld_distr_mat_,info) @@ -276,7 +277,7 @@ subroutine mld_cprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& + case(mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,mld_aggr_kind_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_smoother_pos_,mld_aggr_omega_alg_,& & mld_aggr_eig_,mld_aggr_filter_) @@ -298,8 +299,8 @@ subroutine mld_cprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) #if defined(HAVE_SLU_) call onelev_set_solver(p%precv(nlev_),mld_slu_,info) -#elif defined(HAVE_MUMPS_) - call onelev_set_solver(p%precv(nlev_),mld_mumps_,info) +#elif defined(HAVE_SLU_) + call onelev_set_solver(p%precv(nlev_),mld_slu_,info) #else call onelev_set_solver(p%precv(nlev_),mld_ilu_n_,info) #endif @@ -599,7 +600,6 @@ contains end if end if #endif - case default ! ! Do nothing and hope for the best :) @@ -723,13 +723,26 @@ subroutine mld_cprecsetsv(p,val,info,ilev) do ilev_ = ilmin, ilmax if (allocated(p%precv(ilev_)%sm)) then - if (allocated(p%precv(ilev_)%sm%sv)) & - & deallocate(p%precv(ilev_)%sm%sv) + if (allocated(p%precv(ilev_)%sm%sv)) then + if (.not.same_type_as(p%precv(ilev_)%sm%sv,val)) then + deallocate(p%precv(ilev_)%sm%sv,stat=info) + if (info /= 0) then + info = 3111 + return + end if + end if + if (.not.allocated(p%precv(ilev_)%sm%sv)) then #ifdef HAVE_MOLD - allocate(p%precv(ilev_)%sm%sv,mold=val) + allocate(p%precv(ilev_)%sm%sv,mold=val,stat=info) #else - allocate(p%precv(ilev_)%sm%sv,source=val) + allocate(p%precv(ilev_)%sm%sv,source=val,stat=info) #endif + if (info /= 0) then + info = 3111 + return + end if + end if + end if call p%precv(ilev_)%sm%sv%default() else info = 3111 diff --git a/mlprec/impl/mld_d_dec_map_bld.f90 b/mlprec/impl/mld_d_dec_map_bld.f90 index bc0d40bf..67e4437b 100644 --- a/mlprec/impl/mld_d_dec_map_bld.f90 +++ b/mlprec/impl/mld_d_dec_map_bld.f90 @@ -37,7 +37,7 @@ !!$ !!$ -subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_d_dec_map_bld @@ -45,6 +45,7 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -52,9 +53,11 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), intent(out) :: info ! Local variables - integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:) + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:) real(psb_dpk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + type(psb_d_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl logical :: recovery integer(psb_ipk_) :: debug_level, debug_unit,err_act @@ -75,7 +78,7 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) ncol = desc_a%get_local_cols() nr = a%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info) if(info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& @@ -90,11 +93,20 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + call a%cp_to(acsr) + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = acsr%irp(i+1) - acsr%irp(i) + end do + call acsr%free() + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i @@ -106,7 +118,8 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) nlp = 0 do icnt = 0 - do i=1, nr + do ii=1, nr + i = idxs(ii) if (ilaggr(i) == -(nr+1)) then ! ! 1. Untouched nodes are marked >0 together @@ -124,14 +137,16 @@ subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) end if do k=1, nz - j = icol(k) - ilg = ilaggr(j) - if ((ilg<0).and.(1<=j).and.(j<=nr).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif + j = icol(k) + if ((1<=j).and.(j<=nr)) then + ilg = ilaggr(j) + if ((ilg<0).and.(i /= j)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ilaggr(j) = naggr + else + ilaggr(j) = -naggr + endif + end if end if enddo diff --git a/mlprec/impl/mld_daggrmap_bld.f90 b/mlprec/impl/mld_daggrmap_bld.f90 index 025622f3..c652983a 100644 --- a/mlprec/impl/mld_daggrmap_bld.f90 +++ b/mlprec/impl/mld_daggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod use mld_d_inner_mod, mld_protect_name => mld_daggrmap_bld @@ -87,6 +87,7 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_dpk_), intent(in) :: theta type(psb_dspmat_type), intent(in) :: a @@ -118,7 +119,7 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) select case (aggr_type) case (mld_dec_aggr_) - call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + call mld_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) case (mld_sym_dec_aggr_) nr = a%get_nrows() @@ -126,14 +127,14 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) & rscale=.false.,cscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) - if (info == psb_success_) call atrans%transp(atmp) + if (info == psb_success_) call atmp%transp(atrans) if (info == psb_success_) call atrans%cscnv(info,type='COO') if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') - if (info == psb_success_) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call mld_dec_map_bld(iorder,theta,atmp,desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() case default diff --git a/mlprec/impl/mld_dcoarse_bld.f90 b/mlprec/impl/mld_dcoarse_bld.f90 index 8d8b5e86..b52c4215 100644 --- a/mlprec/impl/mld_dcoarse_bld.f90 +++ b/mlprec/impl/mld_dcoarse_bld.f90 @@ -95,6 +95,8 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) & mld_mult_ml_,is_legal_ml_type) call mld_check_def(p%parms%aggr_alg,'Aggregation',& & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(p%parms%aggr_kind,'Smoother',& & mld_smooth_prol_,is_legal_ml_aggr_kind) call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& @@ -117,7 +119,7 @@ subroutine mld_dcoarse_bld(a,desc_a,p,info) ! aggregation algorithm. This also defines a tentative prolongator from ! the coarse to the fine level. ! - call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_thresh,& + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& & a,desc_a,ilaggr,nlaggr,info) if (info /= psb_success_) then diff --git a/mlprec/impl/mld_dcprecset.F90 b/mlprec/impl/mld_dcprecset.F90 index 13d18ff6..47782d44 100644 --- a/mlprec/impl/mld_dcprecset.F90 +++ b/mlprec/impl/mld_dcprecset.F90 @@ -99,6 +99,7 @@ subroutine mld_dcprecseti(p,what,val,info,ilev) use mld_d_mumps_solver #endif + implicit none ! Arguments @@ -113,6 +114,7 @@ subroutine mld_dcprecseti(p,what,val,info,ilev) character(len=*), parameter :: name='mld_precseti' info = psb_success_ + if (.not.allocated(p%precv)) then info = 3111 write(psb_err_unit,*) name,& @@ -153,9 +155,10 @@ subroutine mld_dcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN') call p%precv(ilev_)%set(what,val,info) @@ -171,9 +174,10 @@ subroutine mld_dcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN',& & 'COARSE_MAT') @@ -287,7 +291,7 @@ subroutine mld_dcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case('ML_TYPE','AGGR_ALG','AGGR_KIND',& + case('ML_TYPE','AGGR_ALG','AGGR_ORD','AGGR_KIND',& & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& & 'SMOOTHER_POS','AGGR_OMEGA_ALG',& & 'AGGR_EIG','AGGR_FILTER') @@ -592,6 +596,27 @@ contains info = -5 end if #endif +#ifdef HAVE_MUMPS_ + case (mld_mumps_) + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_d_mumps_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_d_mumps_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_d_mumps_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if +#endif + #ifdef HAVE_UMF_ case (mld_umf_) if (allocated(level%sm)) then @@ -642,27 +667,6 @@ contains info = -5 end if #endif - -#ifdef HAVE_MUMPS_ - case (mld_mumps_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_d_mumps_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_d_mumps_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_d_mumps_solver_type :: level%sm%sv, stat=info) - endif - if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() - end if -#endif case default ! ! Do nothing and hope for the best :) diff --git a/mlprec/impl/mld_dprecinit.F90 b/mlprec/impl/mld_dprecinit.F90 index bc8e6ced..86be5b05 100644 --- a/mlprec/impl/mld_dprecinit.F90 +++ b/mlprec/impl/mld_dprecinit.F90 @@ -2,9 +2,9 @@ !!$ !!$ MLD2P4 version 2.0 !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) !!$ -!!$ (C) Copyright 2008,2009,2010,2012,2013 +!!$ (C) Copyright 2008, 2010, 2012, 2015 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -61,7 +61,7 @@ ! ! 'ML' - Multilevel hybrid preconditioner (additive on the ! same level and multiplicative through the levels), -! with 2 levels and post-smoothing only. RAS with +! with 2 levels, pre and post-smoothing, RAS with ! overlap 1 and ILU(0) on the local blocks is ! applied as post-smoother at each level, but the ! coarsest one; four sweeps of the block-Jacobi solver, @@ -104,6 +104,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev) use mld_d_slu_solver #endif + implicit none ! Arguments @@ -114,7 +115,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev) ! Local variables integer(psb_ipk_) :: nlev_, ilev_ - real(psb_dpk_) :: thr + real(psb_dpk_) :: thr, scale character(len=*), parameter :: name='mld_precinit' info = psb_success_ @@ -193,7 +194,7 @@ subroutine mld_dprecinit(p,ptype,info,nlev) #if defined(HAVE_UMF_) allocate(mld_d_umf_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #elif defined(HAVE_SLU_) - allocate(mld_d_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) + allocate(mld_d_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #else allocate(mld_d_ilu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #endif @@ -204,10 +205,11 @@ subroutine mld_dprecinit(p,ptype,info,nlev) call p%precv(ilev_)%set(mld_sub_prol_,psb_none_,info) call p%precv(ilev_)%set(mld_sub_ovr_,izero,info) - thr = 0.16d0 + thr = 0.05 + scale = 1.0 do ilev_=1,nlev_ call p%precv(ilev_)%set(mld_aggr_thresh_,thr,info) - thr = thr/2 + call p%precv(ilev_)%set(mld_aggr_scale_,scale,info) end do case default diff --git a/mlprec/impl/mld_dprecset.F90 b/mlprec/impl/mld_dprecset.F90 index 2f328351..1a5fe0ef 100644 --- a/mlprec/impl/mld_dprecset.F90 +++ b/mlprec/impl/mld_dprecset.F90 @@ -155,8 +155,8 @@ subroutine mld_dprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) @@ -173,8 +173,8 @@ subroutine mld_dprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& @@ -285,7 +285,7 @@ subroutine mld_dprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& + case(mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,mld_aggr_kind_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_smoother_pos_,mld_aggr_omega_alg_,& & mld_aggr_eig_,mld_aggr_filter_) @@ -309,8 +309,6 @@ subroutine mld_dprecseti(p,what,val,info,ilev) call onelev_set_solver(p%precv(nlev_),mld_umf_,info) #elif defined(HAVE_SLU_) call onelev_set_solver(p%precv(nlev_),mld_slu_,info) -#elif defined(HAVE_MUMPS_) - call onelev_set_solver(p%precv(nlev_),mld_mumps_,info) #else call onelev_set_solver(p%precv(nlev_),mld_ilu_n_,info) #endif @@ -1002,6 +1000,7 @@ subroutine mld_dprecsetr(p,what,val,info,ilev) case(mld_coarse_iluthrs_) ilev_=nlev_ call p%precv(ilev_)%set(mld_sub_iluthrs_,val,info) + case(mld_aggr_thresh_) thr = val do ilev_ = 2, nlev_ diff --git a/mlprec/impl/mld_s_dec_map_bld.f90 b/mlprec/impl/mld_s_dec_map_bld.f90 index 2fa3a6c3..e980ecf8 100644 --- a/mlprec/impl/mld_s_dec_map_bld.f90 +++ b/mlprec/impl/mld_s_dec_map_bld.f90 @@ -37,7 +37,7 @@ !!$ !!$ -subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_s_dec_map_bld @@ -45,6 +45,7 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta @@ -52,9 +53,11 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), intent(out) :: info ! Local variables - integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:) + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:) real(psb_spk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + type(psb_s_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl logical :: recovery integer(psb_ipk_) :: debug_level, debug_unit,err_act @@ -75,7 +78,7 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) ncol = desc_a%get_local_cols() nr = a%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info) if(info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& @@ -90,11 +93,20 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + call a%cp_to(acsr) + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = acsr%irp(i+1) - acsr%irp(i) + end do + call acsr%free() + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i @@ -106,7 +118,8 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) nlp = 0 do icnt = 0 - do i=1, nr + do ii=1, nr + i = idxs(ii) if (ilaggr(i) == -(nr+1)) then ! ! 1. Untouched nodes are marked >0 together @@ -124,14 +137,16 @@ subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) end if do k=1, nz - j = icol(k) - ilg = ilaggr(j) - if ((ilg<0).and.(1<=j).and.(j<=nr).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif + j = icol(k) + if ((1<=j).and.(j<=nr)) then + ilg = ilaggr(j) + if ((ilg<0).and.(i /= j)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ilaggr(j) = naggr + else + ilaggr(j) = -naggr + endif + end if end if enddo diff --git a/mlprec/impl/mld_saggrmap_bld.f90 b/mlprec/impl/mld_saggrmap_bld.f90 index 4bc55695..b46bb98a 100644 --- a/mlprec/impl/mld_saggrmap_bld.f90 +++ b/mlprec/impl/mld_saggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod use mld_s_inner_mod, mld_protect_name => mld_saggrmap_bld @@ -87,6 +87,7 @@ subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_spk_), intent(in) :: theta type(psb_sspmat_type), intent(in) :: a @@ -118,7 +119,7 @@ subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) select case (aggr_type) case (mld_dec_aggr_) - call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + call mld_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) case (mld_sym_dec_aggr_) nr = a%get_nrows() @@ -126,14 +127,14 @@ subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) & rscale=.false.,cscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) - if (info == psb_success_) call atrans%transp(atmp) + if (info == psb_success_) call atmp%transp(atrans) if (info == psb_success_) call atrans%cscnv(info,type='COO') if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') - if (info == psb_success_) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call mld_dec_map_bld(iorder,theta,atmp,desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() case default diff --git a/mlprec/impl/mld_scoarse_bld.f90 b/mlprec/impl/mld_scoarse_bld.f90 index 344bfc18..acc85154 100644 --- a/mlprec/impl/mld_scoarse_bld.f90 +++ b/mlprec/impl/mld_scoarse_bld.f90 @@ -95,6 +95,8 @@ subroutine mld_scoarse_bld(a,desc_a,p,info) & mld_mult_ml_,is_legal_ml_type) call mld_check_def(p%parms%aggr_alg,'Aggregation',& & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(p%parms%aggr_kind,'Smoother',& & mld_smooth_prol_,is_legal_ml_aggr_kind) call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& @@ -117,7 +119,7 @@ subroutine mld_scoarse_bld(a,desc_a,p,info) ! aggregation algorithm. This also defines a tentative prolongator from ! the coarse to the fine level. ! - call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_thresh,& + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& & a,desc_a,ilaggr,nlaggr,info) if (info /= psb_success_) then @@ -136,6 +138,7 @@ subroutine mld_scoarse_bld(a,desc_a,p,info) call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') goto 9999 end if + ! ! Fix the base_a and base_desc pointers for handling of residuals. ! This is correct because this routine is only called at levels >=2. diff --git a/mlprec/impl/mld_scprecset.F90 b/mlprec/impl/mld_scprecset.F90 index d1c651f6..5b80da49 100644 --- a/mlprec/impl/mld_scprecset.F90 +++ b/mlprec/impl/mld_scprecset.F90 @@ -93,6 +93,7 @@ subroutine mld_scprecseti(p,what,val,info,ilev) use mld_s_mumps_solver #endif + implicit none ! Arguments @@ -148,9 +149,10 @@ subroutine mld_scprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN') call p%precv(ilev_)%set(what,val,info) @@ -166,9 +168,10 @@ subroutine mld_scprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN',& & 'COARSE_MAT') @@ -280,7 +283,7 @@ subroutine mld_scprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case('ML_TYPE','AGGR_ALG','AGGR_KIND',& + case('ML_TYPE','AGGR_ALG','AGGR_ORD','AGGR_KIND',& & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& & 'SMOOTHER_POS','AGGR_OMEGA_ALG',& & 'AGGR_EIG','AGGR_FILTER') @@ -583,7 +586,6 @@ contains info = -5 end if #endif - #ifdef HAVE_MUMPS_ case (mld_mumps_) if (allocated(level%sm%sv)) then diff --git a/mlprec/impl/mld_sprecinit.F90 b/mlprec/impl/mld_sprecinit.F90 index 3ff111ba..24021e6d 100644 --- a/mlprec/impl/mld_sprecinit.F90 +++ b/mlprec/impl/mld_sprecinit.F90 @@ -2,9 +2,9 @@ !!$ !!$ MLD2P4 version 2.0 !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) !!$ -!!$ (C) Copyright 2008,2009,2010,2012,2013 +!!$ (C) Copyright 2008, 2010, 2012, 2015 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -61,7 +61,7 @@ ! ! 'ML' - Multilevel hybrid preconditioner (additive on the ! same level and multiplicative through the levels), -! with 2 levels and post-smoothing only. RAS with +! with 2 levels, pre and post-smoothing, RAS with ! overlap 1 and ILU(0) on the local blocks is ! applied as post-smoother at each level, but the ! coarsest one; four sweeps of the block-Jacobi solver, @@ -97,9 +97,6 @@ subroutine mld_sprecinit(p,ptype,info,nlev) use mld_s_id_solver use mld_s_diag_solver use mld_s_ilu_solver -#if defined(HAVE_UMF_) && 0 - use mld_s_umf_solver -#endif #if defined(HAVE_SLU_) use mld_s_slu_solver #endif @@ -115,7 +112,7 @@ subroutine mld_sprecinit(p,ptype,info,nlev) ! Local variables integer(psb_ipk_) :: nlev_, ilev_ - real(psb_spk_) :: thr + real(psb_spk_) :: thr, scale character(len=*), parameter :: name='mld_precinit' info = psb_success_ @@ -191,10 +188,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev) ilev_ = nlev_ allocate(mld_s_jac_smoother_type :: p%precv(ilev_)%sm, stat=info) if (info /= psb_success_) return -#if defined(HAVE_UMF_) && 0 - allocate(mld_s_umf_solver_type :: p%precv(ilev_)%sm%sv, stat=info) -#elif defined(HAVE_SLU_) - allocate(mld_s_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) +#if defined(HAVE_SLU_) + allocate(mld_s_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #else allocate(mld_s_ilu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #endif @@ -205,10 +200,11 @@ subroutine mld_sprecinit(p,ptype,info,nlev) call p%precv(ilev_)%set(mld_sub_prol_,psb_none_,info) call p%precv(ilev_)%set(mld_sub_ovr_,izero,info) - thr = 0.16d0 + thr = 0.05 + scale = 1.0 do ilev_=1,nlev_ call p%precv(ilev_)%set(mld_aggr_thresh_,thr,info) - thr = thr/2 + call p%precv(ilev_)%set(mld_aggr_scale_,scale,info) end do case default diff --git a/mlprec/impl/mld_sprecset.F90 b/mlprec/impl/mld_sprecset.F90 index c7d277d9..822e1ea2 100644 --- a/mlprec/impl/mld_sprecset.F90 +++ b/mlprec/impl/mld_sprecset.F90 @@ -149,8 +149,8 @@ subroutine mld_sprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) @@ -167,8 +167,8 @@ subroutine mld_sprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& @@ -277,7 +277,7 @@ subroutine mld_sprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& + case(mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,mld_aggr_kind_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_smoother_pos_,mld_aggr_omega_alg_,& & mld_aggr_eig_,mld_aggr_filter_) @@ -299,8 +299,8 @@ subroutine mld_sprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) #if defined(HAVE_SLU_) call onelev_set_solver(p%precv(nlev_),mld_slu_,info) -#elif defined(HAVE_MUMPS_) - call onelev_set_solver(p%precv(nlev_),mld_mumps_,info) +#elif defined(HAVE_SLU_) + call onelev_set_solver(p%precv(nlev_),mld_slu_,info) #else call onelev_set_solver(p%precv(nlev_),mld_ilu_n_,info) #endif @@ -579,7 +579,6 @@ contains info = -5 end if #endif - #ifdef HAVE_MUMPS_ case (mld_mumps_) if (allocated(level%sm%sv)) then @@ -601,8 +600,6 @@ contains end if end if #endif - - case default ! ! Do nothing and hope for the best :) @@ -726,13 +723,26 @@ subroutine mld_sprecsetsv(p,val,info,ilev) do ilev_ = ilmin, ilmax if (allocated(p%precv(ilev_)%sm)) then - if (allocated(p%precv(ilev_)%sm%sv)) & - & deallocate(p%precv(ilev_)%sm%sv) + if (allocated(p%precv(ilev_)%sm%sv)) then + if (.not.same_type_as(p%precv(ilev_)%sm%sv,val)) then + deallocate(p%precv(ilev_)%sm%sv,stat=info) + if (info /= 0) then + info = 3111 + return + end if + end if + if (.not.allocated(p%precv(ilev_)%sm%sv)) then #ifdef HAVE_MOLD - allocate(p%precv(ilev_)%sm%sv,mold=val) + allocate(p%precv(ilev_)%sm%sv,mold=val,stat=info) #else - allocate(p%precv(ilev_)%sm%sv,source=val) + allocate(p%precv(ilev_)%sm%sv,source=val,stat=info) #endif + if (info /= 0) then + info = 3111 + return + end if + end if + end if call p%precv(ilev_)%sm%sv%default() else info = 3111 diff --git a/mlprec/impl/mld_z_dec_map_bld.f90 b/mlprec/impl/mld_z_dec_map_bld.f90 index b578cc78..718856aa 100644 --- a/mlprec/impl/mld_z_dec_map_bld.f90 +++ b/mlprec/impl/mld_z_dec_map_bld.f90 @@ -37,7 +37,7 @@ !!$ !!$ -subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_z_dec_map_bld @@ -45,6 +45,7 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta @@ -52,9 +53,11 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), intent(out) :: info ! Local variables - integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:) + integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& + & ideg(:), idxs(:) complex(psb_dpk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + type(psb_z_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl logical :: recovery integer(psb_ipk_) :: debug_level, debug_unit,err_act @@ -75,7 +78,7 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) ncol = desc_a%get_local_cols() nr = a%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) + allocate(ilaggr(nr),neigh(nr),ideg(nr),idxs(nr),stat=info) if(info /= psb_success_) then info=psb_err_alloc_request_ call psb_errpush(info,name,i_err=(/2*nr,izero,izero,izero,izero/),& @@ -90,11 +93,20 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - + if (iorder == mld_aggr_ord_nat_) then + do i=1, nr + ilaggr(i) = -(nr+1) + idxs(i) = i + end do + else + call a%cp_to(acsr) + do i=1, nr + ilaggr(i) = -(nr+1) + ideg(i) = acsr%irp(i+1) - acsr%irp(i) + end do + call acsr%free() + call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) + end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i @@ -106,7 +118,8 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) nlp = 0 do icnt = 0 - do i=1, nr + do ii=1, nr + i = idxs(ii) if (ilaggr(i) == -(nr+1)) then ! ! 1. Untouched nodes are marked >0 together @@ -124,14 +137,16 @@ subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) end if do k=1, nz - j = icol(k) - ilg = ilaggr(j) - if ((ilg<0).and.(1<=j).and.(j<=nr).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif + j = icol(k) + if ((1<=j).and.(j<=nr)) then + ilg = ilaggr(j) + if ((ilg<0).and.(i /= j)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ilaggr(j) = naggr + else + ilaggr(j) = -naggr + endif + end if end if enddo diff --git a/mlprec/impl/mld_zaggrmap_bld.f90 b/mlprec/impl/mld_zaggrmap_bld.f90 index 41e6f677..2058fd66 100644 --- a/mlprec/impl/mld_zaggrmap_bld.f90 +++ b/mlprec/impl/mld_zaggrmap_bld.f90 @@ -79,7 +79,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) +subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod use mld_z_inner_mod, mld_protect_name => mld_zaggrmap_bld @@ -87,6 +87,7 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) implicit none ! Arguments + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_dpk_), intent(in) :: theta type(psb_zspmat_type), intent(in) :: a @@ -118,7 +119,7 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) select case (aggr_type) case (mld_dec_aggr_) - call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + call mld_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) case (mld_sym_dec_aggr_) nr = a%get_nrows() @@ -126,14 +127,14 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) & rscale=.false.,cscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) - if (info == psb_success_) call atrans%transp(atmp) + if (info == psb_success_) call atmp%transp(atrans) if (info == psb_success_) call atrans%cscnv(info,type='COO') if (info == psb_success_) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) call atmp%set_nrows(nr) call atmp%set_ncols(nr) if (info == psb_success_) call atrans%free() if (info == psb_success_) call atmp%cscnv(info,type='CSR') - if (info == psb_success_) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == psb_success_) call mld_dec_map_bld(iorder,theta,atmp,desc_a,nlaggr,ilaggr,info) if (info == psb_success_) call atmp%free() case default diff --git a/mlprec/impl/mld_zcoarse_bld.f90 b/mlprec/impl/mld_zcoarse_bld.f90 index dc67d86f..644284de 100644 --- a/mlprec/impl/mld_zcoarse_bld.f90 +++ b/mlprec/impl/mld_zcoarse_bld.f90 @@ -39,7 +39,7 @@ ! File: mld_zcoarse_bld.f90 ! ! Subroutine: mld_zcoarse_bld -! Version: complex +! Version: real ! ! This routine builds the matrix associated to the current level of the ! multilevel preconditioner from the matrix associated to the previous level, @@ -95,6 +95,8 @@ subroutine mld_zcoarse_bld(a,desc_a,p,info) & mld_mult_ml_,is_legal_ml_type) call mld_check_def(p%parms%aggr_alg,'Aggregation',& & mld_dec_aggr_,is_legal_ml_aggr_alg) + call mld_check_def(p%parms%aggr_ord,'Ordering',& + & mld_aggr_ord_nat_,is_legal_ml_aggr_ord) call mld_check_def(p%parms%aggr_kind,'Smoother',& & mld_smooth_prol_,is_legal_ml_aggr_kind) call mld_check_def(p%parms%coarse_mat,'Coarse matrix',& @@ -117,7 +119,7 @@ subroutine mld_zcoarse_bld(a,desc_a,p,info) ! aggregation algorithm. This also defines a tentative prolongator from ! the coarse to the fine level. ! - call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_thresh,& + call mld_aggrmap_bld(p%parms%aggr_alg,p%parms%aggr_ord,p%parms%aggr_thresh,& & a,desc_a,ilaggr,nlaggr,info) if (info /= psb_success_) then @@ -136,6 +138,7 @@ subroutine mld_zcoarse_bld(a,desc_a,p,info) call psb_errpush(psb_err_from_subroutine_,name,a_err='mld_aggrmat_asb') goto 9999 end if + ! ! Fix the base_a and base_desc pointers for handling of residuals. ! This is correct because this routine is only called at levels >=2. @@ -147,7 +150,6 @@ subroutine mld_zcoarse_bld(a,desc_a,p,info) return 9999 call psb_error_handler(err_act) - return end subroutine mld_zcoarse_bld diff --git a/mlprec/impl/mld_zcprecset.F90 b/mlprec/impl/mld_zcprecset.F90 index 0c16e732..c75b80ce 100644 --- a/mlprec/impl/mld_zcprecset.F90 +++ b/mlprec/impl/mld_zcprecset.F90 @@ -99,6 +99,7 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) use mld_z_mumps_solver #endif + implicit none ! Arguments @@ -154,9 +155,10 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN') call p%precv(ilev_)%set(what,val,info) @@ -172,9 +174,10 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case('SUB_SOLVE') call onelev_set_solver(p%precv(ilev_),val,info) - case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_KIND',& - & 'SMOOTHER_POS','AGGR_OMEGA_ALG','AGGR_EIG',& - & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& + case('SMOOTHER_SWEEPS','ML_TYPE','AGGR_ALG','AGGR_ORD',& + & 'AGGR_KIND','SMOOTHER_POS','AGGR_OMEGA_ALG',& + & 'AGGR_EIG','SMOOTHER_SWEEPS_PRE',& + & 'SMOOTHER_SWEEPS_POST',& & 'SUB_RESTR','SUB_PROL', & & 'SUB_REN','SUB_OVR','SUB_FILLIN',& & 'COARSE_MAT') @@ -288,7 +291,7 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case('ML_TYPE','AGGR_ALG','AGGR_KIND',& + case('ML_TYPE','AGGR_ALG','AGGR_ORD','AGGR_KIND',& & 'SMOOTHER_SWEEPS_PRE','SMOOTHER_SWEEPS_POST',& & 'SMOOTHER_POS','AGGR_OMEGA_ALG',& & 'AGGR_EIG','AGGR_FILTER') @@ -330,7 +333,6 @@ subroutine mld_zcprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set('COARSE_MAT',mld_distr_mat_,info) - case(mld_jac_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),mld_diag_scale_,info) @@ -594,6 +596,27 @@ contains info = -5 end if #endif +#ifdef HAVE_MUMPS_ + case (mld_mumps_) + if (allocated(level%sm%sv)) then + select type (sv => level%sm%sv) + class is (mld_z_mumps_solver_type) + ! do nothing + class default + call level%sm%sv%free(info) + if (info == 0) deallocate(level%sm%sv) + if (info == 0) allocate(mld_z_mumps_solver_type ::& + & level%sm%sv, stat=info) + end select + else + allocate(mld_z_mumps_solver_type :: level%sm%sv, stat=info) + endif + if (allocated(level%sm)) then + if (allocated(level%sm%sv)) & + & call level%sm%sv%default() + end if +#endif + #ifdef HAVE_UMF_ case (mld_umf_) if (allocated(level%sm)) then @@ -644,28 +667,6 @@ contains info = -5 end if #endif - -#ifdef HAVE_MUMPS_ - case (mld_mumps_) - if (allocated(level%sm%sv)) then - select type (sv => level%sm%sv) - class is (mld_z_mumps_solver_type) - ! do nothing - class default - call level%sm%sv%free(info) - if (info == 0) deallocate(level%sm%sv) - if (info == 0) allocate(mld_z_mumps_solver_type ::& - & level%sm%sv, stat=info) - end select - else - allocate(mld_z_mumps_solver_type :: level%sm%sv, stat=info) - endif - if (allocated(level%sm)) then - if (allocated(level%sm%sv)) & - & call level%sm%sv%default() - end if -#endif - case default ! ! Do nothing and hope for the best :) diff --git a/mlprec/impl/mld_zprecinit.F90 b/mlprec/impl/mld_zprecinit.F90 index 2b04946e..fa0f9337 100644 --- a/mlprec/impl/mld_zprecinit.F90 +++ b/mlprec/impl/mld_zprecinit.F90 @@ -2,9 +2,9 @@ !!$ !!$ MLD2P4 version 2.0 !!$ MultiLevel Domain Decomposition Parallel Preconditioners Package -!!$ based on PSBLAS (Parallel Sparse BLAS version 3.0) +!!$ based on PSBLAS (Parallel Sparse BLAS version 3.3) !!$ -!!$ (C) Copyright 2008,2009,2010,2012,2013 +!!$ (C) Copyright 2008, 2010, 2012, 2015 !!$ !!$ Salvatore Filippone University of Rome Tor Vergata !!$ Alfredo Buttari CNRS-IRIT, Toulouse @@ -61,7 +61,7 @@ ! ! 'ML' - Multilevel hybrid preconditioner (additive on the ! same level and multiplicative through the levels), -! with 2 levels and post-smoothing only. RAS with +! with 2 levels, pre and post-smoothing, RAS with ! overlap 1 and ILU(0) on the local blocks is ! applied as post-smoother at each level, but the ! coarsest one; four sweeps of the block-Jacobi solver, @@ -115,7 +115,7 @@ subroutine mld_zprecinit(p,ptype,info,nlev) ! Local variables integer(psb_ipk_) :: nlev_, ilev_ - real(psb_dpk_) :: thr + real(psb_dpk_) :: thr, scale character(len=*), parameter :: name='mld_precinit' info = psb_success_ @@ -194,7 +194,7 @@ subroutine mld_zprecinit(p,ptype,info,nlev) #if defined(HAVE_UMF_) allocate(mld_z_umf_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #elif defined(HAVE_SLU_) - allocate(mld_z_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) + allocate(mld_z_slu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #else allocate(mld_z_ilu_solver_type :: p%precv(ilev_)%sm%sv, stat=info) #endif @@ -205,10 +205,11 @@ subroutine mld_zprecinit(p,ptype,info,nlev) call p%precv(ilev_)%set(mld_sub_prol_,psb_none_,info) call p%precv(ilev_)%set(mld_sub_ovr_,izero,info) - thr = 0.16d0 + thr = 0.05 + scale = 1.0 do ilev_=1,nlev_ call p%precv(ilev_)%set(mld_aggr_thresh_,thr,info) - thr = thr/2 + call p%precv(ilev_)%set(mld_aggr_scale_,scale,info) end do case default diff --git a/mlprec/impl/mld_zprecset.F90 b/mlprec/impl/mld_zprecset.F90 index f0bdfd51..7b46a90c 100644 --- a/mlprec/impl/mld_zprecset.F90 +++ b/mlprec/impl/mld_zprecset.F90 @@ -98,6 +98,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev) #if defined(HAVE_MUMPS_) use mld_z_mumps_solver #endif + implicit none ! Arguments @@ -154,8 +155,8 @@ subroutine mld_zprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_) @@ -172,8 +173,8 @@ subroutine mld_zprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) case(mld_sub_solve_) call onelev_set_solver(p%precv(ilev_),val,info) - case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& - & mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& + case(mld_smoother_sweeps_,mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,& + & mld_aggr_kind_,mld_smoother_pos_,mld_aggr_omega_alg_,mld_aggr_eig_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_sub_restr_,mld_sub_prol_, & & mld_sub_ren_,mld_sub_ovr_,mld_sub_fillin_,& @@ -215,7 +216,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set(mld_coarse_mat_,mld_repl_mat_,info) - case(mld_sludist_, mld_mumps_) + case(mld_sludist_,mld_mumps_) call onelev_set_smoother(p%precv(nlev_),mld_bjac_,info) call onelev_set_solver(p%precv(nlev_),val,info) call p%precv(nlev_)%set(mld_coarse_mat_,mld_distr_mat_,info) @@ -284,7 +285,7 @@ subroutine mld_zprecseti(p,what,val,info,ilev) call onelev_set_smoother(p%precv(ilev_),val,info) end do - case(mld_ml_type_,mld_aggr_alg_,mld_aggr_kind_,& + case(mld_ml_type_,mld_aggr_alg_,mld_aggr_ord_,mld_aggr_kind_,& & mld_smoother_sweeps_pre_,mld_smoother_sweeps_post_,& & mld_smoother_pos_,mld_aggr_omega_alg_,& & mld_aggr_eig_,mld_aggr_filter_) @@ -308,8 +309,6 @@ subroutine mld_zprecseti(p,what,val,info,ilev) call onelev_set_solver(p%precv(nlev_),mld_umf_,info) #elif defined(HAVE_SLU_) call onelev_set_solver(p%precv(nlev_),mld_slu_,info) -#elif defined(HAVE_MUMPS_) - call onelev_set_solver(p%precv(nlev_),mld_mumps_,info) #else call onelev_set_solver(p%precv(nlev_),mld_ilu_n_,info) #endif @@ -782,13 +781,26 @@ subroutine mld_zprecsetsv(p,val,info,ilev) do ilev_ = ilmin, ilmax if (allocated(p%precv(ilev_)%sm)) then - if (allocated(p%precv(ilev_)%sm%sv)) & - & deallocate(p%precv(ilev_)%sm%sv) + if (allocated(p%precv(ilev_)%sm%sv)) then + if (.not.same_type_as(p%precv(ilev_)%sm%sv,val)) then + deallocate(p%precv(ilev_)%sm%sv,stat=info) + if (info /= 0) then + info = 3111 + return + end if + end if + if (.not.allocated(p%precv(ilev_)%sm%sv)) then #ifdef HAVE_MOLD - allocate(p%precv(ilev_)%sm%sv,mold=val) + allocate(p%precv(ilev_)%sm%sv,mold=val,stat=info) #else - allocate(p%precv(ilev_)%sm%sv,source=val) + allocate(p%precv(ilev_)%sm%sv,source=val,stat=info) #endif + if (info /= 0) then + info = 3111 + return + end if + end if + end if call p%precv(ilev_)%sm%sv%default() else info = 3111 diff --git a/mlprec/mld_base_prec_type.F90 b/mlprec/mld_base_prec_type.F90 index d3e07fbb..142c10a1 100644 --- a/mlprec/mld_base_prec_type.F90 +++ b/mlprec/mld_base_prec_type.F90 @@ -98,7 +98,7 @@ module mld_base_prec_type type mld_ml_parms integer(psb_ipk_) :: sweeps, sweeps_pre, sweeps_post integer(psb_ipk_) :: ml_type, smoother_pos - integer(psb_ipk_) :: aggr_alg, aggr_kind + integer(psb_ipk_) :: aggr_alg, aggr_ord, aggr_kind integer(psb_ipk_) :: aggr_omega_alg, aggr_eig, aggr_filter integer(psb_ipk_) :: coarse_mat, coarse_solve logical :: clean_zeros=.true. @@ -151,18 +151,19 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_smoother_pos_ = 23 integer(psb_ipk_), parameter :: mld_aggr_kind_ = 24 integer(psb_ipk_), parameter :: mld_aggr_alg_ = 25 - integer(psb_ipk_), parameter :: mld_aggr_omega_alg_ = 26 - integer(psb_ipk_), parameter :: mld_aggr_eig_ = 27 - integer(psb_ipk_), parameter :: mld_aggr_filter_ = 28 - integer(psb_ipk_), parameter :: mld_coarse_mat_ = 29 - integer(psb_ipk_), parameter :: mld_coarse_solve_ = 30 - integer(psb_ipk_), parameter :: mld_coarse_sweeps_ = 31 - integer(psb_ipk_), parameter :: mld_coarse_fillin_ = 32 - integer(psb_ipk_), parameter :: mld_coarse_subsolve_ = 33 - integer(psb_ipk_), parameter :: mld_smoother_sweeps_ = 34 - integer(psb_ipk_), parameter :: mld_coarse_aggr_size_ = 35 - integer(psb_ipk_), parameter :: mld_solver_sweeps_ = 36 - integer(psb_ipk_), parameter :: mld_ifpsz_ = 37 + integer(psb_ipk_), parameter :: mld_aggr_ord_ = 26 + integer(psb_ipk_), parameter :: mld_aggr_omega_alg_ = 27 + integer(psb_ipk_), parameter :: mld_aggr_eig_ = 28 + integer(psb_ipk_), parameter :: mld_aggr_filter_ = 29 + integer(psb_ipk_), parameter :: mld_coarse_mat_ = 30 + integer(psb_ipk_), parameter :: mld_coarse_solve_ = 31 + integer(psb_ipk_), parameter :: mld_coarse_sweeps_ = 32 + integer(psb_ipk_), parameter :: mld_coarse_fillin_ = 33 + integer(psb_ipk_), parameter :: mld_coarse_subsolve_ = 34 + integer(psb_ipk_), parameter :: mld_smoother_sweeps_ = 35 + integer(psb_ipk_), parameter :: mld_coarse_aggr_size_ = 36 + integer(psb_ipk_), parameter :: mld_solver_sweeps_ = 37 + integer(psb_ipk_), parameter :: mld_ifpsz_ = 38 ! ! Legal values for entry: mld_smoother_type_ @@ -249,8 +250,13 @@ module mld_base_prec_type integer(psb_ipk_), parameter :: mld_glb_aggr_=2 integer(psb_ipk_), parameter :: mld_new_dec_aggr_=3 integer(psb_ipk_), parameter :: mld_new_glb_aggr_=4 - integer(psb_ipk_), parameter :: mld_max_aggr_alg_=mld_dec_aggr_ - + integer(psb_ipk_), parameter :: mld_max_aggr_alg_=mld_sym_dec_aggr_ + ! + ! Legal values for entry: mld_aggr_ord_ + ! + integer(psb_ipk_), parameter :: mld_aggr_ord_nat_ = 0 + integer(psb_ipk_), parameter :: mld_aggr_ord_desc_deg_ = 1 + integer(psb_ipk_), parameter :: mld_max_aggr_ord_ = mld_aggr_ord_desc_deg_ ! ! Legal values for entry: mld_aggr_omega_alg_ ! @@ -320,6 +326,8 @@ module mld_base_prec_type character(len=18), parameter, private :: & & aggr_names(0:4)=(/'local aggregation ','sym. local aggr. ',& & 'global aggregation', 'new local aggr. ','new global aggr. '/) + character(len=18), parameter, private :: & + & ord_names(0:1)=(/'Natural ordering ','Desc. degree ord. '/) character(len=6), parameter, private :: & & restrict_names(0:4)=(/'none ','halo ',' ',' ',' '/) character(len=12), parameter, private :: & @@ -411,6 +419,10 @@ contains val = mld_dec_aggr_ case('SYMDEC') val = mld_sym_dec_aggr_ + case('NAT','NATURAL') + val = mld_aggr_ord_nat_ + case('DESC','RDEGREE','DEGREE') + val = mld_aggr_ord_desc_deg_ case('GLB') val = mld_glb_aggr_ case('REPL') @@ -467,7 +479,7 @@ contains write(iout,*) 'Sweeps: ',pm%sweeps,pm%sweeps_pre,pm%sweeps_post write(iout,*) 'ML : ',pm%ml_type,pm%smoother_pos - write(iout,*) 'AGGR : ',pm%aggr_alg,pm%aggr_kind + write(iout,*) 'AGGR : ',pm%aggr_alg,pm%aggr_kind, pm%aggr_ord write(iout,*) ' : ',pm%aggr_omega_alg,pm%aggr_eig,pm%aggr_filter write(iout,*) 'COARSE: ',pm%coarse_mat,pm%coarse_solve end subroutine ml_parms_printout @@ -533,6 +545,8 @@ contains end if write(iout,*) ' Aggregation: ', & & aggr_names(pm%aggr_alg) + write(iout,*) ' with initial ordering: ',& + & ord_names(pm%aggr_ord) write(iout,*) ' Aggregation type: ', & & aggr_kinds(pm%aggr_kind) if (pm%aggr_kind /= mld_no_smooth_) then @@ -720,6 +734,14 @@ contains is_legal_ml_aggr_alg = ((ip>=mld_dec_aggr_).and.(ip<=mld_max_aggr_alg_)) return end function is_legal_ml_aggr_alg + function is_legal_ml_aggr_ord(ip) + implicit none + integer(psb_ipk_), intent(in) :: ip + logical :: is_legal_ml_aggr_ord + + is_legal_ml_aggr_ord = ((mld_aggr_ord_nat_<=ip).and.(ip<=mld_max_aggr_ord_)) + return + end function is_legal_ml_aggr_ord function is_legal_ml_aggr_omega_alg(ip) implicit none integer(psb_ipk_), intent(in) :: ip diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index f0840043..a1aa0253 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -109,9 +109,10 @@ module mld_c_inner_mod end interface mld_coarse_bld interface mld_aggrmap_bld - subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) + subroutine mld_caggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_spk_), intent(in) :: theta type(psb_cspmat_type), intent(in) :: a @@ -123,9 +124,10 @@ module mld_c_inner_mod interface mld_dec_map_bld - subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta diff --git a/mlprec/mld_c_onelev_mod.f90 b/mlprec/mld_c_onelev_mod.f90 index ddc0266a..e0115610 100644 --- a/mlprec/mld_c_onelev_mod.f90 +++ b/mlprec/mld_c_onelev_mod.f90 @@ -379,6 +379,7 @@ contains lv%parms%sweeps_post = 1 lv%parms%ml_type = mld_mult_ml_ lv%parms%aggr_alg = mld_dec_aggr_ + lv%parms%aggr_ord = mld_aggr_ord_nat_ lv%parms%aggr_kind = mld_smooth_prol_ lv%parms%coarse_mat = mld_distr_mat_ lv%parms%smoother_pos = mld_twoside_smooth_ diff --git a/mlprec/mld_d_base_solver_mod.f90 b/mlprec/mld_d_base_solver_mod.f90 index 6971fed1..89f6cd3c 100644 --- a/mlprec/mld_d_base_solver_mod.f90 +++ b/mlprec/mld_d_base_solver_mod.f90 @@ -360,6 +360,7 @@ contains integer(psb_long_int_k_) :: val integer(psb_ipk_) :: i val = 0 + return end function d_base_solver_sizeof diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 3c6b8bea..84bf1fc5 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -109,9 +109,10 @@ module mld_d_inner_mod end interface mld_coarse_bld interface mld_aggrmap_bld - subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) + subroutine mld_daggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_dpk_), intent(in) :: theta type(psb_dspmat_type), intent(in) :: a @@ -123,9 +124,10 @@ module mld_d_inner_mod interface mld_dec_map_bld - subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta diff --git a/mlprec/mld_d_onelev_mod.f90 b/mlprec/mld_d_onelev_mod.f90 index 3401115b..ee3b1b7f 100644 --- a/mlprec/mld_d_onelev_mod.f90 +++ b/mlprec/mld_d_onelev_mod.f90 @@ -338,6 +338,7 @@ contains class(mld_d_onelev_type), intent(in) :: lv integer(psb_long_int_k_) :: val integer(psb_ipk_) :: i + val = 0 val = val + lv%desc_ac%sizeof() val = val + lv%ac%sizeof() @@ -378,6 +379,7 @@ contains lv%parms%sweeps_post = 1 lv%parms%ml_type = mld_mult_ml_ lv%parms%aggr_alg = mld_dec_aggr_ + lv%parms%aggr_ord = mld_aggr_ord_nat_ lv%parms%aggr_kind = mld_smooth_prol_ lv%parms%coarse_mat = mld_distr_mat_ lv%parms%smoother_pos = mld_twoside_smooth_ @@ -386,7 +388,7 @@ contains lv%parms%aggr_filter = mld_no_filter_mat_ lv%parms%aggr_omega_val = dzero lv%parms%aggr_thresh = dzero - + if (allocated(lv%sm)) call lv%sm%default() return @@ -402,8 +404,8 @@ contains ! Arguments class(mld_d_onelev_type), target, intent(inout) :: lv class(mld_d_onelev_type), intent(inout) :: lvout - integer(psb_ipk_), intent(out) :: info - + integer(psb_ipk_), intent(out) :: info + info = psb_success_ if (allocated(lv%sm)) then call lv%sm%clone(lvout%sm,info) @@ -430,7 +432,7 @@ contains implicit none type(mld_d_onelev_type), intent(inout) :: a, b integer(psb_ipk_), intent(out) :: info - + call b%free(info) b%parms = a%parms call move_alloc(a%sm,b%sm) diff --git a/mlprec/mld_d_prec_type.f90 b/mlprec/mld_d_prec_type.f90 index 765f6999..134baa15 100644 --- a/mlprec/mld_d_prec_type.f90 +++ b/mlprec/mld_d_prec_type.f90 @@ -262,7 +262,7 @@ contains integer(psb_ipk_), optional :: ilev class(mld_d_base_smoother_type), pointer :: val integer(psb_ipk_) :: ilev_ - + val => null() if (present(ilev)) then ilev_ = ilev @@ -327,7 +327,6 @@ contains class(mld_dprec_type), intent(in) :: prec integer(psb_long_int_k_) :: val integer(psb_ipk_) :: i - val = 0 val = val + psb_sizeof_int if (allocated(prec%precv)) then @@ -530,6 +529,7 @@ contains ! Local variables integer(psb_ipk_) :: me,err_act,i character(len=20) :: name + if(psb_get_errstatus().ne.0) return info=psb_success_ name = 'mld_dprecfree' @@ -799,7 +799,7 @@ contains type(mld_dprec_type), intent(inout), target :: b integer(psb_ipk_), intent(out) :: info integer(psb_ipk_) :: i - + if (allocated(b%precv)) then ! This might not be required if FINAL procedures are available. call mld_precfree(b,info) diff --git a/mlprec/mld_d_umf_solver.F90 b/mlprec/mld_d_umf_solver.F90 index 0f7d76da..e76afb51 100644 --- a/mlprec/mld_d_umf_solver.F90 +++ b/mlprec/mld_d_umf_solver.F90 @@ -133,6 +133,7 @@ contains integer :: ictxt,np,me,i, err_act character :: trans_ character(len=20) :: name='d_umf_solver_apply' + call psb_erractionsave(err_act) info = psb_success_ @@ -216,7 +217,6 @@ contains integer :: err_act character(len=20) :: name='d_umf_solver_apply_vect' - call psb_erractionsave(err_act) call psb_erractionsave(err_act) info = psb_success_ @@ -257,8 +257,6 @@ contains integer :: n_row,n_col, nrow_a, nztota integer :: ictxt,np,me,i, err_act, debug_unit, debug_level character(len=20) :: name='d_umf_solver_bld', ch_err - - call psb_erractionsave(err_act) info=psb_success_ call psb_erractionsave(err_act) diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index e79453bc..74cbcdfb 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -109,9 +109,10 @@ module mld_s_inner_mod end interface mld_coarse_bld interface mld_aggrmap_bld - subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) + subroutine mld_saggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_spk_), intent(in) :: theta type(psb_sspmat_type), intent(in) :: a @@ -123,9 +124,10 @@ module mld_s_inner_mod interface mld_dec_map_bld - subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_spk_), intent(in) :: theta diff --git a/mlprec/mld_s_onelev_mod.f90 b/mlprec/mld_s_onelev_mod.f90 index f8e89935..7122f167 100644 --- a/mlprec/mld_s_onelev_mod.f90 +++ b/mlprec/mld_s_onelev_mod.f90 @@ -379,6 +379,7 @@ contains lv%parms%sweeps_post = 1 lv%parms%ml_type = mld_mult_ml_ lv%parms%aggr_alg = mld_dec_aggr_ + lv%parms%aggr_ord = mld_aggr_ord_nat_ lv%parms%aggr_kind = mld_smooth_prol_ lv%parms%coarse_mat = mld_distr_mat_ lv%parms%smoother_pos = mld_twoside_smooth_ diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index d14123ec..e9c3d360 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -109,9 +109,10 @@ module mld_z_inner_mod end interface mld_coarse_bld interface mld_aggrmap_bld - subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) + subroutine mld_zaggrmap_bld(aggr_type,iorder,theta,a,desc_a,ilaggr,nlaggr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder integer(psb_ipk_), intent(in) :: aggr_type real(psb_dpk_), intent(in) :: theta type(psb_zspmat_type), intent(in) :: a @@ -123,9 +124,10 @@ module mld_z_inner_mod interface mld_dec_map_bld - subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_, psb_ipk_ implicit none + integer(psb_ipk_), intent(in) :: iorder type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a real(psb_dpk_), intent(in) :: theta diff --git a/mlprec/mld_z_onelev_mod.f90 b/mlprec/mld_z_onelev_mod.f90 index 4f60dccb..e94a7d5c 100644 --- a/mlprec/mld_z_onelev_mod.f90 +++ b/mlprec/mld_z_onelev_mod.f90 @@ -379,6 +379,7 @@ contains lv%parms%sweeps_post = 1 lv%parms%ml_type = mld_mult_ml_ lv%parms%aggr_alg = mld_dec_aggr_ + lv%parms%aggr_ord = mld_aggr_ord_nat_ lv%parms%aggr_kind = mld_smooth_prol_ lv%parms%coarse_mat = mld_distr_mat_ lv%parms%smoother_pos = mld_twoside_smooth_ diff --git a/tests/pdegen/ppde2d.f90 b/tests/pdegen/ppde2d.f90 index 907a29e8..c96d4eae 100644 --- a/tests/pdegen/ppde2d.f90 +++ b/tests/pdegen/ppde2d.f90 @@ -159,6 +159,7 @@ program ppde2d integer(psb_ipk_) :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: aggr_ord ! Ordering for aggregation character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing integer(psb_ipk_) :: csize ! aggregation size at which to stop. @@ -241,6 +242,7 @@ program ppde2d call mld_precset(prec,'sub_iluthrs', prectype%thr1, info) call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) + call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) call mld_precset(prec,'ml_type', prectype%mltype, info) call mld_precset(prec,'smoother_pos', prectype%smthpos, info) if (prectype%athres >= dzero) & @@ -383,6 +385,7 @@ contains call read_data(prectype%nlev,psb_inp_unit) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,psb_inp_unit) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,psb_inp_unit) ! local or global aggregation + call read_data(prectype%aggr_ord,psb_inp_unit) ! aggregation ordering call read_data(prectype%mltype,psb_inp_unit) ! additive or multiplicative 2nd level prec call read_data(prectype%smthpos,psb_inp_unit) ! side: pre, post, both smoothing call read_data(prectype%cmat,psb_inp_unit) ! coarse mat @@ -422,6 +425,7 @@ contains call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%aggr_ord) ! aggregation ordering call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat diff --git a/tests/pdegen/ppde3d.f90 b/tests/pdegen/ppde3d.f90 index b8730ce2..47a28a2e 100644 --- a/tests/pdegen/ppde3d.f90 +++ b/tests/pdegen/ppde3d.f90 @@ -171,6 +171,7 @@ program ppde3d integer(psb_ipk_) :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: aggr_ord ! Ordering for aggregation character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing integer(psb_ipk_) :: csize ! aggregation size at which to stop. @@ -253,6 +254,7 @@ program ppde3d call mld_precset(prec,'sub_iluthrs', prectype%thr1, info) call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) + call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) call mld_precset(prec,'ml_type', prectype%mltype, info) call mld_precset(prec,'smoother_pos', prectype%smthpos, info) if (prectype%athres >= dzero) & @@ -396,6 +398,7 @@ contains call read_data(prectype%nlev,psb_inp_unit) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,psb_inp_unit) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,psb_inp_unit) ! local or global aggregation + call read_data(prectype%aggr_ord,psb_inp_unit) ! aggregation ordering call read_data(prectype%mltype,psb_inp_unit) ! additive or multiplicative 2nd level prec call read_data(prectype%smthpos,psb_inp_unit) ! side: pre, post, both smoothing call read_data(prectype%cmat,psb_inp_unit) ! coarse mat @@ -435,6 +438,7 @@ contains call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%aggr_ord) ! aggregation ordering call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat @@ -482,4 +486,3 @@ contains end subroutine pr_usage end program ppde3d - diff --git a/tests/pdegen/runs/ppde.inp b/tests/pdegen/runs/ppde.inp index 7e9993c4..c77e56c6 100644 --- a/tests/pdegen/runs/ppde.inp +++ b/tests/pdegen/runs/ppde.inp @@ -3,15 +3,15 @@ CSR ! Storage format CSR COO JAD 0100 ! IDIM; domain size is idim**3 2 ! ISTOPC 0100 ! ITMAX -1 ! ITRACE +10 ! ITRACE 30 ! IRST (restart for RGMRES and BiCGSTABL) 1.d-6 ! EPS 3L-MUL-RAS-BJAC4-ILU ! Descriptive name for preconditioner (up to 40 chars) -ML ! Preconditioner NONE JACOBI BJAC AS ML +ML ! Preconditioner NONE JACOBI BJAC AS ML 1 ! Number of overlap layers for AS preconditioner at finest level HALO ! Restriction operator NONE HALO NONE ! Prolongation operator NONE SUM AVG -GS ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU +GS ! Subdomain solver DSCALE ILU MILU ILUT UMF SLU 1 ! sweeps for GS 0 ! Level-set N for ILU(N), and P for ILUT 1.d-4 ! Threshold T for ILU(T,P) @@ -20,11 +20,12 @@ BJAC ! Smoother type JACOBI BJAC AS; ignored for non-ML 2 ! Number of levels in a multilevel preconditioner SMOOTHED ! Kind of aggregation: SMOOTHED, NONSMOOTHED, MINENERGY DEC ! Type of aggregation DEC SYMDEC GLB +DEGREE ! Ordering of aggregation NATURAL DEGREE MULT ! Type of multilevel correction: ADD MULT TWOSIDE ! Side of correction PRE POST TWOSIDE (ignored for ADD) REPL ! Coarse level: matrix distribution DIST REPL BJAC ! Coarse level: solver JACOBI BJAC UMF SLU SLUDIST MUMPS -MUMPS ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST MUMPS +MUMPS ! Coarse level: subsolver DSCALE ILU UMF SLU SLUDIST MUMPS 1 ! Coarse level: Level-set N for ILU(N) 1.d-4 ! Coarse level: Threshold T for ILU(T,P) 4 ! Coarse level: Number of Jacobi sweeps diff --git a/tests/pdegen/spde2d.f90 b/tests/pdegen/spde2d.f90 index b6decae2..ccff4dd4 100644 --- a/tests/pdegen/spde2d.f90 +++ b/tests/pdegen/spde2d.f90 @@ -159,6 +159,7 @@ program spde2d integer(psb_ipk_) :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: aggr_ord ! Ordering for aggregation character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing integer(psb_ipk_) :: csize ! aggregation size at which to stop. @@ -241,6 +242,7 @@ program spde2d call mld_precset(prec,'sub_iluthrs', prectype%thr1, info) call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) + call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) call mld_precset(prec,'ml_type', prectype%mltype, info) call mld_precset(prec,'smoother_pos', prectype%smthpos, info) if (prectype%athres >= dzero) & @@ -384,6 +386,7 @@ contains call read_data(prectype%nlev,psb_inp_unit) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,psb_inp_unit) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,psb_inp_unit) ! local or global aggregation + call read_data(prectype%aggr_ord,psb_inp_unit) ! aggregation ordering call read_data(prectype%mltype,psb_inp_unit) ! additive or multiplicative 2nd level prec call read_data(prectype%smthpos,psb_inp_unit) ! side: pre, post, both smoothing call read_data(prectype%cmat,psb_inp_unit) ! coarse mat @@ -423,6 +426,7 @@ contains call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%aggr_ord) ! aggregation ordering call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat diff --git a/tests/pdegen/spde3d.f90 b/tests/pdegen/spde3d.f90 index ce397f94..308cd03b 100644 --- a/tests/pdegen/spde3d.f90 +++ b/tests/pdegen/spde3d.f90 @@ -171,6 +171,7 @@ program spde3d integer(psb_ipk_) :: nlev ! Number of levels in multilevel prec. character(len=16) :: aggrkind ! smoothed/raw aggregatin character(len=16) :: aggr_alg ! local or global aggregation + character(len=16) :: aggr_ord ! Ordering for aggregation character(len=16) :: mltype ! additive or multiplicative 2nd level prec character(len=16) :: smthpos ! side: pre, post, both smoothing integer(psb_ipk_) :: csize ! aggregation size at which to stop. @@ -254,6 +255,7 @@ program spde3d call mld_precset(prec,'sub_iluthrs', prectype%thr1, info) call mld_precset(prec,'aggr_kind', prectype%aggrkind,info) call mld_precset(prec,'aggr_alg', prectype%aggr_alg,info) + call mld_precset(prec,'aggr_ord', prectype%aggr_ord,info) call mld_precset(prec,'ml_type', prectype%mltype, info) call mld_precset(prec,'smoother_pos', prectype%smthpos, info) if (prectype%athres >= dzero) & @@ -397,6 +399,7 @@ contains call read_data(prectype%nlev,psb_inp_unit) ! Number of levels in multilevel prec. call read_data(prectype%aggrkind,psb_inp_unit) ! smoothed/raw aggregatin call read_data(prectype%aggr_alg,psb_inp_unit) ! local or global aggregation + call read_data(prectype%aggr_ord,psb_inp_unit) ! aggregation ordering call read_data(prectype%mltype,psb_inp_unit) ! additive or multiplicative 2nd level prec call read_data(prectype%smthpos,psb_inp_unit) ! side: pre, post, both smoothing call read_data(prectype%cmat,psb_inp_unit) ! coarse mat @@ -436,6 +439,7 @@ contains call psb_bcast(ictxt,prectype%nlev) ! Number of levels in multilevel prec. call psb_bcast(ictxt,prectype%aggrkind) ! smoothed/raw aggregatin call psb_bcast(ictxt,prectype%aggr_alg) ! local or global aggregation + call psb_bcast(ictxt,prectype%aggr_ord) ! aggregation ordering call psb_bcast(ictxt,prectype%mltype) ! additive or multiplicative 2nd level prec call psb_bcast(ictxt,prectype%smthpos) ! side: pre, post, both smoothing call psb_bcast(ictxt,prectype%cmat) ! coarse mat From c592733a8ef53437de4a925b36bc10e7dc452ea9 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Mon, 9 May 2016 16:04:33 +0000 Subject: [PATCH 2/2] mld2p4-2: docs/html/img100.png docs/html/img101.png docs/html/img102.png docs/html/img103.png docs/html/img104.png docs/html/img105.png docs/html/img92.png docs/html/img93.png docs/html/img95.png docs/html/img96.png docs/html/img97.png docs/html/img98.png docs/html/img99.png docs/html/index.html docs/html/node1.html docs/html/node10.html docs/html/node11.html docs/html/node12.html docs/html/node13.html docs/html/node14.html docs/html/node15.html docs/html/node16.html docs/html/node17.html docs/html/node18.html docs/html/node19.html docs/html/node2.html docs/html/node20.html docs/html/node21.html docs/html/node22.html docs/html/node23.html docs/html/node24.html docs/html/node25.html docs/html/node26.html docs/html/node3.html docs/html/node4.html docs/html/node5.html docs/html/node6.html docs/html/node7.html docs/html/node8.html docs/html/node9.html docs/html/userhtml.html docs/mld2p4-2.1-guide.pdf Added aggregation ordering option to docs. --- docs/html/img100.png | Bin 600 -> 611 bytes docs/html/img101.png | Bin 251 -> 245 bytes docs/html/img102.png | Bin 193 -> 195 bytes docs/html/img103.png | Bin 210 -> 210 bytes docs/html/img104.png | Bin 610 -> 614 bytes docs/html/img105.png | Bin 650 -> 645 bytes docs/html/img92.png | Bin 201 -> 201 bytes docs/html/img93.png | Bin 328 -> 328 bytes docs/html/img95.png | Bin 293 -> 293 bytes docs/html/img96.png | Bin 337 -> 340 bytes docs/html/img97.png | Bin 667 -> 676 bytes docs/html/img98.png | Bin 467 -> 471 bytes docs/html/img99.png | Bin 250 -> 250 bytes docs/html/index.html | 81 +- docs/html/node1.html | 44 +- docs/html/node10.html | 23 +- docs/html/node11.html | 76 +- docs/html/node12.html | 46 +- docs/html/node13.html | 40 +- docs/html/node14.html | 38 +- docs/html/node15.html | 100 +- docs/html/node16.html | 80 +- docs/html/node17.html | 50 +- docs/html/node18.html | 188 +-- docs/html/node19.html | 57 +- docs/html/node2.html | 67 +- docs/html/node20.html | 112 +- docs/html/node21.html | 29 +- docs/html/node22.html | 29 +- docs/html/node23.html | 25 +- docs/html/node24.html | 23 +- docs/html/node25.html | 59 +- docs/html/node26.html | 14 +- docs/html/node3.html | 83 +- docs/html/node4.html | 25 +- docs/html/node5.html | 35 +- docs/html/node6.html | 63 +- docs/html/node7.html | 53 +- docs/html/node8.html | 53 +- docs/html/node9.html | 23 +- docs/html/userhtml.html | 81 +- docs/mld2p4-2.1-guide.pdf | 2635 ++++++++++++++++++++----------------- 42 files changed, 2219 insertions(+), 2013 deletions(-) diff --git a/docs/html/img100.png b/docs/html/img100.png index f8eef1828f38f5a8ca528fef4c93f5503a395ed3..bda1e4801baa4c2b4654a79dfb57197868a620c2 100644 GIT binary patch delta 547 zcmV+;0^I%B1mgrDiBL{Q4GJ0x0000DNk~Le0001X0000d1Oos704Q6baFHQQf8j|) zK~zYI?UO%D+dve?KmUnw4Zb*4(*YQ~L8Vl9K*~U(rW@*1m4Vd*AONK{=Qz`_RajA$=$28(W0KFPj&_ulWlcfK?5-^;l_ zx*8-7H5%P30aA7{b{;s~piw!ufBvxY%Vcboo^=;^68^1=Vj`4L38T;g$|dNTOy+UK zs)?FT5oV5gU~_jZb5`<%kk<`tbuf$jP=?13Fi3MLZDBm;2rTag+s$H>!g_%#8cw-R zXrjL|t*5^lsyoM7sRE~P*Jk;v2ctwg@g|S?H_#VC)I~OLv*1%6>Bbouf16U<@E>)w zMYT=j4PIB1D(sEudj}S{Sqxa5Z7P9|=iRJ^?Co6jB=OAH512JfW69yH#ljW&TQ>{x zM;fxzy#B^D>5|w0swLO4haJx1yRU>ID%bd{jutm^SaL;Qh>x}s$nP;ki5z@}R1#U% z2B!3YMQmz!wxg%=;X_DkdB0nD;AzTD?OX{X8>0(2Y9aGg=>0 zay-gdVHKgy?=jVeT6zn5tR>p-*c+&E8%I1FHi@+qhJO7*gKgnOOIGKMe`V2|AgnQ8 zMOY=1d~OEBiUm;K$k?=Gsp!dP6$G+-IU(@>jGVx?hhB_#G1egVym*+;g66XiVCu|& z>yzr8vKC}RJi@u=usBXiHjj$z328pVf*{6Ye9p$ADOutyGaw=CLx}P$p)EU9i#TlA z*-*|*b2n8G)}AuhfKuO!Zgo;W#rO85fbPiPLNptH1)`^XcHN1_*M>l#%VPnd;2u)dIX?m#``iFK(DtS2a(I-qypz=86v#&YIn zcIQ6E3AGK4t_?3G4wo(vTIS57%XNl*#w&+sCUzo833GmYV(d7}q7pev?7-T~o(;w~ lQ{)tGyqCLpqJ9d{ExA(D+xRC&L;zjG;OXk;vd$@?2>_xAMfU&z diff --git a/docs/html/img102.png b/docs/html/img102.png index 06cc29ded3b2065ecb7879d40427a96fc66ea77f..891cce5fd222dfa34c28bcce6df915852e017342 100644 GIT binary patch delta 131 zcmX@ec$iVKGr-TCmrII^fq{Y7)59eQNb>`+5DPPqEZZ>CW}>2Jq>HDEV+hC0^QzHJR-Fyu23Kf}=)5DPPqoLTgVX`-TMxRa-gV+hC0BcRT9-M8%B&uqO!} delta 41 scmcb_c!^Q5Gr-TCmrII^fq{Y7)59eQNb>=)JPR|BRDJXC&P2tH0Irq^3IG5A diff --git a/docs/html/img104.png b/docs/html/img104.png index c34f3c6d3bf96022690bb6c75df6becf390ca3bf..b840012ee9ca9c2a6adaa95448654cc7b67ecf96 100644 GIT binary patch delta 515 zcmV+e0{s2r1m*;gZGYrRL_t(YiS3imD+FO2#-Ev4Yj)XpXT@d_bF*sI+*}sxw#UsM zpi#)d#m$A*CI_Y1qm(T-r5wn0tdb;cdm%}?gp{=Y0^WBr>^mmA#YOh1cYb`|XFlH_ zy#xFUNfjMhg?;&|5@Fe}1-9*uMxajpmSID<6F_l;K1eUh2Y*6L@Fo%}(d%Hq8R20G z*L9#L`&kA0ZFxOMl)G-#W<%W3<&74A`55FtAP_VWE?Xy?XLLx-)I=S+9Wg?euR94c zrj%77BzaeU2Dl%x-rfJ&Q1WJ)xwHm{@M>My5@LL1 zvml1vg8EK|$V-Ndp)ZI}`P3+qE`^Jc6Zm!z$J78rC?uN_4><#m5bzXMBmB8023(4G zgmbTnz;WiIQ&S?FvgegaGZBqTu!dz>vRF_hL!4vjh50{pP(REtdz zRv!ykgS$PBG-;pK%bU-Yn@4w(9<{h1G+G{&sMZUv~_DQH1&Yj$7l^5!@ZL z5=c3xy>yj`=Nooga@dkH)gdfWm9d&?s_75v^vC)I#2d_%+53C$^_kSU|0 z3d!bOmkYrCi2GJT%v@n$Cz-+l`=HrCZ`i6p5Mda`GvwtgkI=^SmYYS4s!@iJ)?h_S zwQHU2v^5rfWPiHxhLT@OYg%#g%o7#yI$4(Ww?&5bD|j zHlTLIuueMW6P3+hRGKFaZNFJOZsQSh0f_a}LlRzUPI!hz$vMaDnVYVOfQdZYl?$)j zdTof=N0U8PJyKe#HbhBkjj_7ws_PHx{KxtQ#5=4K*xujKt+N0C002ovPDHLkV1h=L B?XUm< diff --git a/docs/html/img105.png b/docs/html/img105.png index a6135e0eb306ddbb10b55c0bc383726d680ff202..461a8c822b1240e306263bf1c9b5cf0313119e90 100644 GIT binary patch delta 547 zcmV+;0^I$I1%(BWZGQtvL_t(YiS3fVOB-Pv#-F>3$0g=o??gg_lHO)28LpE{q;V-3 zL;rvhsfA9Wn~QLYlOTv=vGvfcMQ~^-IK>2llp?{Yg|;L}ol1x}m-c-xcP8(f5L0xL zJn-K8>v=!l9~>~y|4F&Y@Z&sDtM;~N-)~^3hQ6iP{MmtMihlw_R<9@sR|ubpngeRP zbfH`t*AO|!za#Mm{XQ|9FEH%N{C$`O-N^%wimgmABSa6EDS0Qz@4cAZFFGP=n%!y! znrWIb4-yjU**~RajQIw96Cuy(&?g;`d7~{8$>CjR%fQ8KuzM3MUVX*my-gQ~?8%6M zWXf)X5TPl?C4c0#Qm{MnXYH9j`f~+O)eM73Td>o91=_>|S|B^FGhH+xaVSdAH1pVx zH*JxV7mN5qSljodM#LTCLKiql#D0i(QbugW3P@SXVabD>DFqF)9gvoED-9=PfEW1; z{V?RHB@bed;+kCI3_L`_@1&Iy{!)+vLPjdDFGf5FPJi6%iap7!-@AXQu9}eO27;!w zOhv5I6ep`JB(d=ke8paZEuCXH2fk@ctF=+Y6K{dnP{*^_CgV%__9j%?&2t}|uvPmM z6CrXH$WGD4W(S{XB;WP8Q*zm}cfEVrCj)b*?5k9|-EH>^S*TlV+URnnx&1%wt002ovPDHLkV1oG-2Q>fy delta 552 zcmV+@0@wY81&RfbZGQ+!L_t(YiS3e2OB+EP#-H6u)7|J!W{XOF!S2lzQrMf9SmVuO z=m#h(w8cX&y?Lp`P&^0)>#(Dm9Od^?P40k?RE}pi zmxp27c3%idDD|v<;AQl16L)y{tU~VcVT}|)@3EFlw4&_VodJGKM1L=$#q=kTfy@`I ziz{h{OxA6IP=BQ<>==4!JNnyf%*Gvy-8!r!9^?f?TZEfyRA3BU=LL#wdL-n9kbV#n zHO<~?#2a;JFh7d)H_`T0t)B9Gq|ku`Nx3($rsdUEJO&NCUR)g(kX;fOm^cS%%5rij zCKzPVcd*tmiHElZgwb+I-gd1I`@dK4Par4mR^T}Gmr#lT0000>+rZ5=QE>$Ts8|QI diff --git a/docs/html/img93.png b/docs/html/img93.png index 80faaac9be44b753f85286c1e72cc9b6db17c623..8ef15cb236fd166f1af443fcfa59fb40c4e2c16b 100644 GIT binary patch delta 217 zcmV;~04D#)0>}c8cYjd8{R5fD_#Z`*pMgQw7|J#Pb5SJOoHznNEI!Et2+4Kt9)KlT z8yG%=1^5t>_QFs}o(Bw)Xp#`4K)g2#7{t*eHQ{!?VPKF!le~;iQsn>t{~sWd3K)`d z2%QfYM9?Ic!X>#MFodE>0(}RSWEEhjW3W#^@n_%ys3g-rMp6L=n+ptV@+ZY$l3emo z2QlzH1JV~@1|b6|k8J@SNk#^2k}S`_loHgg_b8Hl&yaasuNX!FEx{Q)0IwtfvsGOn Tr9&NR00000NkvXXu0mjfWr0<* delta 217 zcmV;~04D#)0>}c8cYj#0zX6#i{{cnv1B0$Hl*tU{qDZniaRh)_59c5x*S&iHmSk;U z_zV$w1DCWHhD!20V30(Ugc=3py;;B@jwY!Iw-YENgC=W@E`>{SKVS$&lLYz>D#HB0;K|wlYgt%fbA~BLa>N{7XwcMgD8U-Bg+CCkN`*8A%^D+yafzGZ8-`o zAQ3jk2|xwhAMY+mkYr#4>tZ)x_`$H|0g&avzyTI{&CGOc(Haoz5Q7<5M39m1zybdJ}OVl7iML)Bpeg07*qo IM6N<$f^>B|d;kCd delta 154 zcmV;L0A>HB0;K|wlYe2wWgv@703;&d#lVvQRK}Pf;Kh6ZD8P|+h~YT{!*vEm4Is+_ zD8j}#fkBjkVZZDx5DTn}-GJc-1J^kqi=}~~0VML8nd#UfE|38{3d{*05kW@20~Z*! zF)%OSVc-DUZra7d^2<^?zy+ip0dW>;K;Z@t2q_o^Llm6B4-ODi005iVXM0PpiVOe% N002ovPDHLkV1h{oQ05AX40?`7HZ2`}bZXkb(3?r0&gAfV8ry%zN1AhYpSFZuvU516&6bN`R@FXyZ zGKev|%fFtb?!*d4S0tTVB90e9E3fLGY0JU>}yt^Pl5-c_U|M&m@KO|%Y%Ufn;Tv;GaXyB2E;nVV1`YBAS2&_3)~3|RtyIi8n7ua?P5f7zQ8PJyMQ6^ ztpc)w#6J@gz^+ta-~t9FP>yW@YCz!xe1sH?f}sk|;0Ff?DgXd+$7M2YP?T~20000< KMNUMnLSTXsFHrIT diff --git a/docs/html/img97.png b/docs/html/img97.png index fa31ee53afb5bfc29ed2a6e4bf4f3023b41e8cec..faee56df4ca9b196bee5ffaeb626e0e19688483e 100644 GIT binary patch delta 578 zcmV-I0=@m41*8R#ZGR$3L_t(YiS3imYZE~f$G>TsZ8pjNs%@2e*^4MxAvZl3(4Is* z1VN9o1-AS_hHgJS*xVFM~CqVX(3EEIZJ(p6OSq~F<@tdmL*axD46 z3^Q-$^WK}8H^9DBeOz!o9=SB-0#CJhr8~SXNgKEbl>CA9F@Kv?htY)Lh01mf`25S8 zpB*L>;0U8M2iCJ}4q*CdU{3kF!}c8cRFM~eB9mz-qb0}{U>o>`L}qEzZRHIx?Up^y zTTrBRy9X@Os{Cz!VYdie0X`s+>6W;~Y6FEcddw3E&4*}JXL1FZ0&Ou7=^ZRQrNFC)0%pPB!75@>^ZGvJy{x z4t9F@#8mzkwIjN9rqPcu9r9MndkS9D4{i_g*-DNBOPM(x1@d2r?=aA!=!!DV$*gbu z%SiEGWINwMNk&-9`cqk7jz*Ac!)r2cL*EgMAcH;oCYorX|LE=HM1)BwsZGRa_L_t(YiS3iWYZO5g$3J^XHg}iY-mWHMA$k`608`j$6HqG= zi_=0ZQf$C<7Wo4NFRc(T!6rrEKoN^TY%PR4G=d`GL<>QbH4sQ~o~QT=8+~tQFFS#Y z#g%ejn3?aJna}Lo*|+dt>SV9rd~Nr&!;o*Jq3u-(88K9X%YU8O2im=c8Qh?f%*}vz zXd$>%+Ch3XGG#c@AFbM-mLIovVAFs?enr2cwtW@>X zxgZF_gnYX$w49xsclV2QPL-jES2Rt!;^CdNG5eXB4?;S!W51d7l`Kkf8fdo1_Zq(> zZOSH{qfbman|~(Rl9Hd`nhaZ4-9J@|D!bkH%pmCT|;(>ptV;+Ug4ruvdC-@+}uaD5|bN@%yOP# zU&)-DR%~0O6w4orbX%>p4RUfZt9Z{Dp3OX?8Lkq4PJiB%EOHOrkqeeoP9!sRQx&7Q z7ntw|YW#{1U5U~{PVHx5#e zRNkl4`a0PLnO2p-N)-FuR1D$Ab|kfz;4l%UzMJGZ`9VdfB6Dy_0_d+NRy^+EK>)00000NkvXX Hu0mjf4=fbE diff --git a/docs/html/img98.png b/docs/html/img98.png index 6e5db9a1d8556cfb44494c07e604f0f567ad0fd7..c3daca1c8672a3dc5e9aa1f66bc9deec2b86c4bf 100644 GIT binary patch delta 352 zcmV-m0iXWU1J?tPZ2?u0ZWn(jJ_d#$BvG~i28M?}1x1BWRlQ+gP(@P3sQ@JQGBAiJ zps4Zzy6^%*)l(qmI>o@iBETUiC@2P1#oP0Mfx!u(Dh?#Jfq{Y10og3wH$YVq5D}(E z1BMMCTPH9uupU5GC6r?Tu?7g%EBs>v8Nm&70e=A`BxFD)@HOxVFfe~Cff!Z9AjU?>EunhC`JkrXirFz|nX zsnSW{JI3_^5?s7U-ezD_VBiG0R~yW+Vc>6I0D6t<93oJ`-~~`2Q2afZBY7`^fngQ{ zo5(@I&;}$$OV2N06F*e=fZ;Bf=K+jThF%7o;2YSD$3P@MPDO0suwe#=BXXRg0&qHD yV1~G*3zsSs2{v4+5^%a@6pUzah6r$gI8^`)pm>Xue@~?V0000}6mO zQ9xG30#wC(0io(C5ObXZiVHLf3JNAbRoyrTB%Bbc;((ZU0|Nu21Cm)x_du#7AR>+rZ5=QSl}Kx7i3s delta 41 scmeyx_={1oGr-TCmrII^fq{Y7)59eQNGkxbJPR|BEUCD(f1=_|0J|ayLjV8( diff --git a/docs/html/index.html b/docs/html/index.html index 9c76bb1d..91af5ad7 100644 --- a/docs/html/index.html +++ b/docs/html/index.html @@ -1,4 +1,4 @@ - + - + next up previous - contents
- Next: Next: Abstract -   Contents
-
+
-MLD2P4 +MLD2P4

-User's and Reference Guide
-
A guide for the Multi-Level Domain Decomposition +User's and Reference Guide
+
A guide for the Multi-Level Domain Decomposition Parallel Preconditioners Package -based on PSBLAS
+based on PSBLAS


@@ -62,9 +61,9 @@ University of Rome ``Tor Vergata'', Italy


-Software version: 2.0 +Software version: 2.1
-Oct. 12, 2015 +Mar. 31, 2016
@@ -77,70 +76,70 @@ Oct. 12, 2015 -