*** empty log message ***

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,6 +107,163 @@ 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
+
+
+  if (.true.) then
+
+    !
+    ! New version. 
+    !
+
+    ! Note: -(nr+1)  Untouched as yet
+    !       -i    1<=i<=nr  Adjacent to aggregate i
+    !        i    1<=i<=nr  Belonging to aggregate i
+
+
+    !
+    ! Phase one: Start with disjoint groups.
+    ! 
+    naggr = 0
+    icnt = 0
+    step1: 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='csget')
+          goto 9999
+        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
+          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 
+          !
+          ! If the whole strongly coupled neighborhood of I is
+          ! as yet unconnected, turn it into the next aggregate
+          !
+          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: join the neighbours
+    ! 
+    step2: 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
+        !
+        ! 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
+
+  
+    !
+    ! 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
+        !
+        ! 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 (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 step3
+    
+
+  else
+    ! Original version, keep it for the time being 
+
     ! Note: -(nr+1)  Untouched as yet
     !       -i    1<=i<=nr  Adjacent to aggregate i
     !        i    1<=i<=nr  Belonging to aggregate i
@@ -274,6 +431,7 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
 
       end if
     end do
+
     if (debug_level >= psb_debug_outer_) then 
       if (recovery) then 
         write(debug_unit,*) me,' ',trim(name),&
@@ -292,6 +450,17 @@ subroutine mld_d_dec_map_bld(iorder,theta,a,desc_a,nlaggr,ilaggr,info)
       write(debug_unit,*) me,' ',trim(name),&
            & 'Leftovers ',count(ilaggr<0), '  in ',nlp,' loops'
     end if
+    deallocate(ils,neigh,stat=info)
+    if (info /= psb_success_) then 
+      info=psb_err_alloc_dealloc_
+      call psb_errpush(info,name)
+      goto 9999
+    end if
+
+
+  end if
+
+
 
   if (count(ilaggr<0) >0) then 
     info=psb_err_internal_error_
@@ -299,12 +468,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_