diff --git a/mlprec/Makefile b/mlprec/Makefile index 1bb42258..ca14987b 100644 --- a/mlprec/Makefile +++ b/mlprec/Makefile @@ -36,22 +36,22 @@ MPCOBJS=mld_sslud_interface.o mld_dslud_interface.o mld_cslud_interface.o mld_zs SINNEROBJS= mld_scoarse_bld.o mld_smlprec_bld.o \ mld_silu0_fact.o mld_siluk_fact.o mld_silut_fact.o mld_saggrmap_bld.o \ - mld_smlprec_aply.o mld_saggrmat_asb.o \ + mld_s_dec_map_bld.o mld_smlprec_aply.o mld_saggrmat_asb.o \ $(SMPFOBJS) DINNEROBJS= mld_dcoarse_bld.o mld_dmlprec_bld.o \ mld_dilu0_fact.o mld_diluk_fact.o mld_dilut_fact.o mld_daggrmap_bld.o \ - mld_dmlprec_aply.o mld_daggrmat_asb.o \ + mld_d_dec_map_bld.o mld_dmlprec_aply.o mld_daggrmat_asb.o \ $(DMPFOBJS) CINNEROBJS= mld_ccoarse_bld.o mld_cmlprec_bld.o \ mld_cilu0_fact.o mld_ciluk_fact.o mld_cilut_fact.o mld_caggrmap_bld.o \ - mld_cmlprec_aply.o mld_caggrmat_asb.o \ + mld_c_dec_map_bld.o mld_cmlprec_aply.o mld_caggrmat_asb.o \ $(CMPFOBJS) ZINNEROBJS= mld_zcoarse_bld.o mld_zmlprec_bld.o \ mld_zilu0_fact.o mld_ziluk_fact.o mld_zilut_fact.o mld_zaggrmap_bld.o \ - mld_zmlprec_aply.o mld_zaggrmat_asb.o \ + mld_z_dec_map_bld.o mld_zmlprec_aply.o mld_zaggrmat_asb.o \ $(ZMPFOBJS) INNEROBJS= $(SINNEROBJS) $(DINNEROBJS) $(CINNEROBJS) $(ZINNEROBJS) diff --git a/mlprec/mld_c_dec_map_bld.F90 b/mlprec/mld_c_dec_map_bld.F90 new file mode 100644 index 00000000..f50aebaf --- /dev/null +++ b/mlprec/mld_c_dec_map_bld.F90 @@ -0,0 +1,300 @@ + +subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_c_inner_mod, mld_protect_name => mld_c_dec_map_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=psb_success_ + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + 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%get_nrows() + allocate(ilaggr(nr),neigh(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + allocate(diag(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='complex(psb_spk_)') + goto 9999 + end if + call a%get_diag(diag,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + do i=1, nr + ilaggr(i) = -(nr+1) + end do + + + ! 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. + ! + naggr = 0 + nlp = 0 + do + icnt = 0 + do i=1, nr + if (ilaggr(i) == -(nr+1)) then + ! + ! 1. Untouched nodes are marked >0 together + ! with their neighbours + ! + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='csget') + goto 9999 + end if + + 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 a%get_neigh(i,neigh,n_ne,info,lev=2) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_neigh') + goto 9999 + end if + + do n = 1, n_ne + m = neigh(n) + if ((1<=m).and.(m<=nr)) then + if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr + endif + enddo + endif + enddo + nlp = nlp + 1 + if (icnt == 0) exit + enddo + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: sweep over leftovers. + ! + allocate(ils(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + do i=1, size(ils) + ils(i) = 0 + end do + do i=1, nr + n = ilaggr(i) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') + goto 9999 + else + ils(n) = ils(n) + 1 + end if + + end if + end do + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: number of aggregates ',naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: nodes aggregated ',sum(ils) + end if + + recovery=.false. + do i=1, nr + if (ilaggr(i) < 0) then + ! + ! Now some silly rule to break ties: + ! Group with adjacent aggregate. + ! + isz = nr+1 + ia = -1 + cpling = szero + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getrow') + goto 9999 + end if + + 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=psb_err_internal_error_ + 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 + end if + enddo + + if (ia == -1) then + ! 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=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') + goto 9999 + end if + ils(ia) = ils(ia) + 1 + endif + + end if + end do + if (debug_level >= psb_debug_outer_) then + if (recovery) then + write(debug_unit,*) me,' ',trim(name),& + & 'Had to recover from strange situation in loc_aggregate.' + write(debug_unit,*) me,' ',trim(name),& + & 'Perhaps an unsymmetric pattern?' + endif + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 2: number of aggregates ',naggr,sum(ils) + do i=1, naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) + enddo + write(debug_unit,*) me,' ',trim(name),& + & maxval(ils(1:naggr)) + write(debug_unit,*) me,' ',trim(name),& + & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' + end if + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + deallocate(ils,neigh,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + 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 + +end subroutine mld_c_dec_map_bld + diff --git a/mlprec/mld_c_inner_mod.f90 b/mlprec/mld_c_inner_mod.f90 index 937c8891..5377ae67 100644 --- a/mlprec/mld_c_inner_mod.f90 +++ b/mlprec/mld_c_inner_mod.f90 @@ -102,6 +102,18 @@ module mld_c_inner_mod end subroutine mld_caggrmap_bld end interface mld_aggrmap_bld + + interface mld_dec_map_bld + subroutine mld_c_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ + 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 + end subroutine mld_c_dec_map_bld + end interface mld_dec_map_bld + interface mld_aggrmat_asb subroutine mld_caggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_cspmat_type, psb_desc_type, psb_spk_ diff --git a/mlprec/mld_caggrmap_bld.f90 b/mlprec/mld_caggrmap_bld.f90 index 21078230..6ed756cd 100644 --- a/mlprec/mld_caggrmap_bld.f90 +++ b/mlprec/mld_caggrmap_bld.f90 @@ -160,305 +160,4 @@ subroutine mld_caggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) end if return -contains - - subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) - - use psb_base_mod - use mld_c_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=psb_success_ - name = 'mld_dec_map_bld' - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ! - 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%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - allocate(diag(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& - & a_err='complex(psb_spk_)') - goto 9999 - end if - call a%get_diag(diag,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getdiag') - goto 9999 - end if - - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - - ! 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. - ! - naggr = 0 - nlp = 0 - do - icnt = 0 - do i=1, nr - if (ilaggr(i) == -(nr+1)) then - ! - ! 1. Untouched nodes are marked >0 together - ! with their neighbours - ! - icnt = icnt + 1 - naggr = naggr + 1 - ilaggr(i) = naggr - - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='csget') - goto 9999 - end if - - 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 a%get_neigh(i,neigh,n_ne,info,lev=2) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_neigh') - goto 9999 - end if - - do n = 1, n_ne - m = neigh(n) - if ((1<=m).and.(m<=nr)) then - if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr - endif - enddo - endif - enddo - nlp = nlp + 1 - if (icnt == 0) exit - enddo - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & ' Check 1:',count(ilaggr == -(nr+1)) - end if - - ! - ! Phase two: sweep over leftovers. - ! - allocate(ils(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - do i=1, size(ils) - ils(i) = 0 - end do - do i=1, nr - n = ilaggr(i) - if (n>0) then - if (n>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') - goto 9999 - else - ils(n) = ils(n) + 1 - end if - - end if - end do - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: number of aggregates ',naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: nodes aggregated ',sum(ils) - end if - - recovery=.false. - do i=1, nr - if (ilaggr(i) < 0) then - ! - ! Now some silly rule to break ties: - ! Group with adjacent aggregate. - ! - isz = nr+1 - ia = -1 - cpling = szero - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getrow') - goto 9999 - end if - - 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=psb_err_internal_error_ - 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 - end if - enddo - - if (ia == -1) then - ! 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=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 - end if - ils(ia) = ils(ia) + 1 - endif - - end if - end do - if (debug_level >= psb_debug_outer_) then - if (recovery) then - write(debug_unit,*) me,' ',trim(name),& - & 'Had to recover from strange situation in loc_aggregate.' - write(debug_unit,*) me,' ',trim(name),& - & 'Perhaps an unsymmetric pattern?' - endif - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 2: number of aggregates ',naggr,sum(ils) - do i=1, naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) - enddo - write(debug_unit,*) me,' ',trim(name),& - & maxval(ils(1:naggr)) - write(debug_unit,*) me,' ',trim(name),& - & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' - end if - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif - - deallocate(ils,neigh,stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - end if - if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') - goto 9999 - end if - - call psb_realloc(ncol,ilaggr,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(nlaggr(np),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - nlaggr(:) = 0 - nlaggr(me+1) = naggr - call psb_sum(ictxt,nlaggr(1:np)) - - 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 - - end subroutine mld_dec_map_bld - end subroutine mld_caggrmap_bld diff --git a/mlprec/mld_d_dec_map_bld.F90 b/mlprec/mld_d_dec_map_bld.F90 new file mode 100644 index 00000000..e82eb2a4 --- /dev/null +++ b/mlprec/mld_d_dec_map_bld.F90 @@ -0,0 +1,300 @@ + +subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_d_inner_mod, mld_protect_name => mld_d_dec_map_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=psb_success_ + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + 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%get_nrows() + allocate(ilaggr(nr),neigh(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + allocate(diag(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='real(psb_dpk_)') + goto 9999 + end if + call a%get_diag(diag,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + do i=1, nr + ilaggr(i) = -(nr+1) + end do + + + ! 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. + ! + naggr = 0 + nlp = 0 + do + icnt = 0 + do i=1, nr + if (ilaggr(i) == -(nr+1)) then + ! + ! 1. Untouched nodes are marked >0 together + ! with their neighbours + ! + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='csget') + goto 9999 + end if + + 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 a%get_neigh(i,neigh,n_ne,info,lev=2) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_neigh') + goto 9999 + end if + + do n = 1, n_ne + m = neigh(n) + if ((1<=m).and.(m<=nr)) then + if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr + endif + enddo + endif + enddo + nlp = nlp + 1 + if (icnt == 0) exit + enddo + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: sweep over leftovers. + ! + allocate(ils(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + do i=1, size(ils) + ils(i) = 0 + end do + do i=1, nr + n = ilaggr(i) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') + goto 9999 + else + ils(n) = ils(n) + 1 + end if + + end if + end do + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: number of aggregates ',naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: nodes aggregated ',sum(ils) + end if + + recovery=.false. + do i=1, nr + if (ilaggr(i) < 0) then + ! + ! Now some silly rule to break ties: + ! Group with adjacent aggregate. + ! + isz = nr+1 + ia = -1 + cpling = dzero + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getrow') + goto 9999 + end if + + 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=psb_err_internal_error_ + 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 + end if + enddo + + if (ia == -1) then + ! 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=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') + goto 9999 + end if + ils(ia) = ils(ia) + 1 + endif + + end if + end do + if (debug_level >= psb_debug_outer_) then + if (recovery) then + write(debug_unit,*) me,' ',trim(name),& + & 'Had to recover from strange situation in loc_aggregate.' + write(debug_unit,*) me,' ',trim(name),& + & 'Perhaps an unsymmetric pattern?' + endif + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 2: number of aggregates ',naggr,sum(ils) + do i=1, naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) + enddo + write(debug_unit,*) me,' ',trim(name),& + & maxval(ils(1:naggr)) + write(debug_unit,*) me,' ',trim(name),& + & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' + end if + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + deallocate(ils,neigh,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + 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 + +end subroutine mld_d_dec_map_bld + diff --git a/mlprec/mld_d_ilu_solver.f90 b/mlprec/mld_d_ilu_solver.f90 index 8e772a63..1ce37117 100644 --- a/mlprec/mld_d_ilu_solver.f90 +++ b/mlprec/mld_d_ilu_solver.f90 @@ -407,8 +407,6 @@ contains if (present(mold)) then call sv%l%cscnv(info,mold=mold) call sv%u%cscnv(info,mold=mold) - Write(0,*) 'Converted L into ',sv%l%get_fmt(),& - &' and U into ',sv%u%get_fmt() end if if (debug_level >= psb_debug_outer_) & diff --git a/mlprec/mld_d_inner_mod.f90 b/mlprec/mld_d_inner_mod.f90 index 4ca75e03..24f5b855 100644 --- a/mlprec/mld_d_inner_mod.f90 +++ b/mlprec/mld_d_inner_mod.f90 @@ -104,6 +104,18 @@ module mld_d_inner_mod end subroutine mld_daggrmap_bld end interface mld_aggrmap_bld + + interface mld_dec_map_bld + subroutine mld_d_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ + 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 + end subroutine mld_d_dec_map_bld + end interface mld_dec_map_bld + interface mld_aggrmat_asb subroutine mld_daggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_dspmat_type, psb_desc_type, psb_dpk_ diff --git a/mlprec/mld_daggrmap_bld.f90 b/mlprec/mld_daggrmap_bld.f90 index 9387d202..6297ae80 100644 --- a/mlprec/mld_daggrmap_bld.f90 +++ b/mlprec/mld_daggrmap_bld.f90 @@ -160,305 +160,4 @@ subroutine mld_daggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) end if return -contains - - subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) - - use psb_base_mod - use mld_d_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=psb_success_ - name = 'mld_dec_map_bld' - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ! - 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%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - allocate(diag(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& - & a_err='real(psb_dpk_)') - goto 9999 - end if - call a%get_diag(diag,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getdiag') - goto 9999 - end if - - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - - ! 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. - ! - naggr = 0 - nlp = 0 - do - icnt = 0 - do i=1, nr - if (ilaggr(i) == -(nr+1)) then - ! - ! 1. Untouched nodes are marked >0 together - ! with their neighbours - ! - icnt = icnt + 1 - naggr = naggr + 1 - ilaggr(i) = naggr - - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='csget') - goto 9999 - end if - - 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 a%get_neigh(i,neigh,n_ne,info,lev=2) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_neigh') - goto 9999 - end if - - do n = 1, n_ne - m = neigh(n) - if ((1<=m).and.(m<=nr)) then - if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr - endif - enddo - endif - enddo - nlp = nlp + 1 - if (icnt == 0) exit - enddo - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & ' Check 1:',count(ilaggr == -(nr+1)) - end if - - ! - ! Phase two: sweep over leftovers. - ! - allocate(ils(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - do i=1, size(ils) - ils(i) = 0 - end do - do i=1, nr - n = ilaggr(i) - if (n>0) then - if (n>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') - goto 9999 - else - ils(n) = ils(n) + 1 - end if - - end if - end do - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: number of aggregates ',naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: nodes aggregated ',sum(ils) - end if - - recovery=.false. - do i=1, nr - if (ilaggr(i) < 0) then - ! - ! Now some silly rule to break ties: - ! Group with adjacent aggregate. - ! - isz = nr+1 - ia = -1 - cpling = dzero - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getrow') - goto 9999 - end if - - 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=psb_err_internal_error_ - 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 - end if - enddo - - if (ia == -1) then - ! 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=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 - end if - ils(ia) = ils(ia) + 1 - endif - - end if - end do - if (debug_level >= psb_debug_outer_) then - if (recovery) then - write(debug_unit,*) me,' ',trim(name),& - & 'Had to recover from strange situation in loc_aggregate.' - write(debug_unit,*) me,' ',trim(name),& - & 'Perhaps an unsymmetric pattern?' - endif - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 2: number of aggregates ',naggr,sum(ils) - do i=1, naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) - enddo - write(debug_unit,*) me,' ',trim(name),& - & maxval(ils(1:naggr)) - write(debug_unit,*) me,' ',trim(name),& - & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' - end if - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif - - deallocate(ils,neigh,stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - end if - if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') - goto 9999 - end if - - call psb_realloc(ncol,ilaggr,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(nlaggr(np),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - nlaggr(:) = 0 - nlaggr(me+1) = naggr - call psb_sum(ictxt,nlaggr(1:np)) - - 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 - - end subroutine mld_dec_map_bld - end subroutine mld_daggrmap_bld diff --git a/mlprec/mld_s_dec_map_bld.F90 b/mlprec/mld_s_dec_map_bld.F90 new file mode 100644 index 00000000..2cbfefeb --- /dev/null +++ b/mlprec/mld_s_dec_map_bld.F90 @@ -0,0 +1,296 @@ + +subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_s_inner_mod, mld_protect_name => mld_s_dec_map_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=psb_success_ + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + 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%get_nrows() + allocate(ilaggr(nr),neigh(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + allocate(diag(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='real(psb_spk_)') + goto 9999 + end if + call a%get_diag(diag,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + do i=1, nr + ilaggr(i) = -(nr+1) + end do + + + ! 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. + ! + naggr = 0 + nlp = 0 + do + icnt = 0 + do i=1, nr + if (ilaggr(i) == -(nr+1)) then + ! + ! 1. Untouched nodes are marked >0 together + ! with their neighbours + ! + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='csget') + goto 9999 + end if + + 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 a%get_neigh(i,neigh,n_ne,info,lev=2) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_neigh') + goto 9999 + end if + + do n = 1, n_ne + m = neigh(n) + if ((1<=m).and.(m<=nr)) then + if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr + endif + enddo + endif + enddo + nlp = nlp + 1 + if (icnt == 0) exit + enddo + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: sweep over leftovers. + ! + allocate(ils(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + do i=1, size(ils) + ils(i) = 0 + end do + do i=1, nr + n = ilaggr(i) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') + goto 9999 + else + ils(n) = ils(n) + 1 + end if + + end if + end do + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: number of aggregates ',naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: nodes aggregated ',sum(ils) + end if + + recovery=.false. + do i=1, nr + if (ilaggr(i) < 0) then + ! + ! Now some silly rule to break ties: + ! Group with adjacent aggregate. + ! + isz = nr+1 + ia = -1 + cpling = szero + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getrow') + goto 9999 + end if + + 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 + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + 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 + ia = n + isz = ils(n) + cpling = abs(val(j)) + endif + endif + endif + end if + enddo + + if (ia == -1) then + ! 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=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') + goto 9999 + end if + ils(ia) = ils(ia) + 1 + endif + + end if + end do + if (debug_level >= psb_debug_outer_) then + if (recovery) then + write(debug_unit,*) me,' ',trim(name),& + & 'Had to recover from strange situation in loc_aggregate.' + write(debug_unit,*) me,' ',trim(name),& + & 'Perhaps an unsymmetric pattern?' + endif + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 2: number of aggregates ',naggr,sum(ils) + do i=1, naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) + enddo + write(debug_unit,*) me,' ',trim(name),& + & maxval(ils(1:naggr)) + write(debug_unit,*) me,' ',trim(name),& + & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' + end if + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + deallocate(ils,neigh,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + 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 + +end subroutine mld_s_dec_map_bld diff --git a/mlprec/mld_s_inner_mod.f90 b/mlprec/mld_s_inner_mod.f90 index d176aa72..2a049ec8 100644 --- a/mlprec/mld_s_inner_mod.f90 +++ b/mlprec/mld_s_inner_mod.f90 @@ -102,6 +102,17 @@ module mld_s_inner_mod end subroutine mld_saggrmap_bld end interface mld_aggrmap_bld + interface mld_dec_map_bld + subroutine mld_s_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ + 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 + end subroutine mld_s_dec_map_bld + end interface mld_dec_map_bld + interface mld_aggrmat_asb subroutine mld_saggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_sspmat_type, psb_desc_type, psb_spk_ diff --git a/mlprec/mld_saggrmap_bld.f90 b/mlprec/mld_saggrmap_bld.f90 index ac8e12eb..4ad92ed7 100644 --- a/mlprec/mld_saggrmap_bld.f90 +++ b/mlprec/mld_saggrmap_bld.f90 @@ -160,305 +160,4 @@ subroutine mld_saggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) end if return -contains - - subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) - - use psb_base_mod - use mld_s_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=psb_success_ - name = 'mld_dec_map_bld' - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ! - 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%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - allocate(diag(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& - & a_err='real(psb_spk_)') - goto 9999 - end if - call a%get_diag(diag,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getdiag') - goto 9999 - end if - - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - - ! 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. - ! - naggr = 0 - nlp = 0 - do - icnt = 0 - do i=1, nr - if (ilaggr(i) == -(nr+1)) then - ! - ! 1. Untouched nodes are marked >0 together - ! with their neighbours - ! - icnt = icnt + 1 - naggr = naggr + 1 - ilaggr(i) = naggr - - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='csget') - goto 9999 - end if - - 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 a%get_neigh(i,neigh,n_ne,info,lev=2) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_neigh') - goto 9999 - end if - - do n = 1, n_ne - m = neigh(n) - if ((1<=m).and.(m<=nr)) then - if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr - endif - enddo - endif - enddo - nlp = nlp + 1 - if (icnt == 0) exit - enddo - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & ' Check 1:',count(ilaggr == -(nr+1)) - end if - - ! - ! Phase two: sweep over leftovers. - ! - allocate(ils(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - do i=1, size(ils) - ils(i) = 0 - end do - do i=1, nr - n = ilaggr(i) - if (n>0) then - if (n>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') - goto 9999 - else - ils(n) = ils(n) + 1 - end if - - end if - end do - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: number of aggregates ',naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: nodes aggregated ',sum(ils) - end if - - recovery=.false. - do i=1, nr - if (ilaggr(i) < 0) then - ! - ! Now some silly rule to break ties: - ! Group with adjacent aggregate. - ! - isz = nr+1 - ia = -1 - cpling = szero - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getrow') - goto 9999 - end if - - 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=psb_err_internal_error_ - 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 - end if - enddo - - if (ia == -1) then - ! 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=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 - end if - ils(ia) = ils(ia) + 1 - endif - - end if - end do - if (debug_level >= psb_debug_outer_) then - if (recovery) then - write(debug_unit,*) me,' ',trim(name),& - & 'Had to recover from strange situation in loc_aggregate.' - write(debug_unit,*) me,' ',trim(name),& - & 'Perhaps an unsymmetric pattern?' - endif - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 2: number of aggregates ',naggr,sum(ils) - do i=1, naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) - enddo - write(debug_unit,*) me,' ',trim(name),& - & maxval(ils(1:naggr)) - write(debug_unit,*) me,' ',trim(name),& - & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' - end if - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif - - deallocate(ils,neigh,stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - end if - if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') - goto 9999 - end if - - call psb_realloc(ncol,ilaggr,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(nlaggr(np),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - nlaggr(:) = 0 - nlaggr(me+1) = naggr - call psb_sum(ictxt,nlaggr(1:np)) - - 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 - - end subroutine mld_dec_map_bld - end subroutine mld_saggrmap_bld diff --git a/mlprec/mld_z_dec_map_bld.F90 b/mlprec/mld_z_dec_map_bld.F90 new file mode 100644 index 00000000..51e29e0e --- /dev/null +++ b/mlprec/mld_z_dec_map_bld.F90 @@ -0,0 +1,296 @@ + +subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + + use psb_base_mod + use mld_z_inner_mod, mld_protect_name => mld_z_dec_map_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=psb_success_ + name = 'mld_dec_map_bld' + call psb_erractionsave(err_act) + debug_unit = psb_get_debug_unit() + debug_level = psb_get_debug_level() + ! + 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%get_nrows() + allocate(ilaggr(nr),neigh(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + allocate(diag(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& + & a_err='complex(psb_dpk_)') + goto 9999 + end if + call a%get_diag(diag,info) + if(info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getdiag') + goto 9999 + end if + + do i=1, nr + ilaggr(i) = -(nr+1) + end do + + + ! 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. + ! + naggr = 0 + nlp = 0 + do + icnt = 0 + do i=1, nr + if (ilaggr(i) == -(nr+1)) then + ! + ! 1. Untouched nodes are marked >0 together + ! with their neighbours + ! + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='csget') + goto 9999 + end if + + 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 a%get_neigh(i,neigh,n_ne,info,lev=2) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_neigh') + goto 9999 + end if + + do n = 1, n_ne + m = neigh(n) + if ((1<=m).and.(m<=nr)) then + if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr + endif + enddo + endif + enddo + nlp = nlp + 1 + if (icnt == 0) exit + enddo + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & ' Check 1:',count(ilaggr == -(nr+1)) + end if + + ! + ! Phase two: sweep over leftovers. + ! + allocate(ils(nr),stat=info) + if(info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + do i=1, size(ils) + ils(i) = 0 + end do + do i=1, nr + n = ilaggr(i) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') + goto 9999 + else + ils(n) = ils(n) + 1 + end if + + end if + end do + if (debug_level >= psb_debug_outer_) then + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: number of aggregates ',naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 1: nodes aggregated ',sum(ils) + end if + + recovery=.false. + do i=1, nr + if (ilaggr(i) < 0) then + ! + ! Now some silly rule to break ties: + ! Group with adjacent aggregate. + ! + isz = nr+1 + ia = -1 + cpling = dzero + call a%csget(i,i,nz,irow,icol,val,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + call psb_errpush(info,name,a_err='psb_sp_getrow') + goto 9999 + end if + + 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 + n = ilaggr(k) + if (n>0) then + if (n>naggr) then + info=psb_err_internal_error_ + 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 + ia = n + isz = ils(n) + cpling = abs(val(j)) + endif + endif + endif + end if + enddo + + if (ia == -1) then + ! 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=psb_err_internal_error_ + call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') + goto 9999 + end if + ils(ia) = ils(ia) + 1 + endif + + end if + end do + if (debug_level >= psb_debug_outer_) then + if (recovery) then + write(debug_unit,*) me,' ',trim(name),& + & 'Had to recover from strange situation in loc_aggregate.' + write(debug_unit,*) me,' ',trim(name),& + & 'Perhaps an unsymmetric pattern?' + endif + write(debug_unit,*) me,' ',trim(name),& + & 'Phase 2: number of aggregates ',naggr,sum(ils) + do i=1, naggr + write(debug_unit,*) me,' ',trim(name),& + & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) + enddo + write(debug_unit,*) me,' ',trim(name),& + & maxval(ils(1:naggr)) + write(debug_unit,*) me,' ',trim(name),& + & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' + end if + + if (count(ilaggr<0) >0) then + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: some leftovers') + goto 9999 + endif + + deallocate(ils,neigh,stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_dealloc_ + call psb_errpush(info,name) + goto 9999 + end if + if (naggr > ncol) then + write(0,*) name,'Error : naggr > ncol',naggr,ncol + info=psb_err_internal_error_ + call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') + goto 9999 + end if + + call psb_realloc(ncol,ilaggr,info) + if (info /= psb_success_) then + info=psb_err_from_subroutine_ + ch_err='psb_realloc' + call psb_errpush(info,name,a_err=ch_err) + goto 9999 + end if + + allocate(nlaggr(np),stat=info) + if (info /= psb_success_) then + info=psb_err_alloc_request_ + call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& + & a_err='integer') + goto 9999 + end if + + nlaggr(:) = 0 + nlaggr(me+1) = naggr + call psb_sum(ictxt,nlaggr(1:np)) + + 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 + +end subroutine mld_z_dec_map_bld diff --git a/mlprec/mld_z_inner_mod.f90 b/mlprec/mld_z_inner_mod.f90 index 4c4efa94..b1662896 100644 --- a/mlprec/mld_z_inner_mod.f90 +++ b/mlprec/mld_z_inner_mod.f90 @@ -102,6 +102,18 @@ module mld_z_inner_mod end subroutine mld_zaggrmap_bld end interface mld_aggrmap_bld + + interface mld_dec_map_bld + subroutine mld_z_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) + use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ + 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 + end subroutine mld_z_dec_map_bld + end interface mld_dec_map_bld + interface mld_aggrmat_asb subroutine mld_zaggrmat_asb(a,desc_a,ilaggr,nlaggr,p,info) use psb_base_mod, only : psb_zspmat_type, psb_desc_type, psb_dpk_ diff --git a/mlprec/mld_zaggrmap_bld.f90 b/mlprec/mld_zaggrmap_bld.f90 index 8dbf1c50..3ca1374d 100644 --- a/mlprec/mld_zaggrmap_bld.f90 +++ b/mlprec/mld_zaggrmap_bld.f90 @@ -160,305 +160,4 @@ subroutine mld_zaggrmap_bld(aggr_type,theta,a,desc_a,ilaggr,nlaggr,info) end if return -contains - - subroutine mld_dec_map_bld(theta,a,desc_a,nlaggr,ilaggr,info) - - use psb_base_mod - use mld_z_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=psb_success_ - name = 'mld_dec_map_bld' - call psb_erractionsave(err_act) - debug_unit = psb_get_debug_unit() - debug_level = psb_get_debug_level() - ! - 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%get_nrows() - allocate(ilaggr(nr),neigh(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/2*nr,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - allocate(diag(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/nr,0,0,0,0/),& - & a_err='complex(psb_dpk_)') - goto 9999 - end if - call a%get_diag(diag,info) - if(info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getdiag') - goto 9999 - end if - - do i=1, nr - ilaggr(i) = -(nr+1) - end do - - - ! 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. - ! - naggr = 0 - nlp = 0 - do - icnt = 0 - do i=1, nr - if (ilaggr(i) == -(nr+1)) then - ! - ! 1. Untouched nodes are marked >0 together - ! with their neighbours - ! - icnt = icnt + 1 - naggr = naggr + 1 - ilaggr(i) = naggr - - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='csget') - goto 9999 - end if - - 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 a%get_neigh(i,neigh,n_ne,info,lev=2) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_neigh') - goto 9999 - end if - - do n = 1, n_ne - m = neigh(n) - if ((1<=m).and.(m<=nr)) then - if (ilaggr(m) == -(nr+1)) ilaggr(m) = -naggr - endif - enddo - endif - enddo - nlp = nlp + 1 - if (icnt == 0) exit - enddo - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & ' Check 1:',count(ilaggr == -(nr+1)) - end if - - ! - ! Phase two: sweep over leftovers. - ! - allocate(ils(nr),stat=info) - if(info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/naggr+10,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - do i=1, size(ils) - ils(i) = 0 - end do - do i=1, nr - n = ilaggr(i) - if (n>0) then - if (n>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr 1 ?') - goto 9999 - else - ils(n) = ils(n) + 1 - end if - - end if - end do - if (debug_level >= psb_debug_outer_) then - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: number of aggregates ',naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 1: nodes aggregated ',sum(ils) - end if - - recovery=.false. - do i=1, nr - if (ilaggr(i) < 0) then - ! - ! Now some silly rule to break ties: - ! Group with adjacent aggregate. - ! - isz = nr+1 - ia = -1 - cpling = dzero - call a%csget(i,i,nz,irow,icol,val,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - call psb_errpush(info,name,a_err='psb_sp_getrow') - goto 9999 - end if - - 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=psb_err_internal_error_ - 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 - end if - enddo - - if (ia == -1) then - ! 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=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 - end if - ils(ia) = ils(ia) + 1 - endif - - end if - end do - if (debug_level >= psb_debug_outer_) then - if (recovery) then - write(debug_unit,*) me,' ',trim(name),& - & 'Had to recover from strange situation in loc_aggregate.' - write(debug_unit,*) me,' ',trim(name),& - & 'Perhaps an unsymmetric pattern?' - endif - write(debug_unit,*) me,' ',trim(name),& - & 'Phase 2: number of aggregates ',naggr,sum(ils) - do i=1, naggr - write(debug_unit,*) me,' ',trim(name),& - & 'Size of aggregate ',i,' :',count(ilaggr == i), ils(i) - enddo - write(debug_unit,*) me,' ',trim(name),& - & maxval(ils(1:naggr)) - write(debug_unit,*) me,' ',trim(name),& - & 'Leftovers ',count(ilaggr<0), ' in ',nlp,' loops' - end if - - if (count(ilaggr<0) >0) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: some leftovers') - goto 9999 - endif - - deallocate(ils,neigh,stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_dealloc_ - call psb_errpush(info,name) - goto 9999 - end if - if (naggr > ncol) then - write(0,*) name,'Error : naggr > ncol',naggr,ncol - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='Fatal error: naggr>ncol') - goto 9999 - end if - - call psb_realloc(ncol,ilaggr,info) - if (info /= psb_success_) then - info=psb_err_from_subroutine_ - ch_err='psb_realloc' - call psb_errpush(info,name,a_err=ch_err) - goto 9999 - end if - - allocate(nlaggr(np),stat=info) - if (info /= psb_success_) then - info=psb_err_alloc_request_ - call psb_errpush(info,name,i_err=(/np,0,0,0,0/),& - & a_err='integer') - goto 9999 - end if - - nlaggr(:) = 0 - nlaggr(me+1) = naggr - call psb_sum(ictxt,nlaggr(1:np)) - - 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 - - end subroutine mld_dec_map_bld - end subroutine mld_zaggrmap_bld