diff --git a/src/match_order.f90 b/src/match_order.f90
index 4be07f2d..6311f4a1 100644
--- a/src/match_order.f90
+++ b/src/match_order.f90
@@ -568,7 +568,7 @@ subroutine mo_match(n,row2,ptr2,val2,scale,flag,stat,perm)
 
     k = 0
     do i = 1, n
-       if (cperm(i) .lt. 0) then
+       if (cperm(i) .eq. 0) then
           ! row i and col j are not part of the matching
           old_to_new(i) = -1
        else
@@ -590,11 +590,11 @@ subroutine mo_match(n,row2,ptr2,val2,scale,flag,stat,perm)
        j1 = j2
        j2 = ptr2(i+1)
        ! skip over unmatched entries
-       if (cperm(i) .lt. 0) cycle
+       if (cperm(i) .eq. 0) cycle
        k = k + 1
        do jlong = j1, j2-1
           jj = row2(jlong)
-          if (cperm(jj) .lt. 0) cycle
+          if (cperm(jj) .eq. 0) cycle
           nne = nne + 1
           row2(nne) = old_to_new(jj)
           val2(nne) = val2(jlong)
@@ -611,7 +611,7 @@ subroutine mo_match(n,row2,ptr2,val2,scale,flag,stat,perm)
 
     do i = 1, n
        j = old_to_new(i)
-       if (j .lt. 0) then
+       if (j .eq. 0) then
           scale(i) = -huge(scale)
        else
           ! Note: we need to subtract col max using old matrix numbering
diff --git a/src/scaling.f90 b/src/scaling.f90
index 90f692d1..2a8e9fcb 100644
--- a/src/scaling.f90
+++ b/src/scaling.f90
@@ -1168,29 +1168,6 @@ subroutine hungarian_match(m,n,ptr,row,val,iperm,num,dualu,dualv,st)
    ! Zero dual row variables for unmatched rows
    where (iperm(1:m) .eq. 0) dualu(1:m) = 0.0
 
-   ! Return if matrix has full structural rank
-   if (num .eq. min(m,n)) return
-
-   ! Otherwise, matrix is structurally singular, complete iperm.
-   ! jperm, out are work arrays
-   jperm(1:n) = 0
-   k = 0
-   do i = 1, m
-      if (iperm(i) .eq. 0) then
-         k = k + 1
-         out(k) = i
-      else
-         j = iperm(i)
-         jperm(j) = i
-      end if
-   end do
-   k = 0
-   do j = 1, n
-      if (jperm(j) .ne. 0) cycle
-      k = k + 1
-      jdum = int(out(k))
-      iperm(jdum) = -j
-   end do
  end subroutine hungarian_match
 
 !**********************************************************************
diff --git a/tests/scaling.f90 b/tests/scaling.f90
index f5e38833..40c7aef6 100644
--- a/tests/scaling.f90
+++ b/tests/scaling.f90
@@ -32,6 +32,8 @@ program main
    call test_equilib_sym_random
    call test_equilib_unsym_random
    call test_hungarian_sym_random
+   call test_hungarian_unsym_singular
+   call test_hungarian_sym_singular
    call test_hungarian_unsym_random
 
    write(*, "(/a)") "=========================="
@@ -597,6 +599,104 @@ end subroutine test_equilib_unsym_random
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+subroutine test_hungarian_unsym_singular
+  integer :: m = 3
+  integer :: n = 5
+  integer :: nz = 6
+  integer :: ising = 3
+  type(matrix_type) :: a
+
+  type(hungarian_options) :: options
+  type(hungarian_inform) :: inform
+
+  integer, allocatable, dimension(:) :: match
+  real(wp), allocatable, dimension(:) :: rscaling, cscaling
+
+  write(*, "(a)")
+  write(*, "(a)") "===================================================="
+  write(*, "(a)") "Testing hungarian_scale_unsym() with singular matrix"
+  write(*, "(a)") "===================================================="
+
+  allocate(a%ptr(n+1))
+  allocate(a%row(nz), a%val(nz))
+  allocate(rscaling(m), cscaling(n), match(m))
+
+  ! Produce warning rather than error
+  options%scale_if_singular = .true.
+
+  a%n = n
+  a%m = m
+
+  a%ptr(1:n+1)          = (/   1,   3,   5,   6,   6,   7 /)
+  a%row(1:a%ptr(n+1)-1) = (/   1,   2,   1,   2,   2,   2 /)
+  a%val(1:a%ptr(n+1)-1) = (/ 2.0, 1.0, 1.0, 4.0, 1.0, 1.0 /)
+
+  call hungarian_scale_unsym(a%m, a%n, a%ptr, a%row, a%val, rscaling, &
+         cscaling, options, inform, match=match)
+
+  if(inform%flag .ne. 1) then
+     write(*, "(a, i5)") "Returned inform%flag = ", inform%flag
+     errors = errors + 1
+  endif
+
+  if(match(ising) .ne. 0) then
+     write(*, "(a, i5, a, i5)") "Singular row ", ising, "matched to ", match(ising)
+     errors = errors + 1
+  endif
+
+end  subroutine test_hungarian_unsym_singular
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine test_hungarian_sym_singular
+  integer :: m = 3
+  integer :: n = 3
+  integer :: nz = 2
+  integer :: ising = 3
+  type(matrix_type) :: a
+
+  type(hungarian_options) :: options
+  type(hungarian_inform) :: inform
+
+  integer, allocatable, dimension(:) :: match
+  real(wp), allocatable, dimension(:) :: scaling
+
+  write(*, "(a)")
+  write(*, "(a)") "===================================================="
+  write(*, "(a)") "Testing hungarian_scale_sym() with singular matrix"
+  write(*, "(a)") "===================================================="
+
+  allocate(a%ptr(n+1))
+  allocate(a%row(nz), a%val(nz))
+  allocate(scaling(n), match(m))
+
+  ! Produce warning rather than error
+  options%scale_if_singular = .true.
+
+  a%n = n
+  a%m = m
+
+  a%ptr(1:n+1)          = (/   1,   2,   3,  3 /)
+  a%row(1:a%ptr(n+1)-1) = (/   1,   2/)
+  a%val(1:a%ptr(n+1)-1) = (/ 2.0, 1.0/)
+
+  call hungarian_scale_sym(a%n, a%ptr, a%row, a%val, scaling, &
+         options, inform, match=match)
+
+  if(inform%flag .ne. 1) then
+     write(*, "(a, i5)") "Returned inform%flag = ", inform%flag
+     errors = errors + 1
+  endif
+
+  if(match(ising) .ne. 0) then
+     write(*, "(a, i5, a, i5)") "Singular column ", ising, " has value ", match(ising)
+     errors = errors + 1
+  endif
+
+end  subroutine test_hungarian_sym_singular
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 subroutine test_hungarian_sym_random
    integer :: maxn = 1000
    integer :: maxnz =  1000000