From 8ff8f7c1bba767e1de10099ca0c7249a60603edb Mon Sep 17 00:00:00 2001
From: Alexis Montoison <alexis.montoison@polymtl.ca>
Date: Thu, 5 Oct 2023 20:32:28 -0500
Subject: [PATCH] Patch libHSL upstream for ma57_get_factors

---
 src/HSL.jl                  |   1 -
 src/hsl_ma57.jl             | 178 ++++++++++++++++++++++++++++++++++++
 src/hsl_ma57_patch.jl       | 177 -----------------------------------
 test/runtests.jl            |   1 -
 test/test_hsl_ma57.jl       |  68 ++++++++++++++
 test/test_hsl_ma57_patch.jl |  69 --------------
 6 files changed, 246 insertions(+), 248 deletions(-)
 delete mode 100644 src/hsl_ma57_patch.jl
 delete mode 100644 test/test_hsl_ma57_patch.jl

diff --git a/src/HSL.jl b/src/HSL.jl
index 74294b5..bd28bad 100644
--- a/src/HSL.jl
+++ b/src/HSL.jl
@@ -35,7 +35,6 @@ include("wrappers.jl")
 
 # Interfaces
 include("hsl_ma57.jl")
-# include("hsl_ma57_patch.jl")
 include("hsl_ma97.jl")
 include("kb07.jl")
 include("mc21.jl")
diff --git a/src/hsl_ma57.jl b/src/hsl_ma57.jl
index fe189b5..de478f2 100644
--- a/src/hsl_ma57.jl
+++ b/src/hsl_ma57.jl
@@ -676,3 +676,181 @@ end
 
 ma57_least_squares(A::Matrix{T}, b::Array{T}) where {T <: Ma57Data} =
   ma57_least_squares(sparse(A), b)
+
+## get factors -----------------------------------------------------------------
+for (fname, T) in ((:__hsl_ma57_single_MOD_ma57lf , :Float32),
+                   (:__hsl_ma57_double_MOD_ma57lfd, :Float64))
+  @eval begin
+    function ma57_get_factors(ma57::Ma57{$T})
+
+      # make room for L factor
+      nebdu = ma57.info.info[14]
+      nzl = nebdu
+      ipl = Vector{Cint}(undef, ma57.n + 1)
+      irn = Vector{Cint}(undef, nzl)
+      fl = Vector{$T}(undef, nzl)
+
+      # make room for D; note that entire 2x2 blocks are stored
+      nzd = 2 * ma57.info.num_2x2_pivots + ma57.n
+      ipd = Vector{Cint}(undef, ma57.n + 1)
+      id = Vector{Cint}(undef, nzd)
+      d = Vector{$T}(undef, nzd)
+
+      ivp = Vector{Cint}(undef, ma57.n)
+      iperm = Vector{Cint}(undef, ma57.n)
+
+      status = 0
+
+      @ccall libhsl.$fname(ma57.n::Ref{Cint}, ma57.__fact::Ptr{$T}, ma57.__lfact::Ref{Cint}, ma57.__ifact::Ptr{Cint},
+                           ma57.__lifact::Ref{Cint}, nebdu::Ref{Cint}, nzl::Ref{Cint}, ipl::Ptr{Cint}, irn::Ptr{Cint},
+                           fl::Ptr{$T}, nzd::Ref{Cint}, ipd::Ptr{Cint}, id::Ptr{Cint}, d::Ptr{$T}, ivp::Ptr{Cint},
+                           iperm::Ptr{Cint}, ma57.info.rank::Ref{Cint}, status::Ref{Cint})::Cvoid
+
+      status < 0 && throw(Ma57Exception("Ma57: Error while retrieving factors", status))
+      s = ma57.control.icntl[15] == 1 ? ma57.__fact[(end - ma57.n):(end - 1)] : ones($T, ma57.n)
+      L = SparseMatrixCSC(ma57.n, ma57.n, ipl, irn, fl)
+      D = SparseMatrixCSC(ma57.n, ma57.n, ipd, id, d)
+
+      return (L, D, s, iperm)
+    end
+  end
+end
+
+## alter D ---------------------------------------------------------------------
+function ma57_alter_d(ma57::Ma57{T}, d::Matrix{T}) where {T <: Ma57Data}
+  ka = 1
+  k2 = ma57.__ifact[1]
+  kd = 0
+  kw = 4
+  for blk = 1:abs(ma57.__ifact[3])
+    ncols = ma57.__ifact[kw]
+    nrows = ma57.__ifact[kw + 1]
+    two = false
+    for i = 1:nrows
+      kd = kd + 1
+      ma57.__fact[ka] = d[1, kd]
+      if ma57.__ifact[kw + 1 + i] < 0
+        two = !two
+      end
+      if two
+        ma57.__fact[k2] = d[2, kd]
+        k2 = k2 + 1
+      else
+        status = ma57.info.info[1]
+        d[2, kd] == 0.0 ||
+          throw(Ma57Exception("Ma57: Erroneous 2x2 block specified in d[2,$kd]", status))
+      end
+      ka = ka + nrows + 1 - i
+    end
+    ka = ka + nrows * (ncols - nrows)
+    kw = kw + ncols + 2
+  end
+end
+
+## additional docstrings -------------------------------------------------------
+"""
+# Retrieve factors, scaling, and permutation.
+
+  L, D, s, iperm = ma57_get_factors(M)
+
+  Function will retrieve a unit triangular matrix L, a block-diagonal matrix D a scaling
+  vector s and a permutation vector p such that
+
+  P * S * A * S * P' = L * D * L'
+
+  where S = spdiagm(s) and P = speye(n)[p,:].
+
+  Note that the numerical factorization must have been performed and have succeeded.
+
+## Input arguments
+
+* `M::Ma57`: factorized `Ma57` object
+
+## Return values
+
+* `L::SparseMatrixCSC{T<:Ma57Data,Int}`: L factor
+* `D::SparseMatrixCSC{T<:Ma57Data,Int}`: D factor
+* `s::Vector{T}`: diagonal components of scaling matrix S
+* `iperm::Vector{Int}`: array representin permutation matrix P
+
+## Example:
+
+```julia
+
+julia> using HSL
+
+julia> T = Float64;
+
+julia> rows = Cint[1, 1, 2, 2, 3, 3, 5]; cols = Cint[1, 2, 3, 5, 3, 4, 5];
+
+julia> vals = T[2, 3, 4, 6, 1, 5, 1];
+
+julia> A = sparse(rows, cols, vals); A = A + triu(A, 1)';
+
+julia> b = T[8, 45, 31, 15, 17]
+
+julia> ϵ = sqrt(eps(eltype(A)))
+
+julia> xexact = T[1, 2, 3, 4, 5]
+
+julia> M = Ma57(A)
+
+julia> ma57_factorize(M)
+
+julia> (L, D, s, p) = ma57_get_factors(M)
+
+julia> S = spdiagm(s)
+
+julia> P = speye(M.n)[p, :]
+
+julia> vecnorm(P * S * A * S * P' - L * D * L') ≤ ϵ * vecnorm(A)
+
+true
+```
+"""
+ma57_get_factors
+
+"""
+# Alter block diagonal matrix `D`
+
+## Input arguments:
+
+* `M::Ma57`: `Ma57` object
+* `F::Matrix{Float64}`: diagonal adjustment
+
+
+## Example:
+
+```julia
+julia> using HSL
+
+julia> T = Float64;
+
+julia> rows = Cint[1, 1, 2, 2, 3, 3, 5]; cols = Cint[1, 2, 3, 5, 3, 4, 5];
+
+julia> vals = T[2, 3, 4, 6, 1, 5, 1];
+
+julia> A = sparse(rows, cols, vals); A = A + triu(A, 1)';
+
+julia> b = T[8, 45, 31, 15, 17]
+
+julia> ϵ = sqrt(eps(eltype(A)))
+
+julia> xexact = T[1, 2, 3, 4, 5]
+
+julia> M = Ma57(A)
+
+julia> ma57_factorize(M)
+
+julia> (L, D, s, p) = ma57_get_factors(M)
+
+julia> d1 = abs.(diag(D))
+
+julia> d2 = [diag(D, 1) ; 0]
+
+julia> F = [full(d1)' ; full(d2)']
+
+julia> ma57_alter_d(M, F)
+```
+"""
+ma57_alter_d
diff --git a/src/hsl_ma57_patch.jl b/src/hsl_ma57_patch.jl
deleted file mode 100644
index cd9f308..0000000
--- a/src/hsl_ma57_patch.jl
+++ /dev/null
@@ -1,177 +0,0 @@
-## get factors -----------------------------------------------------------------
-for (fname, T) in ((:ma57lf_ , :Float32),
-                   (:ma57lfd_, :Float64))
-  @eval begin
-    function ma57_get_factors(ma57::Ma57{$T})
-
-      # make room for L factor
-      nebdu = ma57.info.info[14]
-      nzl = nebdu
-      ipl = Vector{Cint}(undef, ma57.n + 1)
-      irn = Vector{Cint}(undef, nzl)
-      fl = Vector{$T}(undef, nzl)
-
-      # make room for D; note that entire 2x2 blocks are stored
-      nzd = 2 * ma57.info.num_2x2_pivots + ma57.n
-      ipd = Vector{Cint}(undef, ma57.n + 1)
-      id = Vector{Cint}(undef, nzd)
-      d = Vector{$T}(undef, nzd)
-
-      ivp = Vector{Cint}(undef, ma57.n)
-      iperm = Vector{Cint}(undef, ma57.n)
-
-      status = 0
-
-      libhsl.$fname(ma57.n::Ref{Cint}, ma57.__fact::Ptr{$T}, ma57.__lfact::Ref{Cint}, ma57.__ifact::Ptr{Cint},
-                    ma57.__lifact::Ref{Cint}, nebdu::Ref{Cint}, nzl::Ref{Cint}, ipl::Ptr{Cint}, irn::Ptr{Cint},
-                    fl::Ptr{$T}, nzd::Ref{Cint}, ipd::Ptr{Cint}, id::Ptr{Cint}, d::Ptr{$T}, ivp::Ptr{Cint},
-                    iperm::Ptr{Cint}, ma57.info.rank::Ref{Cint}, status::Ref{Cint})::Cvoid
-
-      status < 0 && throw(Ma57Exception("Ma57: Error while retrieving factors", status))
-      s = ma57.control.icntl[15] == 1 ? ma57.__fact[(end - ma57.n):(end - 1)] : ones($T, ma57.n)
-      L = SparseMatrixCSC(ma57.n, ma57.n, ipl, irn, fl)
-      D = SparseMatrixCSC(ma57.n, ma57.n, ipd, id, d)
-
-      return (L, D, s, iperm)
-    end
-  end
-end
-
-## alter D ---------------------------------------------------------------------
-function ma57_alter_d(ma57::Ma57{T}, d::Matrix{T}) where {T <: Ma57Data}
-  ka = 1
-  k2 = ma57.__ifact[1]
-  kd = 0
-  kw = 4
-  for blk = 1:abs(ma57.__ifact[3])
-    ncols = ma57.__ifact[kw]
-    nrows = ma57.__ifact[kw + 1]
-    two = false
-    for i = 1:nrows
-      kd = kd + 1
-      ma57.__fact[ka] = d[1, kd]
-      if ma57.__ifact[kw + 1 + i] < 0
-        two = !two
-      end
-      if two
-        ma57.__fact[k2] = d[2, kd]
-        k2 = k2 + 1
-      else
-        status = ma57.info.info[1]
-        d[2, kd] == 0.0 ||
-          throw(Ma57Exception("Ma57: Erroneous 2x2 block specified in d[2,$kd]", status))
-      end
-      ka = ka + nrows + 1 - i
-    end
-    ka = ka + nrows * (ncols - nrows)
-    kw = kw + ncols + 2
-  end
-end
-
-## additional docstrings -------------------------------------------------------
-"""
-# Retrieve factors, scaling, and permutation.
-
-  L, D, s, iperm = ma57_get_factors(M)
-
-  Function will retrieve a unit triangular matrix L, a block-diagonal matrix D a scaling
-  vector s and a permutation vector p such that
-
-  P * S * A * S * P' = L * D * L'
-
-  where S = spdiagm(s) and P = speye(n)[p,:].
-
-  Note that the numerical factorization must have been performed and have succeeded.
-
-## Input arguments
-
-* `M::Ma57`: factorized `Ma57` object
-
-## Return values
-
-* `L::SparseMatrixCSC{T<:Ma57Data,Int}`: L factor
-* `D::SparseMatrixCSC{T<:Ma57Data,Int}`: D factor
-* `s::Vector{T}`: diagonal components of scaling matrix S
-* `iperm::Vector{Int}`: array representin permutation matrix P
-
-## Example:
-
-```julia
-
-julia> using HSL
-
-julia> T = Float64;
-
-julia> rows = Cint[1, 1, 2, 2, 3, 3, 5]; cols = Cint[1, 2, 3, 5, 3, 4, 5];
-
-julia> vals = T[2, 3, 4, 6, 1, 5, 1];
-
-julia> A = sparse(rows, cols, vals); A = A + triu(A, 1)';
-
-julia> b = T[8, 45, 31, 15, 17]
-
-julia> ϵ = sqrt(eps(eltype(A)))
-
-julia> xexact = T[1, 2, 3, 4, 5]
-
-julia> M = Ma57(A)
-
-julia> ma57_factorize(M)
-
-julia> (L, D, s, p) = ma57_get_factors(M)
-
-julia> S = spdiagm(s)
-
-julia> P = speye(M.n)[p, :]
-
-julia> vecnorm(P * S * A * S * P' - L * D * L') ≤ ϵ * vecnorm(A)
-
-true
-```
-"""
-ma57_get_factors
-
-"""
-# Alter block diagonal matrix `D`
-
-## Input arguments:
-
-* `M::Ma57`: `Ma57` object
-* `F::Matrix{Float64}`: diagonal adjustment
-
-
-## Example:
-
-```julia
-julia> using HSL
-
-julia> T = Float64;
-
-julia> rows = Cint[1, 1, 2, 2, 3, 3, 5]; cols = Cint[1, 2, 3, 5, 3, 4, 5];
-
-julia> vals = T[2, 3, 4, 6, 1, 5, 1];
-
-julia> A = sparse(rows, cols, vals); A = A + triu(A, 1)';
-
-julia> b = T[8, 45, 31, 15, 17]
-
-julia> ϵ = sqrt(eps(eltype(A)))
-
-julia> xexact = T[1, 2, 3, 4, 5]
-
-julia> M = Ma57(A)
-
-julia> ma57_factorize(M)
-
-julia> (L, D, s, p) = ma57_get_factors(M)
-
-julia> d1 = abs.(diag(D))
-
-julia> d2 = [diag(D, 1) ; 0]
-
-julia> F = [full(d1)' ; full(d2)']
-
-julia> ma57_alter_d(M, F)
-```
-"""
-ma57_alter_d
diff --git a/test/runtests.jl b/test/runtests.jl
index cd68902..1e72874 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,6 @@ Random.seed!(666)  # Random tests are diabolical
 
 if LIBHSL_isfunctional()
   include("test_hsl_ma57.jl")
-  # include("test_hsl_ma57_patch.jl")
   include("test_hsl_ma97.jl")
   include("test_kb07.jl")
   include("test_mc21.jl")
diff --git a/test/test_hsl_ma57.jl b/test/test_hsl_ma57.jl
index 1aa796f..232b313 100644
--- a/test/test_hsl_ma57.jl
+++ b/test/test_hsl_ma57.jl
@@ -26,6 +26,52 @@ function test_ma57(A, M, b, xexact)
   @test norm(x_ls - x_ls_star) ≤ ϵ * norm(x_ls_star)
 end
 
+function test_ma57_factors(A, M, b, xexact)
+  ϵ = sqrt(eps(eltype(A)))
+  ma57_factorize!(M)
+
+  # extract factors
+  (L, D, s, p) = ma57_get_factors(M)
+  S = spdiagm(0 => s)
+  P = sparse(1.0I, M.n, M.n)[p, :]
+  L = convert(SparseMatrixCSC{Float64, Int}, L) # Convert to 64 bits because \ is not defined for 32 bits.
+  @test norm(P * S * A * S * P' - L * D * L') ≤ ϵ * norm(A)
+
+  # test partial solves
+  b1 = S \ (P' * (L \ (P * (S * b))))
+  x1 = ma57_solve(M, b, job = :LS)
+  @test norm(x1 - b1) ≤ ϵ * norm(b1)
+
+  b2 = S * (P' * (Matrix(D) \ (P * (S * b))))
+  x2 = ma57_solve(M, b, job = :DS)
+  @test norm(x2 - b2) ≤ ϵ * norm(b2)
+
+  b3 = S * (P' * (L' \ (P * (S \ b))))
+  x3 = ma57_solve(M, b, job = :LPS)
+  @test norm(x3 - b3) ≤ ϵ * norm(b3)
+
+  # Test ma57_alter_d
+  B = inv(Tridiagonal(D))
+  d1 = diag(B)
+  d2 = [diag(B, 1); 0][:]
+  ma57_alter_d(M, [Vector(d1)'; Vector(d2)'])
+  x4 = ma57_solve(M, b)
+  @test norm(x4 - xexact) ≤ ϵ * norm(xexact)
+
+  (L2, D2, s2, p2) = ma57_get_factors(M)
+  @test norm(D - D2) ≤ ϵ
+
+  # alter the D factor
+  Dabs = abs.(D)
+  Babs = inv(Tridiagonal(Dabs))
+  d1 = diag(Babs)
+  d2 = [diag(Babs, 1); 0][:]
+  ma57_alter_d(M, [Vector(d1)'; Vector(d2)'])
+
+  (L3, D3, s3, p3) = ma57_get_factors(M)
+  @test norm(Dabs - D3) ≤ ϵ
+end
+
 @testset "hsl_ma57" begin
   @testset "Data Type: $T" for T in (Float32, Float64)
     n = 5
@@ -60,5 +106,27 @@ end
     b = rand(T, n)
     xexact = lu(A) \ b
     test_ma57(A, M, b, xexact)
+
+    @testset "hsl_ma57_factors" begin
+      rows = Cint[1, 1, 2, 2, 3, 3, 5]
+      cols = Cint[1, 2, 3, 5, 3, 4, 5]
+      vals = T[2, 3, 4, 6, 1, 5, 1]
+      A = sparse(rows, cols, vals)
+      A = A + triu(A, 1)'
+
+      M = ma57_coord(5, rows, cols, vals)
+
+      b = T[8, 45, 31, 15, 17]
+      xexact = T[1, 2, 3, 4, 5]
+      test_ma57_factors(A, M, b, xexact)
+
+      n = 5
+      A = convert(T, 3) * convert(SparseMatrixCSC{T, Cint}, sprand(T, n, n, 0.5))
+      A = A + A' + convert(SparseMatrixCSC{T, Cint}, sparse(T(1) * I, n, n))
+      M = Ma57(A)
+      b = rand(T, n)
+      xexact = lu(A) \ b
+      test_ma57_factors(A, M, b, xexact)
+    end
   end
 end
diff --git a/test/test_hsl_ma57_patch.jl b/test/test_hsl_ma57_patch.jl
deleted file mode 100644
index 179e450..0000000
--- a/test/test_hsl_ma57_patch.jl
+++ /dev/null
@@ -1,69 +0,0 @@
-function test_ma57_patch(A, M, b, xexact)
-  ϵ = sqrt(eps(eltype(A)))
-  ma57_factorize!(M)
-
-  # extract factors
-  (L, D, s, p) = ma57_get_factors(M)
-  S = spdiagm(0 => s)
-  P = sparse(1.0I, M.n, M.n)[p, :]
-  L = convert(SparseMatrixCSC{Float64, Int}, L) # Convert to 64 bits because \ is not defined for 32 bits.
-  @test norm(P * S * A * S * P' - L * D * L') ≤ ϵ * norm(A)
-
-  # test partial solves
-  b1 = S \ (P' * (L \ (P * (S * b))))
-  x1 = ma57_solve(M, b, job = :LS)
-  @test norm(x1 - b1) ≤ ϵ * norm(b1)
-
-  b2 = S * (P' * (Matrix(D) \ (P * (S * b))))
-  x2 = ma57_solve(M, b, job = :DS)
-  @test norm(x2 - b2) ≤ ϵ * norm(b2)
-
-  b3 = S * (P' * (L' \ (P * (S \ b))))
-  x3 = ma57_solve(M, b, job = :LPS)
-  @test norm(x3 - b3) ≤ ϵ * norm(b3)
-
-  # Test ma57_alter_d
-  B = inv(Tridiagonal(D))
-  d1 = diag(B)
-  d2 = [diag(B, 1); 0][:]
-  ma57_alter_d(M, [Vector(d1)'; Vector(d2)'])
-  x4 = ma57_solve(M, b)
-  @test norm(x4 - xexact) ≤ ϵ * norm(xexact)
-
-  (L2, D2, s2, p2) = ma57_get_factors(M)
-  @test norm(D - D2) ≤ ϵ
-
-  # alter the D factor
-  Dabs = abs.(D)
-  Babs = inv(Tridiagonal(Dabs))
-  d1 = diag(Babs)
-  d2 = [diag(Babs, 1); 0][:]
-  ma57_alter_d(M, [Vector(d1)'; Vector(d2)'])
-
-  (L3, D3, s3, p3) = ma57_get_factors(M)
-  @test norm(Dabs - D3) ≤ ϵ
-end
-
-@testset "hsl_ma57_patch" begin
-  @testset "Data Type: $T" for T in (Float32, Float64)
-    rows = Cint[1, 1, 2, 2, 3, 3, 5]
-    cols = Cint[1, 2, 3, 5, 3, 4, 5]
-    vals = T[2, 3, 4, 6, 1, 5, 1]
-    A = sparse(rows, cols, vals)
-    A = A + triu(A, 1)'
-
-    M = ma57_coord(5, rows, cols, vals)
-
-    b = T[8, 45, 31, 15, 17]
-    xexact = T[1, 2, 3, 4, 5]
-    test_ma57_patch(A, M, b, xexact)
-
-    n = 5
-    A = convert(T, 3) * convert(SparseMatrixCSC{T, Cint}, sprand(T, n, n, 0.5))
-    A = A + A' + convert(SparseMatrixCSC{T, Cint}, sparse(T(1) * I, n, n))
-    M = Ma57(A)
-    b = rand(T, n)
-    xexact = lu(A) \ b
-    test_ma57_patch(A, M, b, xexact)
-  end
-end