psblas3:

base/modules/psb_c_base_vect_mod.f90
 base/modules/psb_c_vect_mod.f90
 base/modules/psb_d_base_vect_mod.f90
 base/modules/psb_d_vect_mod.f90
 base/modules/psb_s_base_vect_mod.f90
 base/modules/psb_s_vect_mod.f90
 base/modules/psb_z_base_vect_mod.f90
 base/modules/psb_z_vect_mod.f90
 util/psb_c_renum_impl.F90
 util/psb_renum_mod.f90
 util/psb_s_renum_impl.F90
 util/psb_z_renum_impl.F90

Mods for mld-ainv.
AMD renumbering for all data types.

base/modules/psb_c_base_vect_mod.f90

@@ -7,6 +7,7 @@ module psb_c_base_vect_mod
     complex(psb_spk_), allocatable :: v(:)
   contains
     procedure, pass(x) :: get_nrows => c_base_get_nrows
+    procedure, pass(x) :: sizeof   => c_base_sizeof
     procedure, pass(x) :: dot_v    => c_base_dot_v
     procedure, pass(x) :: dot_a    => c_base_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -145,15 +146,21 @@ contains
 
   end function size_const
     
-
   function c_base_get_nrows(x) result(res)
     implicit none 
     class(psb_c_base_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = size(x%v)
   end function c_base_get_nrows
 
+  function c_base_sizeof(x) result(res)
+    implicit none 
+    class(psb_c_base_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = (2*psb_sizeof_sp)*x%get_nrows()
+  end function c_base_sizeof
+
   function c_base_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_c_base_vect_type), intent(inout) :: x, y

base/modules/psb_c_vect_mod.f90

@@ -6,6 +6,7 @@ module psb_c_vect_mod
     class(psb_c_base_vect_type), allocatable :: v 
   contains
     procedure, pass(x) :: get_nrows => c_vect_get_nrows
+    procedure, pass(x) :: sizeof   => c_vect_sizeof
     procedure, pass(x) :: dot_v    => c_vect_dot_v
     procedure, pass(x) :: dot_a    => c_vect_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -148,15 +149,22 @@ contains
 
   end function size_const
 
-
   function c_vect_get_nrows(x) result(res)
     implicit none 
     class(psb_c_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = x%v%get_nrows()
   end function c_vect_get_nrows
 
+  function c_vect_sizeof(x) result(res)
+    implicit none 
+    class(psb_c_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = 0
+    if (allocated(x%v)) res = x%v%sizeof()
+  end function c_vect_sizeof
+
   function c_vect_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_c_vect_type), intent(inout) :: x, y

base/modules/psb_d_base_vect_mod.f90

@@ -7,6 +7,7 @@ module psb_d_base_vect_mod
     real(psb_dpk_), allocatable :: v(:)
   contains
     procedure, pass(x) :: get_nrows => d_base_get_nrows
+    procedure, pass(x) :: sizeof   => d_base_sizeof
     procedure, pass(x) :: dot_v    => d_base_dot_v
     procedure, pass(x) :: dot_a    => d_base_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -145,15 +146,21 @@ contains
 
   end function size_const
     
-
   function d_base_get_nrows(x) result(res)
     implicit none 
     class(psb_d_base_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = size(x%v)
   end function d_base_get_nrows
 
+  function d_base_sizeof(x) result(res)
+    implicit none 
+    class(psb_d_base_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = psb_sizeof_dp*x%get_nrows()
+  end function d_base_sizeof
+
   function d_base_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_d_base_vect_type), intent(inout) :: x, y

base/modules/psb_d_vect_mod.f90

@@ -6,6 +6,7 @@ module psb_d_vect_mod
     class(psb_d_base_vect_type), allocatable :: v 
   contains
     procedure, pass(x) :: get_nrows => d_vect_get_nrows
+    procedure, pass(x) :: sizeof   => d_vect_sizeof
     procedure, pass(x) :: dot_v    => d_vect_dot_v
     procedure, pass(x) :: dot_a    => d_vect_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -148,15 +149,22 @@ contains
 
   end function size_const
 
-
   function d_vect_get_nrows(x) result(res)
     implicit none 
     class(psb_d_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = x%v%get_nrows()
   end function d_vect_get_nrows
 
+  function d_vect_sizeof(x) result(res)
+    implicit none 
+    class(psb_d_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = 0
+    if (allocated(x%v)) res = x%v%sizeof()
+  end function d_vect_sizeof
+
   function d_vect_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_d_vect_type), intent(inout) :: x, y

base/modules/psb_s_base_vect_mod.f90

@@ -7,6 +7,7 @@ module psb_s_base_vect_mod
     real(psb_spk_), allocatable :: v(:)
   contains
     procedure, pass(x) :: get_nrows => s_base_get_nrows
+    procedure, pass(x) :: sizeof   => s_base_sizeof
     procedure, pass(x) :: dot_v    => s_base_dot_v
     procedure, pass(x) :: dot_a    => s_base_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -145,15 +146,21 @@ contains
 
   end function size_const
     
-
   function s_base_get_nrows(x) result(res)
     implicit none 
     class(psb_s_base_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = size(x%v)
   end function s_base_get_nrows
     
+  function s_base_sizeof(x) result(res)
+    implicit none 
+    class(psb_s_base_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = psb_sizeof_sp*x%get_nrows()
+  end function s_base_sizeof
+
   function s_base_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_s_base_vect_type), intent(inout) :: x, y

base/modules/psb_s_vect_mod.f90

@@ -6,6 +6,7 @@ module psb_s_vect_mod
     class(psb_s_base_vect_type), allocatable :: v 
   contains
     procedure, pass(x) :: get_nrows => s_vect_get_nrows
+    procedure, pass(x) :: sizeof   => s_vect_sizeof
     procedure, pass(x) :: dot_v    => s_vect_dot_v
     procedure, pass(x) :: dot_a    => s_vect_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -148,15 +149,22 @@ contains
 
   end function size_const
 
-
   function s_vect_get_nrows(x) result(res)
     implicit none 
     class(psb_s_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = x%v%get_nrows()
   end function s_vect_get_nrows
 
+  function s_vect_sizeof(x) result(res)
+    implicit none 
+    class(psb_s_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = 0
+    if (allocated(x%v)) res = x%v%sizeof()
+  end function s_vect_sizeof
+
   function s_vect_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_s_vect_type), intent(inout) :: x, y

base/modules/psb_z_base_vect_mod.f90

@@ -7,6 +7,7 @@ module psb_z_base_vect_mod
     complex(psb_dpk_), allocatable :: v(:)
   contains
     procedure, pass(x) :: get_nrows => z_base_get_nrows
+    procedure, pass(x) :: sizeof   => z_base_sizeof
     procedure, pass(x) :: dot_v    => z_base_dot_v
     procedure, pass(x) :: dot_a    => z_base_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -146,15 +147,21 @@ contains
 
   end function size_const
     
-
   function z_base_get_nrows(x) result(res)
     implicit none 
     class(psb_z_base_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = size(x%v)
   end function z_base_get_nrows
 
+  function z_base_sizeof(x) result(res)
+    implicit none 
+    class(psb_z_base_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = (2*psb_sizeof_dp)*x%get_nrows()
+  end function z_base_sizeof
+
   function z_base_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_z_base_vect_type), intent(inout) :: x, y

base/modules/psb_z_vect_mod.f90

@@ -6,6 +6,7 @@ module psb_z_vect_mod
     class(psb_z_base_vect_type), allocatable :: v 
   contains
     procedure, pass(x) :: get_nrows => z_vect_get_nrows
+    procedure, pass(x) :: sizeof   => z_vect_sizeof
     procedure, pass(x) :: dot_v    => z_vect_dot_v
     procedure, pass(x) :: dot_a    => z_vect_dot_a
     generic, public    :: dot      => dot_v, dot_a
@@ -148,15 +149,22 @@ contains
 
   end function size_const
 
-
   function z_vect_get_nrows(x) result(res)
     implicit none 
     class(psb_z_vect_type), intent(in) :: x
     integer :: res
-    res = -1
+    res = 0
     if (allocated(x%v)) res = x%v%get_nrows()
   end function z_vect_get_nrows
 
+  function z_vect_sizeof(x) result(res)
+    implicit none 
+    class(psb_z_vect_type), intent(in) :: x
+    integer(psb_long_int_k_) :: res
+    res = 0
+    if (allocated(x%v)) res = x%v%sizeof()
+  end function z_vect_sizeof
+
   function z_vect_dot_v(n,x,y) result(res)
     implicit none 
     class(psb_z_vect_type), intent(inout) :: x, y

util/psb_c_renum_impl.F90

@@ -1,3 +1,31 @@
+subroutine psb_c_mat_renums(alg,mat,info,perm)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_c_mat_renums
+  implicit none 
+  character(len=*), intent(in) :: alg
+  type(psb_cspmat_type), intent(inout) :: mat
+  integer, intent(out) :: info
+  integer, allocatable, optional, intent(out) :: perm(:)
+  
+  integer            :: err_act, nr, nc, ialg
+  character(len=20)  :: name
+
+  info = psb_success_
+  name = 'mat_renum'
+  call psb_erractionsave(err_act)
+
+  info = psb_success_
+  select case (psb_toupper(alg))
+  case ('GPS')
+    ialg = psb_mat_renum_gps_
+  case ('AMD')
+    ialg = psb_mat_renum_amd_
+  case default
+    ialg = -1
+  end select
+  call psb_mat_renum(ialg,mat,info,perm)
+end subroutine psb_c_mat_renums
+  
 subroutine psb_c_mat_renum(alg,mat,info,perm)
   use psb_base_mod
   use psb_renum_mod, psb_protect_name => psb_c_mat_renum
@@ -7,8 +35,8 @@ subroutine psb_c_mat_renum(alg,mat,info,perm)
   integer, intent(out) :: info
   integer, allocatable, optional, intent(out) :: perm(:)
   
-  integer :: err_act
-  character(len=20)           :: name
+  integer            :: err_act, nr, nc
+  character(len=20)  :: name
 
   info = psb_success_
   name = 'mat_renum'
@@ -16,11 +44,23 @@ subroutine psb_c_mat_renum(alg,mat,info,perm)
 
   info = psb_success_
   
+  nr = mat%get_nrows()
+  nc = mat%get_ncols()
+  if (nr /= nc) then 
+    info = psb_err_rectangular_mat_unsupported_
+    call psb_errpush(info,name,i_err=(/nr,nc,0,0,0/))
+    goto 9999
+  end if
+
   select case (alg)
   case(psb_mat_renum_gps_) 
 
     call psb_mat_renum_gps(mat,info,perm)
 
+  case(psb_mat_renum_amd_) 
+
+    call psb_mat_renum_amd(mat,info,perm)
+
   case default
     info = psb_err_input_value_invalid_i_
     call psb_errpush(info,name,i_err=(/1,alg,0,0,0/))
@@ -138,5 +178,162 @@ contains
     return
   end subroutine psb_mat_renum_gps
 
+
+  subroutine psb_mat_renum_amd(a,info,operm)
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    use iso_c_binding
+#endif
+    use psb_base_mod    
+    implicit none 
+    type(psb_cspmat_type), intent(inout) :: a
+    integer, intent(out) :: info
+    integer, allocatable, optional, intent(out) :: operm(:)
+
+    ! 
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    interface 
+      function psb_amd_order(n,ap,ai,p)&
+           & result(res) bind(c,name='psb_amd_order')
+        use iso_c_binding
+        integer(c_int) :: res
+        integer(c_int), value :: n
+        integer(c_int) :: ap(*), ai(*), p(*)
+      end function psb_amd_order
+    end interface
+#endif
+
+    type(psb_c_csc_sparse_mat)  :: acsc
+    class(psb_c_base_sparse_mat), allocatable :: aa
+    type(psb_c_coo_sparse_mat)  :: acoo
+
+    integer :: err_act
+    character(len=20)           :: name
+    integer :: i, j, k, ideg, nr, ibw, ipf, idpth, nz
+
+    info = psb_success_
+    name = 'mat_renum_amd'
+    call psb_erractionsave(err_act)
+
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+
+    info = psb_success_
+    nr   = a%get_nrows()
+    nz   = a%get_nzeros()
+    allocate(perm(nr),stat=info)
+    if (info /= 0) then 
+      info = psb_err_alloc_dealloc_
+      call psb_errpush(info,name) 
+      goto 9999
+    end if
+
+
+    allocate(aa, mold=a%a)
+    call a%mv_to(acsc)
+
+    acsc%ia(:)  = acsc%ia(:) - 1
+    acsc%icp(:) = acsc%icp(:) - 1
+
+    info = psb_amd_order(nr,acsc%icp,acsc%ia,perm)
+    if (info /= psb_success_) then 
+      info = psb_err_from_subroutine_
+      call psb_errpush(info,name,a_err='psb_amd_order') 
+      goto 9999
+    end if
+
+    perm(:)     = perm(:) + 1
+    acsc%ia(:)  = acsc%ia(:) + 1
+    acsc%icp(:) = acsc%icp(:) + 1
+
+    call acsc%mv_to_coo(acoo,info)
+    do i=1, acoo%get_nzeros()
+      acoo%ia(i) = perm(acoo%ia(i))
+      acoo%ja(i) = perm(acoo%ja(i))
+    end do
+    call acoo%fix(info) 
+
+    ! Get back to where we started from
+    call aa%mv_from_coo(acoo,info)
+    call a%mv_from(aa)
+    if (present(operm)) then 
+      call psb_realloc(nr,operm,info)
+      if (info /= psb_success_) then 
+        info = psb_err_alloc_dealloc_
+        call psb_errpush(info,name) 
+        goto 9999
+      end if
+      operm(1:nr) = perm(1:nr)
+    end if
+
+    deallocate(aa,perm)
+
+#else 
+
+    info = psb_err_missing_aux_lib_
+    call psb_errpush(info, name) 
+    goto 9999
+#endif
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine psb_mat_renum_amd
+
 end subroutine psb_c_mat_renum
 
+
+subroutine psb_c_cmp_bwpf(mat,bwl,bwu,prf,info)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_c_cmp_bwpf
+  implicit none 
+  type(psb_cspmat_type), intent(in) :: mat
+  integer, intent(out) :: bwl, bwu
+  integer, intent(out) :: prf
+  integer, intent(out) :: info
+  !
+  integer, allocatable :: irow(:), icol(:)
+  complex(psb_spk_), allocatable :: val(:)
+  integer :: nz
+  integer :: i, j, lrbu, lrbl
+  
+  info = psb_success_
+  bwl = 0
+  bwu = 0
+  prf = 0
+  select type (aa=>mat%a)
+  class is (psb_c_csr_sparse_mat)
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      do j = aa%irp(i), aa%irp(i+1) - 1
+        lrbl = max(lrbl,i-aa%ja(j))
+        lrbu = max(lrbu,aa%ja(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+    
+  class default
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      call aa%csget(i,i,nz,irow,icol,val,info)
+      if (info /= psb_success_) return
+      do j=1, nz
+        lrbl = max(lrbl,i-icol(j))
+        lrbu = max(lrbu,icol(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+  end select
+  
+end subroutine psb_c_cmp_bwpf

util/psb_renum_mod.f90

@@ -20,6 +20,13 @@ module psb_renum_mod
       integer, intent(out) :: info
       integer, allocatable, optional, intent(out) :: perm(:)
     end subroutine psb_d_mat_renum
+    subroutine psb_s_mat_renums(alg,mat,info,perm)
+      import  psb_sspmat_type
+      character(len=*), intent(in) :: alg
+      type(psb_sspmat_type), intent(inout) :: mat
+      integer, intent(out) :: info
+      integer, allocatable, optional, intent(out) :: perm(:)
+    end subroutine psb_s_mat_renums
     subroutine psb_s_mat_renum(alg,mat,info,perm)
       import  psb_sspmat_type
       integer, intent(in) :: alg
@@ -27,6 +34,13 @@ module psb_renum_mod
       integer, intent(out) :: info
       integer, allocatable, optional, intent(out) :: perm(:)
     end subroutine psb_s_mat_renum
+    subroutine psb_z_mat_renums(alg,mat,info,perm)
+      import  psb_zspmat_type
+      character(len=*), intent(in) :: alg
+      type(psb_zspmat_type), intent(inout) :: mat
+      integer, intent(out) :: info
+      integer, allocatable, optional, intent(out) :: perm(:)
+    end subroutine psb_z_mat_renums
     subroutine psb_z_mat_renum(alg,mat,info,perm)
       import  psb_zspmat_type
       integer, intent(in) :: alg
@@ -34,6 +48,13 @@ module psb_renum_mod
       integer, intent(out) :: info
       integer, allocatable, optional, intent(out) :: perm(:)
     end subroutine psb_z_mat_renum
+    subroutine psb_c_mat_renums(alg,mat,info,perm)
+      import  psb_cspmat_type
+      character(len=*), intent(in) :: alg
+      type(psb_cspmat_type), intent(inout) :: mat
+      integer, intent(out) :: info
+      integer, allocatable, optional, intent(out) :: perm(:)
+    end subroutine psb_c_mat_renums
     subroutine psb_c_mat_renum(alg,mat,info,perm)
       import  psb_cspmat_type
       integer, intent(in) :: alg
@@ -45,6 +66,13 @@ module psb_renum_mod
 
   
   interface psb_cmp_bwpf
+    subroutine psb_s_cmp_bwpf(mat,bwl,bwu,prf,info)
+      import  psb_sspmat_type
+      type(psb_sspmat_type), intent(in) :: mat
+      integer, intent(out) :: bwl, bwu
+      integer, intent(out) :: prf
+      integer, intent(out) :: info
+    end subroutine psb_s_cmp_bwpf
     subroutine psb_d_cmp_bwpf(mat,bwl,bwu,prf,info)
       import  psb_dspmat_type
       type(psb_dspmat_type), intent(in) :: mat
@@ -52,6 +80,20 @@ module psb_renum_mod
       integer, intent(out) :: prf
       integer, intent(out) :: info
     end subroutine psb_d_cmp_bwpf
+    subroutine psb_c_cmp_bwpf(mat,bwl,bwu,prf,info)
+      import  psb_cspmat_type
+      type(psb_cspmat_type), intent(in) :: mat
+      integer, intent(out) :: bwl, bwu
+      integer, intent(out) :: prf
+      integer, intent(out) :: info
+    end subroutine psb_c_cmp_bwpf
+    subroutine psb_z_cmp_bwpf(mat,bwl,bwu,prf,info)
+      import  psb_zspmat_type
+      type(psb_zspmat_type), intent(in) :: mat
+      integer, intent(out) :: bwl, bwu
+      integer, intent(out) :: prf
+      integer, intent(out) :: info
+    end subroutine psb_z_cmp_bwpf
   end interface psb_cmp_bwpf
 
 

util/psb_s_renum_impl.F90

@@ -1,3 +1,31 @@
+subroutine psb_s_mat_renums(alg,mat,info,perm)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_s_mat_renums
+  implicit none 
+  character(len=*), intent(in) :: alg
+  type(psb_sspmat_type), intent(inout) :: mat
+  integer, intent(out) :: info
+  integer, allocatable, optional, intent(out) :: perm(:)
+  
+  integer            :: err_act, nr, nc, ialg
+  character(len=20)  :: name
+
+  info = psb_success_
+  name = 'mat_renum'
+  call psb_erractionsave(err_act)
+
+  info = psb_success_
+  select case (psb_toupper(alg))
+  case ('GPS')
+    ialg = psb_mat_renum_gps_
+  case ('AMD')
+    ialg = psb_mat_renum_amd_
+  case default
+    ialg = -1
+  end select
+  call psb_mat_renum(ialg,mat,info,perm)
+end subroutine psb_s_mat_renums
+  
 subroutine psb_s_mat_renum(alg,mat,info,perm)
   use psb_base_mod
   use psb_renum_mod, psb_protect_name => psb_s_mat_renum
@@ -7,8 +35,8 @@ subroutine psb_s_mat_renum(alg,mat,info,perm)
   integer, intent(out) :: info
   integer, allocatable, optional, intent(out) :: perm(:)
   
-  integer :: err_act
-  character(len=20)           :: name
+  integer            :: err_act, nr, nc
+  character(len=20)  :: name
 
   info = psb_success_
   name = 'mat_renum'
@@ -16,11 +44,23 @@ subroutine psb_s_mat_renum(alg,mat,info,perm)
 
   info = psb_success_
   
+  nr = mat%get_nrows()
+  nc = mat%get_ncols()
+  if (nr /= nc) then 
+    info = psb_err_rectangular_mat_unsupported_
+    call psb_errpush(info,name,i_err=(/nr,nc,0,0,0/))
+    goto 9999
+  end if
+
   select case (alg)
   case(psb_mat_renum_gps_) 
 
     call psb_mat_renum_gps(mat,info,perm)
 
+  case(psb_mat_renum_amd_) 
+
+    call psb_mat_renum_amd(mat,info,perm)
+
   case default
     info = psb_err_input_value_invalid_i_
     call psb_errpush(info,name,i_err=(/1,alg,0,0,0/))
@@ -138,5 +178,162 @@ contains
     return
   end subroutine psb_mat_renum_gps
 
+
+  subroutine psb_mat_renum_amd(a,info,operm)
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    use iso_c_binding
+#endif
+    use psb_base_mod    
+    implicit none 
+    type(psb_sspmat_type), intent(inout) :: a
+    integer, intent(out) :: info
+    integer, allocatable, optional, intent(out) :: operm(:)
+
+    ! 
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    interface 
+      function psb_amd_order(n,ap,ai,p)&
+           & result(res) bind(c,name='psb_amd_order')
+        use iso_c_binding
+        integer(c_int) :: res
+        integer(c_int), value :: n
+        integer(c_int) :: ap(*), ai(*), p(*)
+      end function psb_amd_order
+    end interface
+#endif
+
+    type(psb_s_csc_sparse_mat)  :: acsc
+    class(psb_s_base_sparse_mat), allocatable :: aa
+    type(psb_s_coo_sparse_mat)  :: acoo
+
+    integer :: err_act
+    character(len=20)           :: name
+    integer :: i, j, k, ideg, nr, ibw, ipf, idpth, nz
+
+    info = psb_success_
+    name = 'mat_renum_amd'
+    call psb_erractionsave(err_act)
+
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+
+    info = psb_success_
+    nr   = a%get_nrows()
+    nz   = a%get_nzeros()
+    allocate(perm(nr),stat=info)
+    if (info /= 0) then 
+      info = psb_err_alloc_dealloc_
+      call psb_errpush(info,name) 
+      goto 9999
+    end if
+
+
+    allocate(aa, mold=a%a)
+    call a%mv_to(acsc)
+
+    acsc%ia(:)  = acsc%ia(:) - 1
+    acsc%icp(:) = acsc%icp(:) - 1
+
+    info = psb_amd_order(nr,acsc%icp,acsc%ia,perm)
+    if (info /= psb_success_) then 
+      info = psb_err_from_subroutine_
+      call psb_errpush(info,name,a_err='psb_amd_order') 
+      goto 9999
+    end if
+
+    perm(:)     = perm(:) + 1
+    acsc%ia(:)  = acsc%ia(:) + 1
+    acsc%icp(:) = acsc%icp(:) + 1
+
+    call acsc%mv_to_coo(acoo,info)
+    do i=1, acoo%get_nzeros()
+      acoo%ia(i) = perm(acoo%ia(i))
+      acoo%ja(i) = perm(acoo%ja(i))
+    end do
+    call acoo%fix(info) 
+
+    ! Get back to where we started from
+    call aa%mv_from_coo(acoo,info)
+    call a%mv_from(aa)
+    if (present(operm)) then 
+      call psb_realloc(nr,operm,info)
+      if (info /= psb_success_) then 
+        info = psb_err_alloc_dealloc_
+        call psb_errpush(info,name) 
+        goto 9999
+      end if
+      operm(1:nr) = perm(1:nr)
+    end if
+
+    deallocate(aa,perm)
+
+#else 
+
+    info = psb_err_missing_aux_lib_
+    call psb_errpush(info, name) 
+    goto 9999
+#endif
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine psb_mat_renum_amd
+
 end subroutine psb_s_mat_renum
 
+
+subroutine psb_s_cmp_bwpf(mat,bwl,bwu,prf,info)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_s_cmp_bwpf
+  implicit none 
+  type(psb_sspmat_type), intent(in) :: mat
+  integer, intent(out) :: bwl, bwu
+  integer, intent(out) :: prf
+  integer, intent(out) :: info
+  !
+  integer, allocatable :: irow(:), icol(:)
+  real(psb_spk_), allocatable :: val(:)
+  integer :: nz
+  integer :: i, j, lrbu, lrbl
+  
+  info = psb_success_
+  bwl = 0
+  bwu = 0
+  prf = 0
+  select type (aa=>mat%a)
+  class is (psb_s_csr_sparse_mat)
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      do j = aa%irp(i), aa%irp(i+1) - 1
+        lrbl = max(lrbl,i-aa%ja(j))
+        lrbu = max(lrbu,aa%ja(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+    
+  class default
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      call aa%csget(i,i,nz,irow,icol,val,info)
+      if (info /= psb_success_) return
+      do j=1, nz
+        lrbl = max(lrbl,i-icol(j))
+        lrbu = max(lrbu,icol(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+  end select
+  
+end subroutine psb_s_cmp_bwpf

util/psb_z_renum_impl.F90

@@ -1,3 +1,31 @@
+subroutine psb_z_mat_renums(alg,mat,info,perm)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_z_mat_renums
+  implicit none 
+  character(len=*), intent(in) :: alg
+  type(psb_zspmat_type), intent(inout) :: mat
+  integer, intent(out) :: info
+  integer, allocatable, optional, intent(out) :: perm(:)
+  
+  integer            :: err_act, nr, nc, ialg
+  character(len=20)  :: name
+
+  info = psb_success_
+  name = 'mat_renum'
+  call psb_erractionsave(err_act)
+
+  info = psb_success_
+  select case (psb_toupper(alg))
+  case ('GPS')
+    ialg = psb_mat_renum_gps_
+  case ('AMD')
+    ialg = psb_mat_renum_amd_
+  case default
+    ialg = -1
+  end select
+  call psb_mat_renum(ialg,mat,info,perm)
+end subroutine psb_z_mat_renums
+  
 subroutine psb_z_mat_renum(alg,mat,info,perm)
   use psb_base_mod
   use psb_renum_mod, psb_protect_name => psb_z_mat_renum
@@ -7,8 +35,8 @@ subroutine psb_z_mat_renum(alg,mat,info,perm)
   integer, intent(out) :: info
   integer, allocatable, optional, intent(out) :: perm(:)
   
-  integer :: err_act
-  character(len=20)           :: name
+  integer            :: err_act, nr, nc
+  character(len=20)  :: name
 
   info = psb_success_
   name = 'mat_renum'
@@ -16,11 +44,23 @@ subroutine psb_z_mat_renum(alg,mat,info,perm)
 
   info = psb_success_
   
+  nr = mat%get_nrows()
+  nc = mat%get_ncols()
+  if (nr /= nc) then 
+    info = psb_err_rectangular_mat_unsupported_
+    call psb_errpush(info,name,i_err=(/nr,nc,0,0,0/))
+    goto 9999
+  end if
+
   select case (alg)
   case(psb_mat_renum_gps_) 
 
     call psb_mat_renum_gps(mat,info,perm)
 
+  case(psb_mat_renum_amd_) 
+
+    call psb_mat_renum_amd(mat,info,perm)
+
   case default
     info = psb_err_input_value_invalid_i_
     call psb_errpush(info,name,i_err=(/1,alg,0,0,0/))
@@ -138,4 +178,162 @@ contains
     return
   end subroutine psb_mat_renum_gps
 
+
+  subroutine psb_mat_renum_amd(a,info,operm)
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    use iso_c_binding
+#endif
+    use psb_base_mod    
+    implicit none 
+    type(psb_zspmat_type), intent(inout) :: a
+    integer, intent(out) :: info
+    integer, allocatable, optional, intent(out) :: operm(:)
+
+    ! 
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+    interface 
+      function psb_amd_order(n,ap,ai,p)&
+           & result(res) bind(c,name='psb_amd_order')
+        use iso_c_binding
+        integer(c_int) :: res
+        integer(c_int), value :: n
+        integer(c_int) :: ap(*), ai(*), p(*)
+      end function psb_amd_order
+    end interface
+#endif
+
+    type(psb_z_csc_sparse_mat)  :: acsc
+    class(psb_z_base_sparse_mat), allocatable :: aa
+    type(psb_z_coo_sparse_mat)  :: acoo
+
+    integer :: err_act
+    character(len=20)           :: name
+    integer :: i, j, k, ideg, nr, ibw, ipf, idpth, nz
+
+    info = psb_success_
+    name = 'mat_renum_amd'
+    call psb_erractionsave(err_act)
+
+#if defined(HAVE_AMD) && defined(HAVE_ISO_C_BINDING)
+
+    info = psb_success_
+    nr   = a%get_nrows()
+    nz   = a%get_nzeros()
+    allocate(perm(nr),stat=info)
+    if (info /= 0) then 
+      info = psb_err_alloc_dealloc_
+      call psb_errpush(info,name) 
+      goto 9999
+    end if
+
+
+    allocate(aa, mold=a%a)
+    call a%mv_to(acsc)
+
+    acsc%ia(:)  = acsc%ia(:) - 1
+    acsc%icp(:) = acsc%icp(:) - 1
+
+    info = psb_amd_order(nr,acsc%icp,acsc%ia,perm)
+    if (info /= psb_success_) then 
+      info = psb_err_from_subroutine_
+      call psb_errpush(info,name,a_err='psb_amd_order') 
+      goto 9999
+    end if
+
+    perm(:)     = perm(:) + 1
+    acsc%ia(:)  = acsc%ia(:) + 1
+    acsc%icp(:) = acsc%icp(:) + 1
+
+    call acsc%mv_to_coo(acoo,info)
+    do i=1, acoo%get_nzeros()
+      acoo%ia(i) = perm(acoo%ia(i))
+      acoo%ja(i) = perm(acoo%ja(i))
+    end do
+    call acoo%fix(info) 
+
+    ! Get back to where we started from
+    call aa%mv_from_coo(acoo,info)
+    call a%mv_from(aa)
+    if (present(operm)) then 
+      call psb_realloc(nr,operm,info)
+      if (info /= psb_success_) then 
+        info = psb_err_alloc_dealloc_
+        call psb_errpush(info,name) 
+        goto 9999
+      end if
+      operm(1:nr) = perm(1:nr)
+    end if
+
+    deallocate(aa,perm)
+
+#else 
+
+    info = psb_err_missing_aux_lib_
+    call psb_errpush(info, name) 
+    goto 9999
+#endif
+
+    call psb_erractionrestore(err_act)
+    return
+
+9999 continue
+    call psb_erractionrestore(err_act)
+    if (err_act == psb_act_abort_) then
+      call psb_error()
+      return
+    end if
+    return
+  end subroutine psb_mat_renum_amd
+
 end subroutine psb_z_mat_renum
+
+
+subroutine psb_z_cmp_bwpf(mat,bwl,bwu,prf,info)
+  use psb_base_mod
+  use psb_renum_mod, psb_protect_name => psb_z_cmp_bwpf
+  implicit none 
+  type(psb_zspmat_type), intent(in) :: mat
+  integer, intent(out) :: bwl, bwu
+  integer, intent(out) :: prf
+  integer, intent(out) :: info
+  !
+  integer, allocatable :: irow(:), icol(:)
+  complex(psb_dpk_), allocatable :: val(:)
+  integer :: nz
+  integer :: i, j, lrbu, lrbl
+  
+  info = psb_success_
+  bwl = 0
+  bwu = 0
+  prf = 0
+  select type (aa=>mat%a)
+  class is (psb_z_csr_sparse_mat)
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      do j = aa%irp(i), aa%irp(i+1) - 1
+        lrbl = max(lrbl,i-aa%ja(j))
+        lrbu = max(lrbu,aa%ja(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+    
+  class default
+    do i=1, aa%get_nrows()
+      lrbl = 0
+      lrbu = 0
+      call aa%csget(i,i,nz,irow,icol,val,info)
+      if (info /= psb_success_) return
+      do j=1, nz
+        lrbl = max(lrbl,i-icol(j))
+        lrbu = max(lrbu,icol(j)-i)
+      end do
+      prf = prf + lrbl+lrbu
+      bwu  = max(bwu,lrbu)
+      bwl  = max(bwl,lrbu)
+    end do
+  end select
+  
+end subroutine psb_z_cmp_bwpf