From f481a108b1210b73041f3c82754e434f1307ecdf Mon Sep 17 00:00:00 2001
From: Christoph Hansknecht <christoph.hansknecht@tu-clausthal.de>
Date: Fri, 14 Jun 2024 23:41:09 +0200
Subject: [PATCH] Fix memory fault during scaling of singular matrix

---
 src/scaling.f90   |  2 +-
 tests/scaling.f90 | 50 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 51 insertions(+), 1 deletion(-)

diff --git a/src/scaling.f90 b/src/scaling.f90
index 90f692d1..fedcde57 100644
--- a/src/scaling.f90
+++ b/src/scaling.f90
@@ -1189,7 +1189,7 @@ subroutine hungarian_match(m,n,ptr,row,val,iperm,num,dualu,dualv,st)
       if (jperm(j) .ne. 0) cycle
       k = k + 1
       jdum = int(out(k))
-      iperm(jdum) = -j
+      iperm(jdum) = 0
    end do
  end subroutine hungarian_match
 
diff --git a/tests/scaling.f90 b/tests/scaling.f90
index 704fefdd..53a10484 100644
--- a/tests/scaling.f90
+++ b/tests/scaling.f90
@@ -32,6 +32,7 @@ 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_unsym_random
 
    write(*, "(/a)") "=========================="
@@ -573,6 +574,55 @@ 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 equilib_scaling_unsym() with singular matrix"
+  write(*, "(a)") "===================================================="
+
+  allocate(a%ptr(n+1))
+  allocate(a%row(nz), a%val(nz))
+  allocate(rscaling(m), cscaling(n), match(n))
+
+  ! Produce warning rather than error
+  options%scale_if_singular = .true.
+
+  a%n = n
+  a%m = m
+
+  a%ptr(1:n+1)          = (/   1,   3,   5,   5,   5,   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 column ", ising, "matched to ", match(ising)
+     errors = errors + 1
+  endif
+
+end  subroutine test_hungarian_unsym_singular
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
 subroutine test_hungarian_sym_random
    integer :: maxn = 1000
    integer :: maxnz =  1000000