diff --git a/mlprec/mld_caggrmap_bld.f90 b/mlprec/mld_caggrmap_bld.f90 index fdb01385..4acf8ac5 100644 --- a/mlprec/mld_caggrmap_bld.f90 +++ b/mlprec/mld_caggrmap_bld.f90 @@ -74,7 +74,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_inner_mod, mld_protect_name => mld_caggrmap_bld @@ -83,7 +83,8 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer, intent(in) :: aggr_type - type(psb_cspmat_type), intent(in), target :: a + real(psb_spk_), intent(in) :: theta + type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info @@ -91,9 +92,7 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Local variables integer, allocatable :: ils(:), neigh(:) integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m - type(psb_cspmat_type), target :: atmp, atrans - type(psb_cspmat_type), pointer :: apnt - + type(psb_cspmat_type) :: atmp, atrans logical :: recovery integer :: debug_level, debug_unit integer :: ictxt,np,me,err_act @@ -117,8 +116,93 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ncol = psb_cd_get_local_cols(desc_a) select case (aggr_type) - case (mld_dec_aggr_,mld_sym_dec_aggr_) - + case (mld_dec_aggr_) + call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + case (mld_sym_dec_aggr_) + nr = psb_sp_get_nrows(a) + call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_transp(atmp,atrans,fmt='COO') + if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_sp_free(atrans,info) + if (info == 0) call psb_ipcoo2csr(atmp,info) + if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == 0) call psb_sp_free(atmp,info) + + case default + + info = -1 + call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) + goto 9999 + end select + + if (info/=0) then + info=4001 + call psb_errpush(info,name,a_err='dec_map_bld') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +contains + + + + subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_inner_mod !, mld_protect_name => mld_daggrmap_bld + + implicit none + + ! Arguments + type(psb_cspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + + ! Local variables + integer, allocatable :: ils(:), neigh(:), irow(:), icol(:) + complex(psb_spk_), allocatable :: val(:), diag(:) + integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz + real(psb_spk_) :: cpling, tcl + logical :: recovery + integer :: debug_level, debug_unit + integer :: ictxt,np,me,err_act + integer :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=0 + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ! Note. At the time being we are ignoring aggr_type so + ! that we only have decoupled aggregation. This might + ! change in the future. + ! + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + nr = a%m allocate(ilaggr(nr),neigh(nr),stat=info) if(info /= 0) then @@ -128,36 +212,28 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if + allocate(diag(nr),stat=info) + if(info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + call psb_sp_getdiag(a,diag,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + do i=1, nr ilaggr(i) = -(nr+1) end do - if (aggr_type == mld_dec_aggr_) then - apnt => a - else - call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& - & rscale=.false.,cscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_transp(atmp,atrans,fmt='COO') - if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_sp_free(atrans,info) - if (info == 0) call psb_ipcoo2csr(atmp,info) - apnt => atmp - if (info/=0) then - info=4001 - call psb_errpush(info,name,a_err='init apnt') - goto 9999 - end if - - end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i - ! ! Phase one: group nodes together. ! Very simple minded strategy. @@ -172,26 +248,32 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! 1. Untouched nodes are marked >0 together ! with their neighbours ! - icnt = icnt + 1 - naggr = naggr + 1 + icnt = icnt + 1 + naggr = naggr + 1 ilaggr(i) = naggr - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do k=1, n_ne - j = neigh(k) - if ((1<=j).and.(j<=nr)) then - ilaggr(j) = naggr - endif + + do k=1, nz + j = icol(k) + if ((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 + end if enddo + ! ! 2. Untouched neighbours of these nodes are marked <0. ! - call psb_neigh(apnt,i,neigh,n_ne,info,lev=2) + call psb_neigh(a,i,neigh,n_ne,info,lev=2) if (info/=0) then info=4010 call psb_errpush(info,name,a_err='psb_neigh') @@ -218,7 +300,7 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! ! Phase two: sweep over leftovers. ! - allocate(ils(naggr+10),stat=info) + allocate(ils(nr),stat=info) if(info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& @@ -254,57 +336,55 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then ! ! Now some silly rule to break ties: - ! Group with smallest adjacent aggregate. + ! Group with adjacent aggregate. ! - isz = nr+1 - ia = -1 - - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + isz = nr+1 + ia = -1 + cpling = 0.0 + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do j=1, n_ne - k = neigh(j) - if ((1<=k).and.(k<=nr)) then - n = ilaggr(k) - if (n>0) then - if (n>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') - goto 9999 - end if - - if (ils(n) < isz) then - ia = n - isz = ils(n) + do j=1, nz + k = icol(j) + if ((1<=k).and.(k<=nr).and.(k/=i)) then + tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k))) + if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then +!!$ if (tcl > theta) then + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=4001 + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') + goto 9999 + end if + + if ((abs(val(j))>cpling) .or. & + & ((abs(val(j)) == cpling).and. (ils(n) < isz))) then +!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then + ia = n + isz = ils(n) + cpling = abs(val(j)) +!!$ cpling = tcl + endif endif endif - endif + end if enddo + if (ia == -1) then - if (ilaggr(i) > -(nr+1)) then - ilaggr(i) = abs(ilaggr(i)) - if (ilaggr(I)>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?') - goto 9999 - end if - ils(ilaggr(i)) = ils(ilaggr(i)) + 1 - ! - ! This might happen if the pattern is non symmetric. - ! Need a better handling. - ! - recovery = .true. - else - info=4001 - call psb_errpush(info,name,a_err='Unrecoverable error !!') - goto 9999 - endif + ! At this point, the easiest thing is to start a new aggregate + naggr = naggr + 1 + ilaggr(i) = naggr + ils(ilaggr(i)) = 1 + else + ilaggr(i) = ia + if (ia>naggr) then info=4001 call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ') @@ -312,8 +392,9 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) end if ils(ia) = ils(ia) + 1 endif + end if - enddo + end do if (debug_level >= psb_debug_outer_) then if (recovery) then write(debug_unit,*) me,' ',trim(name),& @@ -345,6 +426,12 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) call psb_errpush(info,name) goto 9999 end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=4001 + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if call psb_realloc(ncol,ilaggr,info) if (info/=0) then @@ -366,26 +453,17 @@ subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) nlaggr(me+1) = naggr call psb_sum(ictxt,nlaggr(1:np)) - if (aggr_type == mld_sym_dec_aggr_) then - call psb_sp_free(atmp,info) - end if - - case default - - info = -1 - call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return + call psb_erractionrestore(err_act) + return 9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if return - end if - return + + end subroutine mld_dec_map_bld end subroutine mld_caggrmap_bld diff --git a/mlprec/mld_cmlprec_bld.f90 b/mlprec/mld_cmlprec_bld.f90 index b164f37c..bf470c74 100644 --- a/mlprec/mld_cmlprec_bld.f90 +++ b/mlprec/mld_cmlprec_bld.f90 @@ -109,6 +109,7 @@ subroutine mld_cmlprec_bld(a,desc_a,p,info) call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',szero,is_legal_s_fact_thrs) end select call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',szero,is_legal_s_omega) + call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',szero,is_legal_s_aggr_thrs) call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',& & 1,is_legal_jac_sweeps) @@ -118,7 +119,8 @@ subroutine mld_cmlprec_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%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) + call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),& + & a,desc_a,p%nlaggr,p%mlia,info) if(info /= 0) then call psb_errpush(4010,name,a_err='mld_aggrmap_bld') goto 9999 diff --git a/mlprec/mld_cprecinit.f90 b/mlprec/mld_cprecinit.f90 index d8d0e3a4..78f4655b 100644 --- a/mlprec/mld_cprecinit.f90 +++ b/mlprec/mld_cprecinit.f90 @@ -189,7 +189,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_ p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_ @@ -204,7 +205,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ @@ -225,7 +227,8 @@ subroutine mld_cprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ diff --git a/mlprec/mld_cprecset.f90 b/mlprec/mld_cprecset.f90 index 0822bc02..875a9e09 100644 --- a/mlprec/mld_cprecset.f90 +++ b/mlprec/mld_cprecset.f90 @@ -595,7 +595,7 @@ subroutine mld_cprecsetr(p,what,val,info,ilev) else if (ilev_ > 1) then select case(what) - case(mld_aggr_damp_,mld_fact_thrs_) + case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_) p%baseprecv(ilev_)%rprcparm(what) = val case default write(0,*) name,': Error: invalid WHAT' @@ -610,7 +610,7 @@ subroutine mld_cprecsetr(p,what,val,info,ilev) select case(what) case(mld_fact_thrs_) - do ilev_=1,nlev_-1 + do ilev_=1,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 @@ -619,7 +619,16 @@ subroutine mld_cprecsetr(p,what,val,info,ilev) p%baseprecv(ilev_)%rprcparm(what) = val end do case(mld_aggr_damp_) - do ilev_=2,nlev_-1 + do ilev_=2,nlev_ + if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then + write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' + info = -1 + return + endif + p%baseprecv(ilev_)%rprcparm(what) = val + end do + case(mld_aggr_thresh_) + do ilev_=2,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 diff --git a/mlprec/mld_daggrmap_bld.f90 b/mlprec/mld_daggrmap_bld.f90 index b8ea8cdc..36aa1a3b 100644 --- a/mlprec/mld_daggrmap_bld.f90 +++ b/mlprec/mld_daggrmap_bld.f90 @@ -74,16 +74,17 @@ ! info - integer, output. ! Error code. ! -subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_inner_mod, mld_protect_name => mld_daggrmap_bld - + implicit none ! Arguments integer, intent(in) :: aggr_type - type(psb_dspmat_type), intent(in), target :: a + real(psb_dpk_), intent(in) :: theta + type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info @@ -91,8 +92,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Local variables integer, allocatable :: ils(:), neigh(:) integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m - type(psb_dspmat_type), target :: atmp, atrans - type(psb_dspmat_type), pointer :: apnt + type(psb_dspmat_type) :: atmp, atrans logical :: recovery integer :: debug_level, debug_unit integer :: ictxt,np,me,err_act @@ -116,8 +116,93 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ncol = psb_cd_get_local_cols(desc_a) select case (aggr_type) - case (mld_dec_aggr_,mld_sym_dec_aggr_) - + case (mld_dec_aggr_) + call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + case (mld_sym_dec_aggr_) + nr = psb_sp_get_nrows(a) + call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_transp(atmp,atrans,fmt='COO') + if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_sp_free(atrans,info) + if (info == 0) call psb_ipcoo2csr(atmp,info) + if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == 0) call psb_sp_free(atmp,info) + + case default + + info = -1 + call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) + goto 9999 + end select + + if (info/=0) then + info=4001 + call psb_errpush(info,name,a_err='dec_map_bld') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +contains + + + + subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_inner_mod !, mld_protect_name => mld_daggrmap_bld + + implicit none + + ! Arguments + type(psb_dspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + + ! Local variables + integer, allocatable :: ils(:), neigh(:), irow(:), icol(:) + real(psb_dpk_), allocatable :: val(:), diag(:) + integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz + real(psb_dpk_) :: cpling, tcl + logical :: recovery + integer :: debug_level, debug_unit + integer :: ictxt,np,me,err_act + integer :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=0 + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ! Note. At the time being we are ignoring aggr_type so + ! that we only have decoupled aggregation. This might + ! change in the future. + ! + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + nr = a%m allocate(ilaggr(nr),neigh(nr),stat=info) if(info /= 0) then @@ -127,36 +212,28 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if + allocate(diag(nr),stat=info) + if(info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + call psb_sp_getdiag(a,diag,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + do i=1, nr ilaggr(i) = -(nr+1) end do - if (aggr_type == mld_dec_aggr_) then - apnt => a - else - call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& - & rscale=.false.,cscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_transp(atmp,atrans,fmt='COO') - if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_sp_free(atrans,info) - if (info == 0) call psb_ipcoo2csr(atmp,info) - apnt => atmp - if (info/=0) then - info=4001 - call psb_errpush(info,name,a_err='init apnt') - goto 9999 - end if - - end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i - ! ! Phase one: group nodes together. ! Very simple minded strategy. @@ -171,26 +248,32 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! 1. Untouched nodes are marked >0 together ! with their neighbours ! - icnt = icnt + 1 - naggr = naggr + 1 + icnt = icnt + 1 + naggr = naggr + 1 ilaggr(i) = naggr - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do k=1, n_ne - j = neigh(k) - if ((1<=j).and.(j<=nr)) then - ilaggr(j) = naggr - endif + + do k=1, nz + j = icol(k) + if ((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 + end if enddo + ! ! 2. Untouched neighbours of these nodes are marked <0. ! - call psb_neigh(apnt,i,neigh,n_ne,info,lev=2) + call psb_neigh(a,i,neigh,n_ne,info,lev=2) if (info/=0) then info=4010 call psb_errpush(info,name,a_err='psb_neigh') @@ -217,7 +300,7 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! ! Phase two: sweep over leftovers. ! - allocate(ils(naggr+10),stat=info) + allocate(ils(nr),stat=info) if(info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& @@ -253,57 +336,55 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then ! ! Now some silly rule to break ties: - ! Group with smallest adjacent aggregate. + ! Group with adjacent aggregate. ! - isz = nr+1 - ia = -1 - - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + isz = nr+1 + ia = -1 + cpling = 0.0d0 + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do j=1, n_ne - k = neigh(j) - if ((1<=k).and.(k<=nr)) then - n = ilaggr(k) - if (n>0) then - if (n>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') - goto 9999 - end if - - if (ils(n) < isz) then - ia = n - isz = ils(n) + do j=1, nz + k = icol(j) + if ((1<=k).and.(k<=nr).and.(k/=i)) then + tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k))) + if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then +!!$ if (tcl > theta) then + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=4001 + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') + goto 9999 + end if + + if ((abs(val(j))>cpling) .or. & + & ((abs(val(j)) == cpling).and. (ils(n) < isz))) then +!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then + ia = n + isz = ils(n) + cpling = abs(val(j)) +!!$ cpling = tcl + endif endif endif - endif + end if enddo + if (ia == -1) then - if (ilaggr(i) > -(nr+1)) then - ilaggr(i) = abs(ilaggr(i)) - if (ilaggr(I)>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?') - goto 9999 - end if - ils(ilaggr(i)) = ils(ilaggr(i)) + 1 - ! - ! This might happen if the pattern is non symmetric. - ! Need a better handling. - ! - recovery = .true. - else - info=4001 - call psb_errpush(info,name,a_err='Unrecoverable error !!') - goto 9999 - endif + ! At this point, the easiest thing is to start a new aggregate + naggr = naggr + 1 + ilaggr(i) = naggr + ils(ilaggr(i)) = 1 + else + ilaggr(i) = ia + if (ia>naggr) then info=4001 call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ') @@ -311,8 +392,9 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) end if ils(ia) = ils(ia) + 1 endif + end if - enddo + end do if (debug_level >= psb_debug_outer_) then if (recovery) then write(debug_unit,*) me,' ',trim(name),& @@ -344,6 +426,12 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) call psb_errpush(info,name) goto 9999 end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=4001 + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if call psb_realloc(ncol,ilaggr,info) if (info/=0) then @@ -365,26 +453,17 @@ subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) nlaggr(me+1) = naggr call psb_sum(ictxt,nlaggr(1:np)) - if (aggr_type == mld_sym_dec_aggr_) then - call psb_sp_free(atmp,info) - end if - - case default - - info = -1 - call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return + call psb_erractionrestore(err_act) + return 9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if return - end if - return + + end subroutine mld_dec_map_bld end subroutine mld_daggrmap_bld diff --git a/mlprec/mld_dmlprec_bld.f90 b/mlprec/mld_dmlprec_bld.f90 index 356a8bcc..89a33ef4 100644 --- a/mlprec/mld_dmlprec_bld.f90 +++ b/mlprec/mld_dmlprec_bld.f90 @@ -109,6 +109,7 @@ subroutine mld_dmlprec_bld(a,desc_a,p,info) call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',dzero,is_legal_fact_thrs) end select call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',dzero,is_legal_omega) + call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',dzero,is_legal_aggr_thrs) call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',& & 1,is_legal_jac_sweeps) @@ -118,7 +119,8 @@ subroutine mld_dmlprec_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%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) + call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),& + & a,desc_a,p%nlaggr,p%mlia,info) if(info /= 0) then call psb_errpush(4010,name,a_err='mld_aggrmap_bld') goto 9999 diff --git a/mlprec/mld_dprecinit.f90 b/mlprec/mld_dprecinit.f90 index e0d97dad..7192a526 100644 --- a/mlprec/mld_dprecinit.f90 +++ b/mlprec/mld_dprecinit.f90 @@ -189,7 +189,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_ p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_ @@ -204,7 +205,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ @@ -225,7 +227,8 @@ subroutine mld_dprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ diff --git a/mlprec/mld_dprecset.f90 b/mlprec/mld_dprecset.f90 index 05b68476..60906a53 100644 --- a/mlprec/mld_dprecset.f90 +++ b/mlprec/mld_dprecset.f90 @@ -595,7 +595,7 @@ subroutine mld_dprecsetr(p,what,val,info,ilev) else if (ilev_ > 1) then select case(what) - case(mld_aggr_damp_,mld_fact_thrs_) + case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_) p%baseprecv(ilev_)%rprcparm(what) = val case default write(0,*) name,': Error: invalid WHAT' @@ -610,7 +610,7 @@ subroutine mld_dprecsetr(p,what,val,info,ilev) select case(what) case(mld_fact_thrs_) - do ilev_=1,nlev_-1 + do ilev_=1,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 @@ -619,7 +619,16 @@ subroutine mld_dprecsetr(p,what,val,info,ilev) p%baseprecv(ilev_)%rprcparm(what) = val end do case(mld_aggr_damp_) - do ilev_=2,nlev_-1 + do ilev_=2,nlev_ + if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then + write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' + info = -1 + return + endif + p%baseprecv(ilev_)%rprcparm(what) = val + end do + case(mld_aggr_thresh_) + do ilev_=2,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 diff --git a/mlprec/mld_inner_mod.f90 b/mlprec/mld_inner_mod.f90 index 83cf1e4d..b9728f91 100644 --- a/mlprec/mld_inner_mod.f90 +++ b/mlprec/mld_inner_mod.f90 @@ -384,38 +384,42 @@ module mld_inner_mod end interface interface mld_aggrmap_bld - subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ use mld_prec_type, only : mld_sbaseprc_type integer, intent(in) :: aggr_type - type(psb_sspmat_type), intent(in), target :: a + real(psb_spk_), intent(in) :: theta + type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info end subroutine mld_saggrmap_bld - subroutine mld_daggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ use mld_prec_type, only : mld_dbaseprc_type integer, intent(in) :: aggr_type - type(psb_dspmat_type), intent(in), target :: a + real(psb_dpk_), intent(in) :: theta + type(psb_dspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info end subroutine mld_daggrmap_bld - subroutine mld_caggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ use mld_prec_type, only : mld_cbaseprc_type integer, intent(in) :: aggr_type - type(psb_cspmat_type), intent(in), target :: a + real(psb_spk_), intent(in) :: theta + type(psb_cspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info end subroutine mld_caggrmap_bld - subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) + subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ use mld_prec_type, only : mld_zbaseprc_type integer, intent(in) :: aggr_type - type(psb_zspmat_type), intent(in), target :: a + real(psb_dpk_), intent(in) :: theta + type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info diff --git a/mlprec/mld_prec_type.f90 b/mlprec/mld_prec_type.f90 index 367cbdc4..4ee2a513 100644 --- a/mlprec/mld_prec_type.f90 +++ b/mlprec/mld_prec_type.f90 @@ -66,7 +66,9 @@ module mld_prec_type & psb_dspmat_type, psb_zspmat_type,& & psb_sspmat_type, psb_cspmat_type,& & psb_desc_type, psb_inter_desc_type, psb_sizeof, psb_dpk_, psb_spk_,& - & psb_sp_free, psb_cdfree + & psb_sp_free, psb_cdfree, psb_halo_, psb_none_, & + & psb_nohalo_, psb_square_root_, & + & psb_sizeof_int, psb_sizeof_sp, psb_sizeof_dp, psb_sizeof ! ! Type: mld_dprec_type, mld_zprec_type @@ -718,6 +720,8 @@ contains & aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_)) write(iout,*) 'Aggregation smoothing: ', & & aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_)) + write(iout,*) 'Aggregation threshold: ', & + & p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_) if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then write(iout,*) 'Damping omega: ', & & p%baseprecv(ilev)%rprcparm(mld_aggr_damp_) @@ -822,6 +826,8 @@ contains & aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_)) write(iout,*) 'Aggregation smoothing: ', & & aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_)) + write(iout,*) 'Aggregation threshold: ', & + & p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_) if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then write(iout,*) 'Damping omega: ', & & p%baseprecv(ilev)%rprcparm(mld_aggr_damp_) @@ -944,8 +950,10 @@ contains if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then write(iout,*) 'Multilevel aggregation: ', & & aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_)) - write(iout,*) 'Smoother: ', & + write(iout,*) 'Aggregation smoothing: ', & & aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_)) + write(iout,*) 'Aggregation threshold: ', & + & p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_) if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then write(iout,*) 'Smoothing omega: ', & & p%baseprecv(ilev)%rprcparm(mld_aggr_damp_) @@ -1047,8 +1055,10 @@ contains if (p%baseprecv(ilev)%iprcparm(mld_ml_type_)>mld_no_ml_) then write(iout,*) 'Multilevel aggregation: ', & & aggr_names(p%baseprecv(ilev)%iprcparm(mld_aggr_alg_)) - write(iout,*) 'Smoother: ', & + write(iout,*) 'Aggregation smoothing: ', & & aggr_kinds(p%baseprecv(ilev)%iprcparm(mld_aggr_kind_)) + write(iout,*) 'Aggregation threshold: ', & + & p%baseprecv(ilev)%rprcparm(mld_aggr_thresh_) if (p%baseprecv(ilev)%iprcparm(mld_aggr_kind_) /= mld_no_smooth_) then write(iout,*) 'Smoothing omega: ', & & p%baseprecv(ilev)%rprcparm(mld_aggr_damp_) @@ -1205,6 +1215,13 @@ contains is_legal_fact_thrs = (ip>=0.0d0) return end function is_legal_fact_thrs + function is_legal_aggr_thrs(ip) + real(psb_dpk_), intent(in) :: ip + logical :: is_legal_aggr_thrs + + is_legal_aggr_thrs = (ip>=0.0d0) + return + end function is_legal_aggr_thrs function is_legal_s_omega(ip) real(psb_spk_), intent(in) :: ip @@ -1219,6 +1236,13 @@ contains is_legal_s_fact_thrs = (ip>=0.0) return end function is_legal_s_fact_thrs + function is_legal_s_aggr_thrs(ip) + real(psb_spk_), intent(in) :: ip + logical :: is_legal_s_aggr_thrs + + is_legal_s_aggr_thrs = (ip>=0.0) + return + end function is_legal_s_aggr_thrs subroutine mld_icheck_def(ip,name,id,is_legal) diff --git a/mlprec/mld_saggrmap_bld.f90 b/mlprec/mld_saggrmap_bld.f90 index f534b541..1ea6ec50 100644 --- a/mlprec/mld_saggrmap_bld.f90 +++ b/mlprec/mld_saggrmap_bld.f90 @@ -74,7 +74,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_inner_mod, mld_protect_name => mld_saggrmap_bld @@ -83,7 +83,8 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer, intent(in) :: aggr_type - type(psb_sspmat_type), intent(in), target :: a + real(psb_spk_), intent(in) :: theta + type(psb_sspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info @@ -91,8 +92,7 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Local variables integer, allocatable :: ils(:), neigh(:) integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m - type(psb_sspmat_type), target :: atmp, atrans - type(psb_sspmat_type), pointer :: apnt + type(psb_sspmat_type) :: atmp, atrans logical :: recovery integer :: debug_level, debug_unit integer :: ictxt,np,me,err_act @@ -116,8 +116,93 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ncol = psb_cd_get_local_cols(desc_a) select case (aggr_type) - case (mld_dec_aggr_,mld_sym_dec_aggr_) - + case (mld_dec_aggr_) + call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + case (mld_sym_dec_aggr_) + nr = psb_sp_get_nrows(a) + call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_transp(atmp,atrans,fmt='COO') + if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_sp_free(atrans,info) + if (info == 0) call psb_ipcoo2csr(atmp,info) + if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == 0) call psb_sp_free(atmp,info) + + case default + + info = -1 + call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) + goto 9999 + end select + + if (info/=0) then + info=4001 + call psb_errpush(info,name,a_err='dec_map_bld') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +contains + + + + subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_inner_mod !, mld_protect_name => mld_daggrmap_bld + + implicit none + + ! Arguments + type(psb_sspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_spk_), intent(in) :: theta + integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + + ! Local variables + integer, allocatable :: ils(:), neigh(:), irow(:), icol(:) + real(psb_spk_), allocatable :: val(:), diag(:) + integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz + real(psb_spk_) :: cpling, tcl + logical :: recovery + integer :: debug_level, debug_unit + integer :: ictxt,np,me,err_act + integer :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=0 + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ! Note. At the time being we are ignoring aggr_type so + ! that we only have decoupled aggregation. This might + ! change in the future. + ! + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + nr = a%m allocate(ilaggr(nr),neigh(nr),stat=info) if(info /= 0) then @@ -127,36 +212,28 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if + allocate(diag(nr),stat=info) + if(info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + call psb_sp_getdiag(a,diag,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + do i=1, nr ilaggr(i) = -(nr+1) end do - if (aggr_type == mld_dec_aggr_) then - apnt => a - else - call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& - & rscale=.false.,cscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_transp(atmp,atrans,fmt='COO') - if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_sp_free(atrans,info) - if (info == 0) call psb_ipcoo2csr(atmp,info) - apnt => atmp - if (info/=0) then - info=4001 - call psb_errpush(info,name,a_err='init apnt') - goto 9999 - end if - - end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i - ! ! Phase one: group nodes together. ! Very simple minded strategy. @@ -171,26 +248,32 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! 1. Untouched nodes are marked >0 together ! with their neighbours ! - icnt = icnt + 1 - naggr = naggr + 1 + icnt = icnt + 1 + naggr = naggr + 1 ilaggr(i) = naggr - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do k=1, n_ne - j = neigh(k) - if ((1<=j).and.(j<=nr)) then - ilaggr(j) = naggr - endif + + do k=1, nz + j = icol(k) + if ((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 + end if enddo + ! ! 2. Untouched neighbours of these nodes are marked <0. ! - call psb_neigh(apnt,i,neigh,n_ne,info,lev=2) + call psb_neigh(a,i,neigh,n_ne,info,lev=2) if (info/=0) then info=4010 call psb_errpush(info,name,a_err='psb_neigh') @@ -217,7 +300,7 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! ! Phase two: sweep over leftovers. ! - allocate(ils(naggr+10),stat=info) + allocate(ils(nr),stat=info) if(info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& @@ -253,57 +336,55 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then ! ! Now some silly rule to break ties: - ! Group with smallest adjacent aggregate. + ! Group with adjacent aggregate. ! - isz = nr+1 - ia = -1 - - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + isz = nr+1 + ia = -1 + cpling = 0.0 + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do j=1, n_ne - k = neigh(j) - if ((1<=k).and.(k<=nr)) then - n = ilaggr(k) - if (n>0) then - if (n>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') - goto 9999 - end if - - if (ils(n) < isz) then - ia = n - isz = ils(n) + do j=1, nz + k = icol(j) + if ((1<=k).and.(k<=nr).and.(k/=i)) then + tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k))) + if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then +!!$ if (tcl > theta) then + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=4001 + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') + goto 9999 + end if + + if ((abs(val(j))>cpling) .or. & + & ((abs(val(j)) == cpling).and. (ils(n) < isz))) then +!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then + ia = n + isz = ils(n) + cpling = abs(val(j)) +!!$ cpling = tcl + endif endif endif - endif + end if enddo + if (ia == -1) then - if (ilaggr(i) > -(nr+1)) then - ilaggr(i) = abs(ilaggr(i)) - if (ilaggr(I)>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?') - goto 9999 - end if - ils(ilaggr(i)) = ils(ilaggr(i)) + 1 - ! - ! This might happen if the pattern is non symmetric. - ! Need a better handling. - ! - recovery = .true. - else - info=4001 - call psb_errpush(info,name,a_err='Unrecoverable error !!') - goto 9999 - endif + ! At this point, the easiest thing is to start a new aggregate + naggr = naggr + 1 + ilaggr(i) = naggr + ils(ilaggr(i)) = 1 + else + ilaggr(i) = ia + if (ia>naggr) then info=4001 call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ') @@ -311,8 +392,9 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) end if ils(ia) = ils(ia) + 1 endif + end if - enddo + end do if (debug_level >= psb_debug_outer_) then if (recovery) then write(debug_unit,*) me,' ',trim(name),& @@ -344,6 +426,12 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) call psb_errpush(info,name) goto 9999 end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=4001 + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if call psb_realloc(ncol,ilaggr,info) if (info/=0) then @@ -365,26 +453,17 @@ subroutine mld_saggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) nlaggr(me+1) = naggr call psb_sum(ictxt,nlaggr(1:np)) - if (aggr_type == mld_sym_dec_aggr_) then - call psb_sp_free(atmp,info) - end if - - case default - - info = -1 - call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return + call psb_erractionrestore(err_act) + return 9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if return - end if - return + + end subroutine mld_dec_map_bld end subroutine mld_saggrmap_bld diff --git a/mlprec/mld_smlprec_bld.f90 b/mlprec/mld_smlprec_bld.f90 index 63185b66..dfebec10 100644 --- a/mlprec/mld_smlprec_bld.f90 +++ b/mlprec/mld_smlprec_bld.f90 @@ -109,6 +109,7 @@ subroutine mld_smlprec_bld(a,desc_a,p,info) call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',szero,is_legal_s_fact_thrs) end select call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',szero,is_legal_s_omega) + call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',szero,is_legal_s_aggr_thrs) call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',& & 1,is_legal_jac_sweeps) @@ -118,7 +119,8 @@ subroutine mld_smlprec_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%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) + call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),& + & a,desc_a,p%nlaggr,p%mlia,info) if(info /= 0) then call psb_errpush(4010,name,a_err='mld_aggrmap_bld') goto 9999 diff --git a/mlprec/mld_sprecinit.f90 b/mlprec/mld_sprecinit.f90 index bb3a41fe..d92248f1 100644 --- a/mlprec/mld_sprecinit.f90 +++ b/mlprec/mld_sprecinit.f90 @@ -189,7 +189,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_ p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_ @@ -204,7 +205,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ @@ -225,7 +227,8 @@ subroutine mld_sprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = szero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ diff --git a/mlprec/mld_sprecset.f90 b/mlprec/mld_sprecset.f90 index 613661bf..0a9a2dcb 100644 --- a/mlprec/mld_sprecset.f90 +++ b/mlprec/mld_sprecset.f90 @@ -595,7 +595,7 @@ subroutine mld_sprecsetr(p,what,val,info,ilev) else if (ilev_ > 1) then select case(what) - case(mld_aggr_damp_,mld_fact_thrs_) + case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_) p%baseprecv(ilev_)%rprcparm(what) = val case default write(0,*) name,': Error: invalid WHAT' @@ -610,7 +610,7 @@ subroutine mld_sprecsetr(p,what,val,info,ilev) select case(what) case(mld_fact_thrs_) - do ilev_=1,nlev_-1 + do ilev_=1,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 @@ -619,7 +619,16 @@ subroutine mld_sprecsetr(p,what,val,info,ilev) p%baseprecv(ilev_)%rprcparm(what) = val end do case(mld_aggr_damp_) - do ilev_=2,nlev_-1 + do ilev_=2,nlev_ + if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then + write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' + info = -1 + return + endif + p%baseprecv(ilev_)%rprcparm(what) = val + end do + case(mld_aggr_thresh_) + do ilev_=2,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 diff --git a/mlprec/mld_zaggrmap_bld.f90 b/mlprec/mld_zaggrmap_bld.f90 index a289767d..4349fa57 100644 --- a/mlprec/mld_zaggrmap_bld.f90 +++ b/mlprec/mld_zaggrmap_bld.f90 @@ -74,7 +74,7 @@ ! info - integer, output. ! Error code. ! -subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) +subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,nlaggr,ilaggr,info) use psb_base_mod use mld_inner_mod, mld_protect_name => mld_zaggrmap_bld @@ -83,7 +83,8 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Arguments integer, intent(in) :: aggr_type - type(psb_zspmat_type), intent(in), target :: a + real(psb_dpk_), intent(in) :: theta + type(psb_zspmat_type), intent(in) :: a type(psb_desc_type), intent(in) :: desc_a integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) integer, intent(out) :: info @@ -91,9 +92,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! Local variables integer, allocatable :: ils(:), neigh(:) integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m - type(psb_zspmat_type), target :: atmp, atrans - type(psb_zspmat_type), pointer :: apnt - + type(psb_zspmat_type) :: atmp, atrans logical :: recovery integer :: debug_level, debug_unit integer :: ictxt,np,me,err_act @@ -117,8 +116,93 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ncol = psb_cd_get_local_cols(desc_a) select case (aggr_type) - case (mld_dec_aggr_,mld_sym_dec_aggr_) - + case (mld_dec_aggr_) + call mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + case (mld_sym_dec_aggr_) + nr = psb_sp_get_nrows(a) + call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& + & rscale=.false.,cscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_transp(atmp,atrans,fmt='COO') + if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) + atmp%m=nr + atmp%k=nr + if (info == 0) call psb_sp_free(atrans,info) + if (info == 0) call psb_ipcoo2csr(atmp,info) + if (info == 0) call mld_dec_map_bld(theta,atmp,desc_a,nlaggr,ilaggr,info) + if (info == 0) call psb_sp_free(atmp,info) + + case default + + info = -1 + call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) + goto 9999 + end select + + if (info/=0) then + info=4001 + call psb_errpush(info,name,a_err='dec_map_bld') + goto 9999 + end if + + call psb_erractionrestore(err_act) + return + +9999 continue + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if + return + +contains + + + + subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_inner_mod !, mld_protect_name => mld_daggrmap_bld + + implicit none + + ! Arguments + type(psb_zspmat_type), intent(in) :: a + type(psb_desc_type), intent(in) :: desc_a + real(psb_dpk_), intent(in) :: theta + integer, allocatable, intent(out) :: ilaggr(:),nlaggr(:) + integer, intent(out) :: info + + ! Local variables + integer, allocatable :: ils(:), neigh(:), irow(:), icol(:) + complex(psb_dpk_), allocatable :: val(:), diag(:) + integer :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz + real(psb_dpk_) :: cpling, tcl + logical :: recovery + integer :: debug_level, debug_unit + integer :: ictxt,np,me,err_act + integer :: nrow, ncol, n_ne + character(len=20) :: name, ch_err + + if(psb_get_errstatus() /= 0) return + info=0 + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + ! Note. At the time being we are ignoring aggr_type so + ! that we only have decoupled aggregation. This might + ! change in the future. + ! + ictxt=psb_cd_get_context(desc_a) + call psb_info(ictxt,me,np) + nrow = psb_cd_get_local_rows(desc_a) + ncol = psb_cd_get_local_cols(desc_a) + nr = a%m allocate(ilaggr(nr),neigh(nr),stat=info) if(info /= 0) then @@ -128,36 +212,28 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) goto 9999 end if + allocate(diag(nr),stat=info) + if(info /= 0) then + info=4025 + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + call psb_sp_getdiag(a,diag,info) + if(info /= 0) then + info=4010 + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + do i=1, nr ilaggr(i) = -(nr+1) end do - if (aggr_type == mld_dec_aggr_) then - apnt => a - else - call psb_sp_clip(a,atmp,info,imax=nr,jmax=nr,& - & rscale=.false.,cscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_transp(atmp,atrans,fmt='COO') - if (info == 0) call psb_rwextd(nr,atmp,info,b=atrans,rowscale=.false.) - atmp%m=nr - atmp%k=nr - if (info == 0) call psb_sp_free(atrans,info) - if (info == 0) call psb_ipcoo2csr(atmp,info) - apnt => atmp - if (info/=0) then - info=4001 - call psb_errpush(info,name,a_err='init apnt') - goto 9999 - end if - - end if ! Note: -(nr+1) Untouched as yet ! -i 1<=i<=nr Adjacent to aggregate i ! i 1<=i<=nr Belonging to aggregate i - ! ! Phase one: group nodes together. ! Very simple minded strategy. @@ -172,26 +248,32 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! 1. Untouched nodes are marked >0 together ! with their neighbours ! - icnt = icnt + 1 - naggr = naggr + 1 + icnt = icnt + 1 + naggr = naggr + 1 ilaggr(i) = naggr - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do k=1, n_ne - j = neigh(k) - if ((1<=j).and.(j<=nr)) then - ilaggr(j) = naggr - endif + + do k=1, nz + j = icol(k) + if ((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 + end if enddo + ! ! 2. Untouched neighbours of these nodes are marked <0. ! - call psb_neigh(apnt,i,neigh,n_ne,info,lev=2) + call psb_neigh(a,i,neigh,n_ne,info,lev=2) if (info/=0) then info=4010 call psb_errpush(info,name,a_err='psb_neigh') @@ -218,7 +300,7 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) ! ! Phase two: sweep over leftovers. ! - allocate(ils(naggr+10),stat=info) + allocate(ils(nr),stat=info) if(info /= 0) then info=4025 call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& @@ -254,57 +336,55 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) if (ilaggr(i) < 0) then ! ! Now some silly rule to break ties: - ! Group with smallest adjacent aggregate. + ! Group with adjacent aggregate. ! - isz = nr+1 - ia = -1 - - call psb_neigh(apnt,i,neigh,n_ne,info,lev=1) + isz = nr+1 + ia = -1 + cpling = 0.0d0 + call psb_sp_getrow(i,a,nz,irow,icol,val,info) if (info/=0) then info=4010 - call psb_errpush(info,name,a_err='psb_neigh') + call psb_errpush(info,name,a_err='psb_sp_getrow') goto 9999 end if - do j=1, n_ne - k = neigh(j) - if ((1<=k).and.(k<=nr)) then - n = ilaggr(k) - if (n>0) then - if (n>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') - goto 9999 - end if - - if (ils(n) < isz) then - ia = n - isz = ils(n) + do j=1, nz + k = icol(j) + if ((1<=k).and.(k<=nr).and.(k/=i)) then + tcl = abs(val(j)) / sqrt(abs(diag(i)*diag(k))) + if (abs(val(j)) > theta*sqrt(abs(diag(i)*diag(k)))) then +!!$ if (tcl > theta) then + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=4001 + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 2 ?') + goto 9999 + end if + + if ((abs(val(j))>cpling) .or. & + & ((abs(val(j)) == cpling).and. (ils(n) < isz))) then +!!$ if ((tcl > cpling) .or. ((tcl == cpling).and. (ils(n) < isz))) then + ia = n + isz = ils(n) + cpling = abs(val(j)) +!!$ cpling = tcl + endif endif endif - endif + end if enddo + if (ia == -1) then - if (ilaggr(i) > -(nr+1)) then - ilaggr(i) = abs(ilaggr(i)) - if (ilaggr(I)>naggr) then - info=4001 - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 3 ?') - goto 9999 - end if - ils(ilaggr(i)) = ils(ilaggr(i)) + 1 - ! - ! This might happen if the pattern is non symmetric. - ! Need a better handling. - ! - recovery = .true. - else - info=4001 - call psb_errpush(info,name,a_err='Unrecoverable error !!') - goto 9999 - endif + ! At this point, the easiest thing is to start a new aggregate + naggr = naggr + 1 + ilaggr(i) = naggr + ils(ilaggr(i)) = 1 + else + ilaggr(i) = ia + if (ia>naggr) then info=4001 call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 4? ') @@ -312,8 +392,9 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) end if ils(ia) = ils(ia) + 1 endif + end if - enddo + end do if (debug_level >= psb_debug_outer_) then if (recovery) then write(debug_unit,*) me,' ',trim(name),& @@ -345,6 +426,12 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) call psb_errpush(info,name) goto 9999 end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=4001 + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if call psb_realloc(ncol,ilaggr,info) if (info/=0) then @@ -366,26 +453,17 @@ subroutine mld_zaggrmap_bld(aggr_type,a,desc_a,nlaggr,ilaggr,info) nlaggr(me+1) = naggr call psb_sum(ictxt,nlaggr(1:np)) - if (aggr_type == mld_sym_dec_aggr_) then - call psb_sp_free(atmp,info) - end if - - case default - - info = -1 - call psb_errpush(30,name,i_err=(/1,aggr_type,0,0,0/)) - goto 9999 - end select - - call psb_erractionrestore(err_act) - return + call psb_erractionrestore(err_act) + return 9999 continue - call psb_erractionrestore(err_act) - if (err_act.eq.psb_act_abort_) then - call psb_error() + call psb_erractionrestore(err_act) + if (err_act.eq.psb_act_abort_) then + call psb_error() + return + end if return - end if - return + + end subroutine mld_dec_map_bld end subroutine mld_zaggrmap_bld diff --git a/mlprec/mld_zmlprec_bld.f90 b/mlprec/mld_zmlprec_bld.f90 index 57af6de3..12b13728 100644 --- a/mlprec/mld_zmlprec_bld.f90 +++ b/mlprec/mld_zmlprec_bld.f90 @@ -109,6 +109,7 @@ subroutine mld_zmlprec_bld(a,desc_a,p,info) call mld_check_def(p%rprcparm(mld_fact_thrs_),'Eps',dzero,is_legal_fact_thrs) end select call mld_check_def(p%rprcparm(mld_aggr_damp_),'Omega',dzero,is_legal_omega) + call mld_check_def(p%rprcparm(mld_aggr_thresh_),'Aggr_Thresh',dzero,is_legal_aggr_thrs) call mld_check_def(p%iprcparm(mld_smooth_sweeps_),'Jacobi sweeps',& & 1,is_legal_jac_sweeps) @@ -118,7 +119,8 @@ subroutine mld_zmlprec_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%iprcparm(mld_aggr_alg_),a,desc_a,p%nlaggr,p%mlia,info) + call mld_aggrmap_bld(p%iprcparm(mld_aggr_alg_),p%rprcparm(mld_aggr_thresh_),& + & a,desc_a,p%nlaggr,p%mlia,info) if(info /= 0) then call psb_errpush(4010,name,a_err='mld_aggrmap_bld') goto 9999 diff --git a/mlprec/mld_zprecinit.f90 b/mlprec/mld_zprecinit.f90 index 5b2370eb..0f6b0d2c 100644 --- a/mlprec/mld_zprecinit.f90 +++ b/mlprec/mld_zprecinit.f90 @@ -189,7 +189,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_as_ p%baseprecv(ilev_)%iprcparm(mld_sub_solve_) = mld_ilu_n_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_halo_ @@ -204,7 +205,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ @@ -225,7 +227,8 @@ subroutine mld_zprecinit(p,ptype,info,nlev) if (info == 0) call psb_realloc(mld_ifpsz_,p%baseprecv(ilev_)%iprcparm,info) if (info == 0) call psb_realloc(mld_rfpsz_,p%baseprecv(ilev_)%rprcparm,info) if (info /= 0) return - p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%iprcparm(:) = 0 + p%baseprecv(ilev_)%rprcparm(:) = dzero p%baseprecv(ilev_)%iprcparm(mld_prec_type_) = mld_bjac_ p%baseprecv(ilev_)%iprcparm(mld_sub_restr_) = psb_none_ p%baseprecv(ilev_)%iprcparm(mld_sub_prol_) = psb_none_ diff --git a/mlprec/mld_zprecset.f90 b/mlprec/mld_zprecset.f90 index e999b2a6..963d4827 100644 --- a/mlprec/mld_zprecset.f90 +++ b/mlprec/mld_zprecset.f90 @@ -595,7 +595,7 @@ subroutine mld_zprecsetr(p,what,val,info,ilev) else if (ilev_ > 1) then select case(what) - case(mld_aggr_damp_,mld_fact_thrs_) + case(mld_aggr_damp_,mld_aggr_thresh_,mld_fact_thrs_) p%baseprecv(ilev_)%rprcparm(what) = val case default write(0,*) name,': Error: invalid WHAT' @@ -610,7 +610,7 @@ subroutine mld_zprecsetr(p,what,val,info,ilev) select case(what) case(mld_fact_thrs_) - do ilev_=1,nlev_-1 + do ilev_=1,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 @@ -619,7 +619,16 @@ subroutine mld_zprecsetr(p,what,val,info,ilev) p%baseprecv(ilev_)%rprcparm(what) = val end do case(mld_aggr_damp_) - do ilev_=2,nlev_-1 + do ilev_=2,nlev_ + if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then + write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' + info = -1 + return + endif + p%baseprecv(ilev_)%rprcparm(what) = val + end do + case(mld_aggr_thresh_) + do ilev_=2,nlev_ if (.not.allocated(p%baseprecv(ilev_)%rprcparm)) then write(0,*) name,': Error: Uninitialized preconditioner component, should call MLD_PRECINIT' info = -1 diff --git a/test/fileread/cf_sample.f90 b/test/fileread/cf_sample.f90 index 162b1b55..bcdbab83 100644 --- a/test/fileread/cf_sample.f90 +++ b/test/fileread/cf_sample.f90 @@ -67,6 +67,7 @@ program cf_sample real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_spk_) :: omega ! smoother omega + real(psb_spk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prec_choice @@ -242,7 +243,7 @@ program cf_sample call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info) call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv) call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info) call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) if (prec_choice%omega>=0.0) then call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv) @@ -392,6 +393,7 @@ contains call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prec%cjswp,5) ! Jacobi sweeps call read_data(prec%omega,5) ! smoother omega + call read_data(prec%athres,5) ! smoother aggr thresh end if end if @@ -426,6 +428,7 @@ contains call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps call psb_bcast(icontxt,prec%omega) ! smoother omega + call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh end if end subroutine get_parms diff --git a/test/fileread/df_bench.f90 b/test/fileread/df_bench.f90 index b22d4e17..4a8999e9 100644 --- a/test/fileread/df_bench.f90 +++ b/test/fileread/df_bench.f90 @@ -29,6 +29,7 @@ program df_bench real(psb_dpk_) :: thr2 ! Threshold for fact. 1 ILU(T) integer :: jswp ! Jacobi sweeps real(psb_dpk_) :: omega ! smoother omega + real(psb_dpk_) :: athres ! smoother aggregation character(len=40) :: descr ! verbose description of the prec end type precdata type(precdata), allocatable :: precs(:) diff --git a/test/fileread/df_sample.f90 b/test/fileread/df_sample.f90 index a1d350a7..19a197b2 100644 --- a/test/fileread/df_sample.f90 +++ b/test/fileread/df_sample.f90 @@ -67,6 +67,7 @@ program df_sample real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_dpk_) :: omega ! smoother omega + real(psb_dpk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prec_choice @@ -242,7 +243,7 @@ program df_sample call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info) call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv) call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info) call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) if (prec_choice%omega>=0.0) then call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv) @@ -392,6 +393,7 @@ contains call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prec%cjswp,5) ! Jacobi sweeps call read_data(prec%omega,5) ! smoother omega + call read_data(prec%athres,5) ! smoother aggr thresh end if end if @@ -426,6 +428,7 @@ contains call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps call psb_bcast(icontxt,prec%omega) ! smoother omega + call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh end if end subroutine get_parms diff --git a/test/fileread/runs/cfs.inp b/test/fileread/runs/cfs.inp index b8df9d7c..6d764c6e 100644 --- a/test/fileread/runs/cfs.inp +++ b/test/fileread/runs/cfs.inp @@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST 1.d-4 ! Coarse level: Threshold T for ILU(T,P) 4 ! Coarse level: Number of Jacobi sweeps -1.0d0 ! Smoother Omega: if < 0 means library choice. +0.01d0 ! Smoother Aggregation Threshold: >= 0.0 \ No newline at end of file diff --git a/test/fileread/runs/dfs.inp b/test/fileread/runs/dfs.inp index 2a9f9528..34b16a76 100644 --- a/test/fileread/runs/dfs.inp +++ b/test/fileread/runs/dfs.inp @@ -1,5 +1,5 @@ -thm1000x600.mtx ! les_t4.mtx ! young1r.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or -NONE !les_t4.rhs ! rhs.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html +matphi_140x33x45.mtx !A_1M_gps.mtx !thm1000x600.mtx ! les_t4.mtx ! young1r.mtx This (and others) from: http://math.nist.gov/MatrixMarket/ or +tnoto_phi.mtx !NONE !les_t4.rhs ! rhs.mtx http://www.cise.ufl.edu/research/sparse/matrices/index.html BICGSTAB ! Iterative method: BiCGSTAB BiCG CGS RGMRES BiCGSTABL CG CSR ! Storage format CSR COO JAD 0 ! IPART: Partition method 0: BLK 2: graph (with Metis) @@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST 1.d-4 ! Coarse level: Threshold T for ILU(T,P) 4 ! Coarse level: Number of Jacobi sweeps -1.0d0 ! Smoother Omega: if < 0 means library choice. +0.01d0 ! Smoother Aggregation Threshold: >= 0.0 diff --git a/test/fileread/runs/drt.sh b/test/fileread/runs/drt.sh index e88d3a5c..15c642a0 100755 --- a/test/fileread/runs/drt.sh +++ b/test/fileread/runs/drt.sh @@ -16,7 +16,7 @@ do # 3rd batch: 14bis,64bis 4 sweeps echo "mpirun -np $np -machinefile locm df_bench >>log.part$part.ren$renum.${np}p" #/usr/local/mpich-gcc42/bin/mpirun -np $np -machinefile locm df_bench <>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <>run.gfc.kiva.log.${np}p.$date 2>err.gfc.kiva.log.${np}p.$date <= 0.0 \ No newline at end of file diff --git a/test/fileread/runs/zfs.inp b/test/fileread/runs/zfs.inp index b8df9d7c..6d764c6e 100644 --- a/test/fileread/runs/zfs.inp +++ b/test/fileread/runs/zfs.inp @@ -27,3 +27,4 @@ ILU ! Coarse level: solver ILU ILUT UMF SLU SLUDIST 1.d-4 ! Coarse level: Threshold T for ILU(T,P) 4 ! Coarse level: Number of Jacobi sweeps -1.0d0 ! Smoother Omega: if < 0 means library choice. +0.01d0 ! Smoother Aggregation Threshold: >= 0.0 \ No newline at end of file diff --git a/test/fileread/sf_sample.f90 b/test/fileread/sf_sample.f90 index 88ab6a0d..6145ec5e 100644 --- a/test/fileread/sf_sample.f90 +++ b/test/fileread/sf_sample.f90 @@ -67,6 +67,7 @@ program sf_sample real(psb_spk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_spk_) :: omega ! smoother omega + real(psb_spk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prec_choice @@ -242,7 +243,7 @@ program sf_sample call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info) call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv) call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info) call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) if (prec_choice%omega>=0.0) then call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv) @@ -392,6 +393,7 @@ contains call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prec%cjswp,5) ! Jacobi sweeps call read_data(prec%omega,5) ! smoother omega + call read_data(prec%athres,5) ! smoother aggr thresh end if end if @@ -426,6 +428,7 @@ contains call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps call psb_bcast(icontxt,prec%omega) ! smoother omega + call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh end if end subroutine get_parms diff --git a/test/fileread/zf_sample.f90 b/test/fileread/zf_sample.f90 index 7e12667d..a4838947 100644 --- a/test/fileread/zf_sample.f90 +++ b/test/fileread/zf_sample.f90 @@ -67,6 +67,7 @@ program zf_sample real(psb_dpk_) :: cthres ! Threshold for fact. 1 ILU(T) integer :: cjswp ! Jacobi sweeps real(psb_dpk_) :: omega ! smoother omega + real(psb_dpk_) :: athres ! smoother aggregation threshold end type precdata type(precdata) :: prec_choice @@ -242,7 +243,7 @@ program zf_sample call mld_precset(prec,mld_coarse_solve_,prec_choice%csolve,info) call mld_precset(prec,mld_sub_fill_in_,prec_choice%cfill,info,ilev=nlv) call mld_precset(prec,mld_fact_thrs_,prec_choice%cthres,info,ilev=nlv) - call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) + call mld_precset(prec,mld_aggr_thresh_,prec_choice%athres,info) call mld_precset(prec,mld_smooth_sweeps_,prec_choice%cjswp,info,ilev=nlv) if (prec_choice%omega>=0.0) then call mld_precset(prec,mld_aggr_damp_,prec_choice%omega,info,ilev=nlv) @@ -392,6 +393,7 @@ contains call read_data(prec%cthres,5) ! Threshold for fact. 1 ILU(T) call read_data(prec%cjswp,5) ! Jacobi sweeps call read_data(prec%omega,5) ! smoother omega + call read_data(prec%athres,5) ! smoother aggr thresh end if end if @@ -426,6 +428,7 @@ contains call psb_bcast(icontxt,prec%cthres) ! Threshold for fact. 1 ILU(T) call psb_bcast(icontxt,prec%cjswp) ! Jacobi sweeps call psb_bcast(icontxt,prec%omega) ! smoother omega + call psb_bcast(icontxt,prec%athres) ! smoother aggr thresh end if end subroutine get_parms