From f7df243b6fee121c18ef27ea6b42b2dd2d9c97e2 Mon Sep 17 00:00:00 2001 From: Salvatore Filippone Date: Fri, 19 Aug 2016 15:44:10 +0000 Subject: [PATCH] mld2p4-2: mlprec/impl/mld_c_dec_map_bld.f90 mlprec/impl/mld_d_dec_map_bld.f90 mlprec/impl/mld_s_dec_map_bld.f90 mlprec/impl/mld_z_dec_map_bld.f90 Reimplemented Vanek-Brezina decoupled aggregation. --- mlprec/impl/mld_c_dec_map_bld.f90 | 279 ++++++++++++------------------ mlprec/impl/mld_d_dec_map_bld.f90 | 277 ++++++++++++----------------- mlprec/impl/mld_s_dec_map_bld.f90 | 277 ++++++++++++----------------- mlprec/impl/mld_z_dec_map_bld.f90 | 279 ++++++++++++------------------ 4 files changed, 458 insertions(+), 654 deletions(-) diff --git a/mlprec/impl/mld_c_dec_map_bld.f90 b/mlprec/impl/mld_c_dec_map_bld.f90 index 7d0c5f94..35452b80 100644 --- a/mlprec/impl/mld_c_dec_map_bld.f90 +++ b/mlprec/impl/mld_c_dec_map_bld.f90 @@ -56,10 +56,10 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& & ideg(:), idxs(:) complex(psb_spk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_c_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl - logical :: recovery + logical :: recovery, candidate integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -107,191 +107,148 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call acsr%free() call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) end if - ! Note: -(nr+1) Untouched as yet - ! -i 1<=i<=nr Adjacent to aggregate i - ! i 1<=i<=nr Belonging to aggregate i + + ! - ! Phase one: group nodes together. - ! Very simple minded strategy. + ! Phase one: Start with disjoint groups. ! naggr = 0 - nlp = 0 - do - icnt = 0 - do ii=1, nr - i = idxs(ii) - 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 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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)) then - ilg = ilaggr(j) - if ((ilg<0).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif - end if + ! + ! Build the set of all strongly coupled nodes + ! + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ip = ip + 1 + icol(ip) = icol(k) end if - enddo + end if + enddo + if (ip < 1) then + write(0,*) "Should at least contain the node itself ! " + cycle step1 + end if + candidate = .true. + do k=1, ip + candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1)) + end do + if (candidate) then ! - ! 2. Untouched neighbours of these nodes are marked <0. + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate ! - call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_) - 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 + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, ip + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 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,izero,izero,izero,izero/),& - & 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 + ! Phase two: join the neighbours + ! + step2: do ii=1,nr + i = idxs(ii) - 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 = czero + if (ilaggr(i) == -(nr+1)) then 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = szero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then + ip = k + cpling = abs(val(k)) + end if end if enddo + if (ip > 0) then + ilaggr(i) = ilaggr(icol(ip)) + end if + end if + end do step2 - 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 + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) - if (ia>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 + if (ilaggr(i) == -(nr+1)) then + 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = szero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) < 0)) then + ip = ip + 1 + icol(ip) = icol(k) + end if end if - ils(ia) = ils(ia) + 1 - endif - + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + else + ! + ! This should not happen: we did not even connect with ourselves. + ! Create an isolate anyway. + ! + naggr = naggr + 1 + ilaggr(i) = naggr + end if 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 + end do step3 + if (count(ilaggr<0) >0) then info=psb_err_internal_error_ @@ -299,12 +256,6 @@ subroutine mld_c_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_ diff --git a/mlprec/impl/mld_d_dec_map_bld.f90 b/mlprec/impl/mld_d_dec_map_bld.f90 index 67e4437b..fe263c2a 100644 --- a/mlprec/impl/mld_d_dec_map_bld.f90 +++ b/mlprec/impl/mld_d_dec_map_bld.f90 @@ -56,10 +56,10 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& & ideg(:), idxs(:) real(psb_dpk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_d_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl - logical :: recovery + logical :: recovery, candidate integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -107,191 +107,148 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call acsr%free() call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) end if - ! Note: -(nr+1) Untouched as yet - ! -i 1<=i<=nr Adjacent to aggregate i - ! i 1<=i<=nr Belonging to aggregate i + + ! - ! Phase one: group nodes together. - ! Very simple minded strategy. + ! Phase one: Start with disjoint groups. ! naggr = 0 - nlp = 0 - do - icnt = 0 - do ii=1, nr - i = idxs(ii) - 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 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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)) then - ilg = ilaggr(j) - if ((ilg<0).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif - end if + ! + ! Build the set of all strongly coupled nodes + ! + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ip = ip + 1 + icol(ip) = icol(k) end if - enddo + end if + enddo + if (ip < 1) then + write(0,*) "Should at least contain the node itself ! " + cycle step1 + end if + candidate = .true. + do k=1, ip + candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1)) + end do + if (candidate) then ! - ! 2. Untouched neighbours of these nodes are marked <0. + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate ! - call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_) - 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 + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, ip + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 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,izero,izero,izero,izero/),& - & a_err='integer') - goto 9999 - end if + ! Phase two: join the neighbours + ! + step2: do ii=1,nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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 - - 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. + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate ! - isz = nr+1 - ia = -1 cpling = dzero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then + ip = k + cpling = abs(val(k)) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(icol(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = dzero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) < 0)) then + ip = ip + 1 + icol(ip) = icol(k) + end if 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 - + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do 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 - + ! + ! This should not happen: we did not even connect with ourselves. + ! Create an isolate anyway. + ! + naggr = naggr + 1 + ilaggr(i) = naggr + end if 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 + end do step3 + if (count(ilaggr<0) >0) then info=psb_err_internal_error_ @@ -299,12 +256,6 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_ diff --git a/mlprec/impl/mld_s_dec_map_bld.f90 b/mlprec/impl/mld_s_dec_map_bld.f90 index e980ecf8..10be26e0 100644 --- a/mlprec/impl/mld_s_dec_map_bld.f90 +++ b/mlprec/impl/mld_s_dec_map_bld.f90 @@ -56,10 +56,10 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& & ideg(:), idxs(:) real(psb_spk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_s_csr_sparse_mat) :: acsr real(psb_spk_) :: cpling, tcl - logical :: recovery + logical :: recovery, candidate integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -107,191 +107,148 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call acsr%free() call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) end if - ! Note: -(nr+1) Untouched as yet - ! -i 1<=i<=nr Adjacent to aggregate i - ! i 1<=i<=nr Belonging to aggregate i + + ! - ! Phase one: group nodes together. - ! Very simple minded strategy. + ! Phase one: Start with disjoint groups. ! naggr = 0 - nlp = 0 - do - icnt = 0 - do ii=1, nr - i = idxs(ii) - 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 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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)) then - ilg = ilaggr(j) - if ((ilg<0).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif - end if + ! + ! Build the set of all strongly coupled nodes + ! + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ip = ip + 1 + icol(ip) = icol(k) end if - enddo + end if + enddo + if (ip < 1) then + write(0,*) "Should at least contain the node itself ! " + cycle step1 + end if + candidate = .true. + do k=1, ip + candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1)) + end do + if (candidate) then ! - ! 2. Untouched neighbours of these nodes are marked <0. + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate ! - call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_) - 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 + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, ip + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 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,izero,izero,izero,izero/),& - & a_err='integer') - goto 9999 - end if + ! Phase two: join the neighbours + ! + step2: do ii=1,nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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 - - 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. + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate ! - isz = nr+1 - ia = -1 cpling = szero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then + ip = k + cpling = abs(val(k)) + end if + end if + enddo + if (ip > 0) then + ilaggr(i) = ilaggr(icol(ip)) + end if + end if + end do step2 + + + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) + + if (ilaggr(i) == -(nr+1)) then 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = szero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) < 0)) then + ip = ip + 1 + icol(ip) = icol(k) + end if 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 - + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do 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 - + ! + ! This should not happen: we did not even connect with ourselves. + ! Create an isolate anyway. + ! + naggr = naggr + 1 + ilaggr(i) = naggr + end if 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 + end do step3 + if (count(ilaggr<0) >0) then info=psb_err_internal_error_ @@ -299,12 +256,6 @@ subroutine mld_s_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_ diff --git a/mlprec/impl/mld_z_dec_map_bld.f90 b/mlprec/impl/mld_z_dec_map_bld.f90 index 718856aa..876b8a7e 100644 --- a/mlprec/impl/mld_z_dec_map_bld.f90 +++ b/mlprec/impl/mld_z_dec_map_bld.f90 @@ -56,10 +56,10 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) integer(psb_ipk_), allocatable :: ils(:), neigh(:), irow(:), icol(:),& & ideg(:), idxs(:) complex(psb_dpk_), allocatable :: val(:), diag(:) - integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii + integer(psb_ipk_) :: icnt,nlp,k,n,ia,isz,nr, naggr,i,j,m, nz, ilg, ii, ip type(psb_z_csr_sparse_mat) :: acsr real(psb_dpk_) :: cpling, tcl - logical :: recovery + logical :: recovery, candidate integer(psb_ipk_) :: debug_level, debug_unit,err_act integer(psb_ipk_) :: ictxt,np,me integer(psb_ipk_) :: nrow, ncol, n_ne @@ -107,191 +107,148 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) call acsr%free() call psb_msort(ideg,ix=idxs,dir=psb_sort_down_) end if - ! Note: -(nr+1) Untouched as yet - ! -i 1<=i<=nr Adjacent to aggregate i - ! i 1<=i<=nr Belonging to aggregate i + + ! - ! Phase one: group nodes together. - ! Very simple minded strategy. + ! Phase one: Start with disjoint groups. ! naggr = 0 - nlp = 0 - do - icnt = 0 - do ii=1, nr - i = idxs(ii) - 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 + icnt = 0 + step1: do ii=1, nr + i = idxs(ii) - 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 + if (ilaggr(i) == -(nr+1)) then + 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)) then - ilg = ilaggr(j) - if ((ilg<0).and.(i /= j)) then - if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then - ilaggr(j) = naggr - else - ilaggr(j) = -naggr - endif - end if + ! + ! Build the set of all strongly coupled nodes + ! + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if (abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j)))) then + ip = ip + 1 + icol(ip) = icol(k) end if - enddo + end if + enddo + if (ip < 1) then + write(0,*) "Should at least contain the node itself ! " + cycle step1 + end if + candidate = .true. + do k=1, ip + candidate = candidate .and. (ilaggr(icol(k)) == -(nr+1)) + end do + if (candidate) then ! - ! 2. Untouched neighbours of these nodes are marked <0. + ! If the whole strongly coupled neighborhood of I is + ! as yet unconnected, turn it into the next aggregate ! - call a%get_neigh(i,neigh,n_ne,info,lev=2_psb_ipk_) - 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 + icnt = icnt + 1 + naggr = naggr + 1 + do k=1, ip + ilaggr(icol(k)) = naggr + end do + ilaggr(i) = naggr + end if + endif + enddo step1 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,izero,izero,izero,izero/),& - & 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 + ! Phase two: join the neighbours + ! + step2: do ii=1,nr + i = idxs(ii) - 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 = zzero + if (ilaggr(i) == -(nr+1)) then 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = dzero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) > 0).and. (abs(val(k)) > cpling)) then + ip = k + cpling = abs(val(k)) + end if end if enddo + if (ip > 0) then + ilaggr(i) = ilaggr(icol(ip)) + end if + end if + end do step2 - 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 + ! + ! Phase three: sweep over leftovers, if any + ! + step3: do ii=1,nr + i = idxs(ii) - if (ia>naggr) then - info=psb_err_internal_error_ - call psb_errpush(info,name,a_err='loc_Aggregate: n > naggr ? ') - goto 9999 + if (ilaggr(i) == -(nr+1)) then + 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 + ! + ! Find the most strongly connected neighbour that is + ! already aggregated, if any, and join its aggregate + ! + cpling = dzero + ip = 0 + do k=1, nz + j = icol(k) + if ((1<=j).and.(j<=nr)) then + if ((abs(val(k)) > theta*sqrt(abs(diag(i)*diag(j))))& + & .and. (ilaggr(j) < 0)) then + ip = ip + 1 + icol(ip) = icol(k) + end if end if - ils(ia) = ils(ia) + 1 - endif - + enddo + if (ip > 0) then + icnt = icnt + 1 + naggr = naggr + 1 + ilaggr(i) = naggr + do k=1, ip + ilaggr(icol(k)) = naggr + end do + else + ! + ! This should not happen: we did not even connect with ourselves. + ! Create an isolate anyway. + ! + naggr = naggr + 1 + ilaggr(i) = naggr + end if 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 + end do step3 + if (count(ilaggr<0) >0) then info=psb_err_internal_error_ @@ -299,12 +256,6 @@ subroutine mld_z_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info) 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_