diff --git a/base/modules/psb_realloc_mod.F90 b/base/modules/psb_realloc_mod.F90
index a7f4039b..5e4fa670 100644
--- a/base/modules/psb_realloc_mod.F90
+++ b/base/modules/psb_realloc_mod.F90
@@ -456,7 +456,7 @@ Contains
   end function psb_zsize2d
 
 
-  Subroutine psb_icksz1d(len,v,info,pad)
+  Subroutine psb_icksz1d(len,v,info,pad,addsz,newsz)
     use psb_error_mod
 
     ! ...Subroutine Arguments  
@@ -464,6 +464,7 @@ Contains
     Integer,allocatable, intent(inout) :: v(:)
     integer         :: info
     integer, optional, intent(in) :: pad
+    integer, optional, intent(in) :: addsz,newsz
     ! ...Local Variables
     character(len=20)  :: name
     logical, parameter :: debug=.false.
@@ -476,7 +477,15 @@ Contains
     info=0
     
     If (len > psb_size(v)) Then
-      isz = max((3*psb_size(v))/2,(len+1))
+      if (present(newsz)) then 
+        isz = (max(len+1,newsz))
+      else
+        if (present(addsz)) then 
+          isz = len+max(1,addsz)
+        else
+          isz = len+1
+        endif
+      endif
       call psb_realloc(isz,v,info,pad=pad)
       
       if (info /= 0) then
@@ -503,13 +512,14 @@ Contains
   End Subroutine psb_icksz1d
 
 
-  Subroutine psb_dcksz1d(len,v,info,pad)
+  Subroutine psb_dcksz1d(len,v,info,pad,addsz,newsz)
     use psb_error_mod
 
     ! ...Subroutine Arguments  
     Integer,Intent(in)                 :: len
     real(kind(1.d0)),allocatable, intent(inout) :: v(:)
     integer         :: info
+    integer, optional, intent(in)          :: addsz,newsz
     real(kind(1.d0)), optional, intent(in) :: pad
     ! ...Local Variables
     character(len=20)  :: name
@@ -523,7 +533,16 @@ Contains
     info=0
     
     If (len > psb_size(v)) Then
-      isz = max((3*psb_size(v))/2,(len+1))
+      if (present(newsz)) then 
+        isz = (max(len+1,newsz))
+      else
+        if (present(addsz)) then 
+          isz = len+max(1,addsz)
+        else
+          isz = len+1
+        endif
+      endif
+
       call psb_realloc(isz,v,info,pad=pad)
       if (info /= 0) then
         info=4010
@@ -549,14 +568,15 @@ Contains
   End Subroutine psb_dcksz1d
 
 
-  Subroutine psb_zcksz1d(len,v,info,pad)
+  Subroutine psb_zcksz1d(len,v,info,pad,addsz,newsz)
     use psb_error_mod
 
     ! ...Subroutine Arguments  
-    Integer,Intent(in)                 :: len
+    Integer,Intent(in)                             :: len
     complex(kind(1.d0)),allocatable, intent(inout) :: v(:)
-    integer         :: info
-    complex(kind(1.d0)), optional, intent(in) :: pad
+    integer                                        :: info
+    integer, optional, intent(in)                  :: addsz,newsz
+    complex(kind(1.d0)), optional, intent(in)      :: pad
     ! ...Local Variables
     character(len=20)  :: name
     logical, parameter :: debug=.false.
@@ -569,7 +589,15 @@ Contains
     info=0
     
     If (len > psb_size(v)) Then
-      isz = max((3*psb_size(v))/2,(len+1))
+      if (present(newsz)) then 
+        isz = (max(len+1,newsz))
+      else
+        if (present(addsz)) then 
+          isz = len+max(1,addsz)
+        else
+          isz = len+1
+        endif
+      endif
       call psb_realloc(isz,v,info,pad=pad)
       if (info /= 0) then
         info=4010
diff --git a/base/modules/psb_sort_mod.f90 b/base/modules/psb_sort_mod.f90
index af292346..e0a15e10 100644
--- a/base/modules/psb_sort_mod.f90
+++ b/base/modules/psb_sort_mod.f90
@@ -42,6 +42,21 @@ module psb_sort_mod
     integer              :: last, dir
     integer, allocatable :: keys(:)
   end type psb_int_heap
+  type psb_int_idx_heap
+    integer              :: last, dir
+    integer, allocatable :: keys(:)
+    integer, allocatable :: idxs(:)
+  end type psb_int_idx_heap
+  type psb_double_idx_heap
+    integer              :: last, dir
+    real(kind(1.d0)), allocatable :: keys(:)
+    integer, allocatable          :: idxs(:)
+  end type psb_double_idx_heap
+  type psb_dcomplex_idx_heap
+    integer              :: last, dir
+    complex(kind(1.d0)), allocatable :: keys(:)
+    integer, allocatable             :: idxs(:)
+  end type psb_dcomplex_idx_heap
 
   interface psb_msort
     module procedure imsort
@@ -56,24 +71,28 @@ module psb_sort_mod
   end interface
 
   interface psb_init_heap
-    module procedure psb_init_int_heap
+    module procedure psb_init_int_heap, psb_init_int_idx_heap,&
+         & psb_init_double_idx_heap, psb_init_dcomplex_idx_heap
   end interface
 
   interface psb_dump_heap
-    module procedure psb_dump_int_heap
+    module procedure psb_dump_int_heap, psb_dump_int_idx_heap,&
+         & psb_dump_double_idx_heap, psb_dump_dcomplex_idx_heap
   end interface
 
   interface psb_howmany_heap
-    module procedure psb_howmany_int_heap
+    module procedure psb_howmany_int_heap, psb_howmany_int_idx_heap,&
+         & psb_howmany_double_idx_heap, psb_howmany_dcomplex_idx_heap
   end interface
 
-
   interface psb_insert_heap
-    module procedure psb_insert_int_heap
+    module procedure psb_insert_int_heap, psb_insert_int_idx_heap,&
+         & psb_insert_double_idx_heap, psb_insert_dcomplex_idx_heap
   end interface
 
   interface psb_heap_get_first
-    module procedure psb_int_heap_get_first
+    module procedure psb_int_heap_get_first, psb_int_idx_heap_get_first,&
+         & psb_double_idx_heap_get_first, psb_dcomplex_idx_heap_get_first
   end interface
 
 
@@ -442,7 +461,7 @@ contains
       heap%dir = psb_sort_up_
     endif
     select case(heap%dir) 
-    case (psb_sort_up_,psb_sort_down_)
+    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
       ! ok, do nothing
     case default
       write(0,*) 'Invalid direction, defaulting to psb_sort_up_'
@@ -485,7 +504,8 @@ contains
     integer, intent(in)               :: key
     type(psb_int_heap), intent(inout) :: heap
     integer, intent(out)              :: info
-    integer                           :: i, i2, itemp
+    integer                           :: i, i2
+    integer                           :: temp
     info = 0
     if (heap%last < 0) then 
       write(0,*) 'Invalid last in heap ',heap%last
@@ -493,25 +513,27 @@ contains
       return
     endif
 
+    heap%last = heap%last + 1
+    call psb_ensure_size(heap%last,heap%keys,info,addsz=psb_heap_resize)
+    if (info /= 0) then 
+      write(0,*) 'Memory allocation failure in heap_insert'
+      info = -5
+      return
+    end if
+
+    i            = heap%last
+    heap%keys(i) = key
+    
     select case(heap%dir)
     case (psb_sort_up_)
-      heap%last = heap%last + 1
-      if (heap%last > size(heap%keys)) then
-        call psb_realloc(heap%last+psb_heap_resize,heap%keys,info)
-        if (info /= 0) then
-          write(0,*) 'Memory allocation failure'
-          return
-        end if
-      end if
-      i = heap%last
-      heap%keys(i) = key
+
       do 
         if (i<=1) exit
         i2 = i/2
         if (heap%keys(i) < heap%keys(i2)) then 
-          itemp         = heap%keys(i)
+          temp          = heap%keys(i)
           heap%keys(i)  = heap%keys(i2)
-          heap%keys(i2) = itemp
+          heap%keys(i2) = temp
           i             = i2
         else
           exit
@@ -520,23 +542,45 @@ contains
       
           
     case (psb_sort_down_)
-      heap%last = heap%last + 1
-      if (heap%last > size(heap%keys)) then
-        call psb_realloc(heap%last+psb_heap_resize,heap%keys,info)
-        if (info /= 0) then
-          write(0,*) 'Memory allocation failure'
-          return
-        end if
-      end if
-      i = heap%last
-      heap%keys(i) = key
+
       do 
         if (i<=1) exit
         i2 = i/2
         if (heap%keys(i) > heap%keys(i2)) then 
-          itemp         = heap%keys(i)
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp
+          i             = i2
+        else
+          exit
+        end if
+      end do
+
+    case (psb_asort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) < abs(heap%keys(i2))) then 
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) > abs(heap%keys(i2))) then 
+          temp          = heap%keys(i)
           heap%keys(i)  = heap%keys(i2)
-          heap%keys(i2) = itemp
+          heap%keys(i2) = temp
           i             = i2
         else
           exit
@@ -558,7 +602,8 @@ contains
     type(psb_int_heap), intent(inout) :: heap
     integer, intent(out)              :: key,info
     
-    integer                           :: i, i2, itemp,last,j
+    integer                           :: i, i2, last,j
+    integer                           :: temp
 
     
     info = 0
@@ -568,10 +613,11 @@ contains
       return
     endif
 
-    key = heap%keys(1)
+    key          = heap%keys(1)
     heap%keys(1) = heap%keys(heap%last)
     heap%last    = heap%last - 1
     last         = heap%last
+
     select case(heap%dir)
     case (psb_sort_up_)
 
@@ -586,9 +632,9 @@ contains
         end if
 
         if (heap%keys(i) > heap%keys(j)) then 
-          itemp         = heap%keys(i)
+          temp          = heap%keys(i)
           heap%keys(i)  = heap%keys(j)
-          heap%keys(j)  = itemp
+          heap%keys(j)  = temp
           i             = j 
         else
           exit
@@ -609,9 +655,54 @@ contains
         end if
 
         if (heap%keys(i) < heap%keys(j)) then 
-          itemp         = heap%keys(i)
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case (psb_asort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) < abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) > abs(heap%keys(j))) then 
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) > abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+        
+        if (abs(heap%keys(i)) < abs(heap%keys(j))) then 
+          temp          = heap%keys(i)
           heap%keys(i)  = heap%keys(j)
-          heap%keys(j)  = itemp
+          heap%keys(j)  = temp
           i             = j 
         else
           exit
@@ -627,4 +718,944 @@ contains
 
 
 
+  function  psb_howmany_double_idx_heap(heap)
+    implicit none 
+    type(psb_double_idx_heap), intent(in) :: heap
+    integer :: psb_howmany_double_idx_heap
+    psb_howmany_double_idx_heap = heap%last
+  end function psb_howmany_double_idx_heap
+
+  subroutine psb_init_double_idx_heap(heap,info,dir)
+    use psb_realloc_mod
+    implicit none 
+    type(psb_double_idx_heap), intent(inout) :: heap
+    integer, intent(out)            :: info
+    integer, intent(in), optional   :: dir
+    
+    info = 0
+    heap%last=0
+    if (present(dir)) then 
+      heap%dir = dir
+    else
+      heap%dir = psb_sort_up_
+    endif
+    select case(heap%dir) 
+    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
+      ! ok, do nothing
+    case default
+      write(0,*) 'Invalid direction, defaulting to psb_sort_up_'
+      heap%dir = psb_sort_up_
+    end select
+    
+    call psb_ensure_size(psb_heap_resize,heap%keys,info)
+    call psb_ensure_size(psb_heap_resize,heap%idxs,info)
+    return
+  end subroutine psb_init_double_idx_heap
+
+  subroutine psb_dump_double_idx_heap(iout,heap,info)
+    implicit none 
+    type(psb_double_idx_heap), intent(in) :: heap
+    integer, intent(out)           :: info
+    integer, intent(in)            :: iout
+    
+    info = 0
+    if (iout < 0) then
+      write(0,*) 'Invalid file '
+      info =-1
+      return
+    end if
+    
+    write(iout,*) 'Heap direction ',heap%dir
+    write(iout,*) 'Heap size      ',heap%last
+    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
+         & (size(heap%keys)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else    if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
+         & (size(heap%idxs)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else
+      if (heap%last > 0) then 
+        write(iout,*) heap%keys(1:heap%last)
+        write(iout,*) heap%idxs(1:heap%last)
+      end if
+    end if
+  end subroutine psb_dump_double_idx_heap
+
+  subroutine psb_insert_double_idx_heap(key,index,heap,info)
+    use psb_realloc_mod
+    implicit none 
+
+    real(kind(1.d0)), intent(in)      :: key
+    integer, intent(in)               :: index
+    type(psb_double_idx_heap), intent(inout) :: heap
+    integer, intent(out)              :: info
+    integer                           :: i, i2, itemp
+    real(kind(1.d0))                  :: temp 
+    info = 0
+    if (heap%last < 0) then 
+      write(0,*) 'Invalid last in heap ',heap%last
+      info = heap%last
+      return
+    endif
+
+    heap%last = heap%last + 1
+    call psb_ensure_size(heap%last,heap%keys,info,addsz=psb_heap_resize)
+    if (info == 0) &
+         & call psb_ensure_size(heap%last,heap%idxs,info,addsz=psb_heap_resize)
+    if (info /= 0) then 
+      write(0,*) 'Memory allocation failure in heap_insert'
+      info = -5
+      return
+    end if
+    
+    i            = heap%last
+    heap%keys(i) = key
+    heap%idxs(i) = index
+
+    select case(heap%dir)
+    case (psb_sort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (heap%keys(i) < heap%keys(i2)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_sort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (heap%keys(i) > heap%keys(i2)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+    case (psb_asort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) < abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) > abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_insert_double_idx_heap
+
+  subroutine psb_double_idx_heap_get_first(key,index,heap,info)
+    implicit none 
+
+    type(psb_double_idx_heap), intent(inout) :: heap
+    integer, intent(out)              :: index,info
+    real(kind(1.d0)), intent(out)     :: key
+    
+    integer                           :: i, i2, last,j,itemp
+    real(kind(1.d0))                  :: temp
+
+    
+    info = 0
+    if (heap%last <= 0) then 
+      key  = 0
+      info = -1
+      return
+    endif
+
+    key          = heap%keys(1)
+    index        = heap%idxs(1)
+    heap%keys(1) = heap%keys(heap%last)
+    heap%idxs(1) = heap%idxs(heap%last)
+    heap%last    = heap%last - 1
+    last         = heap%last
+
+    select case(heap%dir)
+    case (psb_sort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (heap%keys(2*i) < heap%keys(2*i+1)) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (heap%keys(i) > heap%keys(j)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_sort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (heap%keys(2*i) > heap%keys(2*i+1)) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (heap%keys(i) < heap%keys(j)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case (psb_asort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) < abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) > abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) > abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) < abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_double_idx_heap_get_first
+
+  function  psb_howmany_int_idx_heap(heap)
+    implicit none 
+    type(psb_int_idx_heap), intent(in) :: heap
+    integer :: psb_howmany_int_idx_heap
+    psb_howmany_int_idx_heap = heap%last
+  end function psb_howmany_int_idx_heap
+
+  subroutine psb_init_int_idx_heap(heap,info,dir)
+    use psb_realloc_mod
+    implicit none 
+    type(psb_int_idx_heap), intent(inout) :: heap
+    integer, intent(out)            :: info
+    integer, intent(in), optional   :: dir
+    
+    info = 0
+    heap%last=0
+    if (present(dir)) then 
+      heap%dir = dir
+    else
+      heap%dir = psb_sort_up_
+    endif
+    select case(heap%dir) 
+    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
+      ! ok, do nothing
+    case default
+      write(0,*) 'Invalid direction, defaulting to psb_sort_up_'
+      heap%dir = psb_sort_up_
+    end select
+    
+    call psb_ensure_size(psb_heap_resize,heap%keys,info)
+    call psb_ensure_size(psb_heap_resize,heap%idxs,info)
+    return
+  end subroutine psb_init_int_idx_heap
+
+  subroutine psb_dump_int_idx_heap(iout,heap,info)
+    implicit none 
+    type(psb_int_idx_heap), intent(in) :: heap
+    integer, intent(out)           :: info
+    integer, intent(in)            :: iout
+    
+    info = 0
+    if (iout < 0) then
+      write(0,*) 'Invalid file '
+      info =-1
+      return
+    end if
+    
+    write(iout,*) 'Heap direction ',heap%dir
+    write(iout,*) 'Heap size      ',heap%last
+    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
+         & (size(heap%keys)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else    if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
+         & (size(heap%idxs)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else
+      if (heap%last > 0) then 
+        write(iout,*) heap%keys(1:heap%last)
+        write(iout,*) heap%idxs(1:heap%last)
+      end if
+    end if
+  end subroutine psb_dump_int_idx_heap
+
+  subroutine psb_insert_int_idx_heap(key,index,heap,info)
+    use psb_realloc_mod
+    implicit none 
+
+    integer, intent(in)                   :: key
+    integer, intent(in)                   :: index
+    type(psb_int_idx_heap), intent(inout) :: heap
+    integer, intent(out)                  :: info
+    integer                               :: i, i2, itemp
+    integer                               :: temp 
+    info = 0
+    if (heap%last < 0) then 
+      write(0,*) 'Invalid last in heap ',heap%last
+      info = heap%last
+      return
+    endif
+
+    heap%last = heap%last + 1
+    call psb_ensure_size(heap%last,heap%keys,info,addsz=psb_heap_resize)
+    if (info == 0) &
+         & call psb_ensure_size(heap%last,heap%idxs,info,addsz=psb_heap_resize)
+    if (info /= 0) then 
+      write(0,*) 'Memory allocation failure in heap_insert'
+      info = -5
+      return
+    end if
+    
+    i            = heap%last
+    heap%keys(i) = key
+    heap%idxs(i) = index
+
+    select case(heap%dir)
+    case (psb_sort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (heap%keys(i) < heap%keys(i2)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_sort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (heap%keys(i) > heap%keys(i2)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+    case (psb_asort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) < abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) > abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_insert_int_idx_heap
+
+  subroutine psb_int_idx_heap_get_first(key,index,heap,info)
+    implicit none 
+
+    type(psb_int_idx_heap), intent(inout) :: heap
+    integer, intent(out)                  :: index,info
+    integer, intent(out)                  :: key
+    
+    integer                               :: i, i2, last,j,itemp
+    integer                               :: temp
+
+    
+    info = 0
+    if (heap%last <= 0) then 
+      key  = 0
+      info = -1
+      return
+    endif
+
+    key          = heap%keys(1)
+    index        = heap%idxs(1)
+    heap%keys(1) = heap%keys(heap%last)
+    heap%idxs(1) = heap%idxs(heap%last)
+    heap%last    = heap%last - 1
+    last         = heap%last
+
+    select case(heap%dir)
+    case (psb_sort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (heap%keys(2*i) < heap%keys(2*i+1)) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (heap%keys(i) > heap%keys(j)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_sort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (heap%keys(2*i) > heap%keys(2*i+1)) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (heap%keys(i) < heap%keys(j)) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case (psb_asort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) < abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) > abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) > abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) < abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_int_idx_heap_get_first
+
+
+
+  function  psb_howmany_dcomplex_idx_heap(heap)
+    implicit none 
+    type(psb_dcomplex_idx_heap), intent(in) :: heap
+    integer :: psb_howmany_dcomplex_idx_heap
+    psb_howmany_dcomplex_idx_heap = heap%last
+  end function psb_howmany_dcomplex_idx_heap
+
+  subroutine psb_init_dcomplex_idx_heap(heap,info,dir)
+    use psb_realloc_mod
+    implicit none 
+    type(psb_dcomplex_idx_heap), intent(inout) :: heap
+    integer, intent(out)            :: info
+    integer, intent(in), optional   :: dir
+    
+    info = 0
+    heap%last=0
+    if (present(dir)) then 
+      heap%dir = dir
+    else
+      heap%dir = psb_sort_up_
+    endif
+    select case(heap%dir) 
+!!$    case (psb_sort_up_,psb_sort_down_,psb_asort_up_,psb_asort_down_)
+    case (psb_asort_up_,psb_asort_down_)
+      ! ok, do nothing
+    case default
+      write(0,*) 'Invalid direction, defaulting to psb_sort_up_'
+      heap%dir = psb_asort_up_
+    end select
+    
+    call psb_ensure_size(psb_heap_resize,heap%keys,info)
+    call psb_ensure_size(psb_heap_resize,heap%idxs,info)
+    return
+  end subroutine psb_init_dcomplex_idx_heap
+
+  subroutine psb_dump_dcomplex_idx_heap(iout,heap,info)
+    implicit none 
+    type(psb_dcomplex_idx_heap), intent(in) :: heap
+    integer, intent(out)           :: info
+    integer, intent(in)            :: iout
+    
+    info = 0
+    if (iout < 0) then
+      write(0,*) 'Invalid file '
+      info =-1
+      return
+    end if
+    
+    write(iout,*) 'Heap direction ',heap%dir
+    write(iout,*) 'Heap size      ',heap%last
+    if ((heap%last > 0).and.((.not.allocated(heap%keys)).or.&
+         & (size(heap%keys)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else    if ((heap%last > 0).and.((.not.allocated(heap%idxs)).or.&
+         & (size(heap%idxs)<heap%last))) then
+      write(iout,*) 'Inconsistent size/allocation status!!'
+    else
+      if (heap%last > 0) then 
+        write(iout,*) heap%keys(1:heap%last)
+        write(iout,*) heap%idxs(1:heap%last)
+      end if
+    end if
+  end subroutine psb_dump_dcomplex_idx_heap
+
+  subroutine psb_insert_dcomplex_idx_heap(key,index,heap,info)
+    use psb_realloc_mod
+    implicit none 
+
+    complex(kind(1.d0)), intent(in)            :: key
+    integer, intent(in)                        :: index
+    type(psb_dcomplex_idx_heap), intent(inout) :: heap
+    integer, intent(out)                       :: info
+    integer                                    :: i, i2, itemp
+    complex(kind(1.d0))                        :: temp 
+    info = 0
+    if (heap%last < 0) then 
+      write(0,*) 'Invalid last in heap ',heap%last
+      info = heap%last
+      return
+    endif
+
+    heap%last = heap%last + 1
+    call psb_ensure_size(heap%last,heap%keys,info,addsz=psb_heap_resize)
+    if (info == 0) &
+         & call psb_ensure_size(heap%last,heap%idxs,info,addsz=psb_heap_resize)
+    if (info /= 0) then 
+      write(0,*) 'Memory allocation failure in heap_insert'
+      info = -5
+      return
+    end if
+    
+    i            = heap%last
+    heap%keys(i) = key
+    heap%idxs(i) = index
+
+    select case(heap%dir)
+!!$    case (psb_sort_up_)
+!!$
+!!$      do 
+!!$        if (i<=1) exit
+!!$        i2 = i/2
+!!$        if (heap%keys(i) < heap%keys(i2)) then 
+!!$          itemp         = heap%idxs(i)
+!!$          heap%idxs(i)  = heap%idxs(i2)
+!!$          heap%idxs(i2) = itemp
+!!$          temp          = heap%keys(i)
+!!$          heap%keys(i)  = heap%keys(i2)
+!!$          heap%keys(i2) = temp 
+!!$          i             = i2
+!!$        else
+!!$          exit
+!!$        end if
+!!$      end do
+!!$      
+!!$          
+!!$    case (psb_sort_down_)
+!!$
+!!$      do 
+!!$        if (i<=1) exit
+!!$        i2 = i/2
+!!$        if (heap%keys(i) > heap%keys(i2)) then 
+!!$          itemp         = heap%idxs(i)
+!!$          heap%idxs(i)  = heap%idxs(i2)
+!!$          heap%idxs(i2) = itemp
+!!$          temp          = heap%keys(i)
+!!$          heap%keys(i)  = heap%keys(i2)
+!!$          heap%keys(i2) = temp 
+!!$          i             = i2
+!!$        else
+!!$          exit
+!!$        end if
+!!$      end do
+      
+    case (psb_asort_up_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) < abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      do 
+        if (i<=1) exit
+        i2 = i/2
+        if (abs(heap%keys(i)) > abs(heap%keys(i2))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(i2)
+          heap%idxs(i2) = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(i2)
+          heap%keys(i2) = temp 
+          i             = i2
+        else
+          exit
+        end if
+      end do
+      
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_insert_dcomplex_idx_heap
+
+  subroutine psb_dcomplex_idx_heap_get_first(key,index,heap,info)
+    implicit none 
+
+    type(psb_dcomplex_idx_heap), intent(inout) :: heap
+    integer, intent(out)                       :: index,info
+    complex(kind(1.d0)), intent(out)           :: key
+    
+    integer                                    :: i, i2, last,j,itemp
+    complex(kind(1.d0))                        :: temp
+
+    
+    info = 0
+    if (heap%last <= 0) then 
+      key  = 0
+      info = -1
+      return
+    endif
+
+    key          = heap%keys(1)
+    index        = heap%idxs(1)
+    heap%keys(1) = heap%keys(heap%last)
+    heap%idxs(1) = heap%idxs(heap%last)
+    heap%last    = heap%last - 1
+    last         = heap%last
+
+    select case(heap%dir)
+!!$    case (psb_sort_up_)
+!!$
+!!$      i = 1
+!!$      do 
+!!$        if (i > (last/2)) exit
+!!$        if ( (heap%keys(2*i) < heap%keys(2*i+1)) .or.&
+!!$             & (2*i == last)) then 
+!!$          j = 2*i
+!!$        else
+!!$          j = 2*i + 1
+!!$        end if
+!!$
+!!$        if (heap%keys(i) > heap%keys(j)) then 
+!!$          itemp         = heap%idxs(i)
+!!$          heap%idxs(i)  = heap%idxs(j)
+!!$          heap%idxs(j)  = itemp
+!!$          temp          = heap%keys(i)
+!!$          heap%keys(i)  = heap%keys(j)
+!!$          heap%keys(j)  = temp 
+!!$          i             = j 
+!!$        else
+!!$          exit
+!!$        end if
+!!$      end do
+!!$      
+!!$          
+!!$    case (psb_sort_down_)
+!!$
+!!$      i = 1
+!!$      do 
+!!$        if (i > (last/2)) exit
+!!$        if ( (heap%keys(2*i) > heap%keys(2*i+1)) .or.&
+!!$             & (2*i == last)) then 
+!!$          j = 2*i
+!!$        else
+!!$          j = 2*i + 1
+!!$        end if
+!!$
+!!$        if (heap%keys(i) < heap%keys(j)) then 
+!!$          itemp         = heap%idxs(i)
+!!$          heap%idxs(i)  = heap%idxs(j)
+!!$          heap%idxs(j)  = itemp
+!!$          temp          = heap%keys(i)
+!!$          heap%keys(i)  = heap%keys(j)
+!!$          heap%keys(j)  = temp 
+!!$          i             = j 
+!!$        else
+!!$          exit
+!!$        end if
+!!$      end do
+
+    case (psb_asort_up_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) < abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) > abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+      
+          
+    case (psb_asort_down_)
+
+      i = 1
+      do 
+        if (i > (last/2)) exit
+        if ( (abs(heap%keys(2*i)) > abs(heap%keys(2*i+1))) .or.&
+             & (2*i == last)) then 
+          j = 2*i
+        else
+          j = 2*i + 1
+        end if
+
+        if (abs(heap%keys(i)) < abs(heap%keys(j))) then 
+          itemp         = heap%idxs(i)
+          heap%idxs(i)  = heap%idxs(j)
+          heap%idxs(j)  = itemp
+          temp          = heap%keys(i)
+          heap%keys(i)  = heap%keys(j)
+          heap%keys(j)  = temp 
+          i             = j 
+        else
+          exit
+        end if
+      end do
+
+    case default
+      write(0,*) 'Invalid direction in heap ',heap%dir
+    end select
+
+    return
+  end subroutine psb_dcomplex_idx_heap_get_first
+
+
 end module psb_sort_mod