diff --git a/base/internals/psi_bld_glb_dep_list.F90 b/base/internals/psi_bld_glb_dep_list.F90
index fb4f9f32..bf28e49b 100644
--- a/base/internals/psi_bld_glb_dep_list.F90
+++ b/base/internals/psi_bld_glb_dep_list.F90
@@ -181,7 +181,7 @@ subroutine psi_i_bld_glb_csr_dep_list(ictxt,loc_dl,length_dl,c_dep_list,dl_ptr,i
     goto 9999      
   end if
   icomm = psb_get_mpi_comm(iictxt)
-  call mpi_allgather(loc_dl,myld,psb_mpi_ipk_,&
+  call mpi_allgatherv(loc_dl,myld,psb_mpi_ipk_,&
        & c_dep_list,length_dl,dl_ptr,psb_mpi_ipk_,icomm,minfo)
   
   info = minfo  
diff --git a/base/internals/psi_crea_index.f90 b/base/internals/psi_crea_index.f90
index eb22893e..68c4615b 100644
--- a/base/internals/psi_crea_index.f90
+++ b/base/internals/psi_crea_index.f90
@@ -154,8 +154,8 @@ subroutine psi_i_crea_index(desc_a,index_in,index_out,nxch,nsnd,nrcv,info)
 !!$         & '  avg:',dlavg, choose_sorting(dlmax,dlavg,np)
 
     if (choose_sorting(dlmax,dlavg,np)) then 
-      if (.true.) then 
-        call psi_bld_glb_dep_list(ictxt,&
+      if (.false.) then 
+        call psi_bld_glb_dep_list(ictxt,& 
              & loc_dl,length_dl,dep_list,dl_lda,info)
 
         if (info /= psb_success_) then
@@ -196,6 +196,9 @@ subroutine psi_i_crea_index(desc_a,index_in,index_out,nxch,nsnd,nrcv,info)
 !!$
 !!$        ! ....now i can sort dependency lists.
         call psi_sort_dl(dl_ptr,c_dep_list,length_dl,np,info)
+        ldl    = length_dl(me)
+        loc_dl = c_dep_list(dl_ptr(me):dl_ptr(me)+ldl-1)
+        
 !!$        if(info /= psb_success_) then
 !!$          call psb_errpush(psb_err_from_subroutine_,name,a_err='psi_sort_dl')
 !!$          goto 9999
@@ -215,6 +218,7 @@ subroutine psi_i_crea_index(desc_a,index_in,index_out,nxch,nsnd,nrcv,info)
   if (do_timings) call psb_tic(idx_phase3)
   if(debug_level >= psb_debug_inner_)&
        & write(debug_unit,*) me,' ',trim(name),': calling psi_desc_index',ldl,':',loc_dl(1:ldl)
+
   ! Do the actual format conversion. 
   call psi_desc_index(desc_a,index_in,loc_dl,ldl,nsnd,nrcv,index_out,info)
   if(debug_level >= psb_debug_inner_) &
@@ -249,6 +253,7 @@ contains
     logical                       :: val
 
     val = .not.(((dlmax>(26*4)).or.((dlavg>=(26*2)).and.(np>=128))))
+    val = .true.
   end function choose_sorting
 
 end subroutine psi_i_crea_index
diff --git a/base/internals/psi_sort_dl.f90 b/base/internals/psi_sort_dl.f90
index 60212d30..29a55b49 100644
--- a/base/internals/psi_sort_dl.f90
+++ b/base/internals/psi_sort_dl.f90
@@ -69,7 +69,6 @@ subroutine psi_i_sort_dl(dep_list,l_dep_list,np,info)
   isz = iich + ndgmx
   if (debug_level >= psb_debug_inner_)&
        & write(debug_unit,*) name,': ndgmx ',ndgmx,isz
-
   allocate(work(isz))
   ! call srtlist(dep_list, dl_lda, l_dep_list, np, info)
   call srtlist(dep_list,size(dep_list,1,kind=psb_ipk_),l_dep_list,np,work(idg),&
@@ -89,19 +88,193 @@ subroutine psi_i_sort_dl(dep_list,l_dep_list,np,info)
 
 end subroutine psi_i_sort_dl
 
+!**********************************************************************
+!                                                                     *
+!               The communication step among processors at each       *
+!               matrix-vector  product is a variable all-to-all       *
+!               collective communication that we reimplement          *
+!               in terms of point-to-point communications.            *
+!               The data in input is a list of dependencies:          *
+!               for each node a list of all the nodes it has to       *
+!               communicate with. The lists are guaranteed to be      *
+!               symmetric, i.e. for each pair (I,J) there is a        *
+!               pair (J,I). The idea is to organize the ordering      *
+!               so   that at each communication step as many          *
+!               processors as possible are communicating at the       *
+!               same time, i.e. a step is defined by the fact         *
+!               that all edges (I,J) in it have no common node.       *
+!                                                                     *
+!               Formulation of the problem is:                        *
+!                 Given an undirected graph (forest):                 *
+!                 Find the shortest series of steps to cancel all     *
+!                 graph edges, where at each step all edges belonging *
+!                 to a matching in the graph are canceled.            *
+!                                                                     *
+!               An obvious lower bound to the optimum number of steps *
+!               is the largest degree of any node in the graph.       *
+!                                                                     *
+!               The algorithm proceeds as follows:                    *
+!               1. Build a list of all edges, e.g. copy the           *
+!                  dependencies lists keeping only (I,J) with I<J     *
+!               2. Compute an auxiliary vector with the degree of     *
+!                  each node of the graph.                            *
+!               3. While there are edges in the graph do:             *
+!               4. Weight the edges with the sum of the degrees       *
+!                  of their nodes and sort them into descending order *
+!               5. Scan the list of edges; if neither node of the     *
+!                  edge has been marked yet, cancel the edge and mark *
+!                  the two nodes                                      *
+!               6. If no edge was chosen but the graph is nonempty    *
+!                  raise an error condition                           *
+!               7. Queue the edges in the matchin to the output       *
+!                  sequence;                                          *
+!               8. Decrease by 1 the degree of all marked nodes,      *
+!                  then clear all marks                               *
+!               9. Cycle to 3.                                        *
+!              10. For each node: scan the edge sequence; if an       *
+!                  edge has the node as an endpoint, queue the other  *
+!                  node in the dependency list for the current one    *
+!                                                                     *
+!**********************************************************************
 subroutine psi_i_csr_sort_dl(dl_ptr,c_dep_list,l_dep_list,np,info)
   use psi_mod, psb_protect_name => psi_i_csr_sort_dl
   use psb_const_mod
   use psb_error_mod
+  use psb_sort_mod
   implicit none
   
-  integer(psb_ipk_), intent(in)    :: c_dep_list(:), dl_ptr(0:)
-  integer(psb_ipk_), intent(inout) :: l_dep_list(0:)
+  integer(psb_ipk_), intent(inout) :: c_dep_list(:), dl_ptr(0:), l_dep_list(0:)
   integer(psb_ipk_), intent(in)    :: np
   integer(psb_ipk_), intent(out)   :: info
   ! Local variables
+  integer(psb_ipk_), allocatable  ::  dg(:), dgp(:),&
+       &  idx(:), upd(:), edges(:,:), ich(:)
+  integer(psb_ipk_) ::  i, j, nedges, ip1, ip2, nch, ip, iedge,&
+       &  i1, ix, ist, iswap(2)
 
+  nedges = size(c_dep_list)
   
-end subroutine psi_i_csr_sort_dl
+  allocate(dg(0:np-1),dgp(nedges),edges(2,nedges),upd(0:np-1),&
+       & idx(nedges),ich(nedges),stat = info)
+  
+  if (info /= 0) then
+    info = -9
+    return
+  end if
+  !
+  !  1. Compute an auxiliary vector with the degree of 
+  !     each node of the graph.
+  dg(0:np-1) = l_dep_list(0:np-1)
+  !
+  !  2. Build a list of all edges, e.g. copy the      
+  !     dependencies lists keeping only (I,J) with I<J
+  !
+  nedges  = 0
+  do i = 0, np-1
+    do j = dl_ptr(i),dl_ptr(i+1) - 1
+      ip = c_dep_list(j)
+      if (i<ip) then
+        nedges = nedges + 1
+        edges(1,nedges) = i
+        edges(2,nedges) = ip
+      end if
+    end do
+  end do
 
-      
+  !
+  !   3. Loop over all edges
+  !   
+  ist = 1 
+  do while (ist <= nedges)
+    !
+    !  4. Weight the edges with the sum of the degrees       
+    !     of their nodes and sort them into descending order 
+    upd(:) = 0
+    do i = ist, nedges
+      dgp(i) = (dg(edges(1,i)) + dg(edges(2,i)))
+    end do
+    call psb_msort(dgp(ist:nedges),ix=idx(ist:nedges),dir=psb_sort_down_)
+    
+    !  5. Scan the list of edges; if neither node of the     
+    !     edge has been marked yet, take out the edge and mark 
+    !     the two nodes                                      
+    i1 = ist         
+    nch = 0
+    do i = ist, nedges
+      ix  = idx(i)+ist-1
+      ip1 = edges(1,ix)
+      ip2 = edges(2,ix)
+      if ((upd(ip1)==0).and.(upd(ip2)==0)) then 
+        upd(ip1) = -1
+        upd(ip2) = -1
+        nch      = nch + 1
+        ich(nch) = ix
+      end if
+    end do
+    !
+    !  6. If no edge was chosen but the graph is nonempty 
+    !     raise an error condition                        
+    if (nch == 0) then
+      write(psb_err_unit,*)&
+           & 'srtlist ?????? impossible error !!!!!?????',&
+           &  nedges,ist
+      do i=ist, nedges
+        ix = idx(i)+ist-1
+        write(psb_err_unit,*)&
+             & 'SRTLIST: Edge:',ix,edges(1,ix),&
+             & edges(2,ix),dgp(ix)
+      end do
+      info = psb_err_input_value_invalid_i_
+      return
+    end if
+    !
+    ! 7. Queue the edges in the matching to the output  
+    !    sequence; decrease by 1 the degree of all marked
+    !    nodes, then clear all marks
+    !    
+    call psb_msort(ich(1:nch))
+    do i=1, nch
+      iswap(1:2)        = edges(1:2,ist)
+      edges(1:2,ist)    = edges(1:2,ich(i))
+      edges(1:2,ich(i)) = iswap(1:2)
+      ist = ist + 1 
+    end do
+    do i=0, np-1
+      dg(i) = dg(i) + upd(i) 
+    end do
+  end do
+
+  do i=0, np-1
+    if (dg(i) /= 0) then 
+      write(psb_err_unit,*)&
+           &  'SRTLIST Error on exit:',i,dg(i)
+    end if
+    dg(i) = 0
+  end do
+  
+  !
+  ! 10. Scan the edge sequence;
+  !     for each edge, take each one of its
+  !     endpoints and queue the other
+  !     node in the endpoint dependency list
+  !     
+  do j=1,nedges
+    i     = edges(1,j)
+    ix    = dl_ptr(i)
+    c_dep_list(ix+dg(i)) = edges(2,j)
+    dg(i) = dg(i)+1
+    
+    i     = edges(2,j)
+    ix    = dl_ptr(i)
+    c_dep_list(ix+dg(i)) = edges(1,j)
+    dg(i) = dg(i)+1
+  end do
+
+  do i=0, np-1
+    if (dg(i) /= l_dep_list(i)) then 
+      write(psb_err_unit,*) &
+           &   'SRTLIST Mismatch on output',i,dg(i),l_dep_list(i)
+    end if
+  end do
+  
+end subroutine psi_i_csr_sort_dl
diff --git a/base/internals/psi_srtlist.f90 b/base/internals/psi_srtlist.f90
index fb9e2b7d..47a0e2b4 100644
--- a/base/internals/psi_srtlist.f90
+++ b/base/internals/psi_srtlist.f90
@@ -124,7 +124,6 @@ subroutine srtlist(dep_list,dl_lda,ldl,np,dg,dgp,upd, edges,idx,ich,info)
 
   ist = 1
   do while (ist.le.nedges)         
-
     do i=1, np
       upd(i) = 0      
     enddo