JuliaSmoothOptimizers · amontoison · May 22, 2025 · Jul 5, 2025
diff --git a/gen/Project.toml b/gen/Project.toml
@@ -4,5 +4,5 @@ HSL_jll = "017b0a0e-03f4-516a-9b91-836bbd1904dd"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 
 [compat]
-julia = "1.6"
-HSL_jll = "=2024.12.10"
+julia = "1.9"
+HSL_jll = "2025.1.19"
diff --git a/gen/analyzer.jl b/gen/analyzer.jl
@@ -4,7 +4,8 @@ using JuliaFormatter
 
 # libHSL
 name_hsl_package = "libhsl"
-release = "2024.11.28"
+# release = "2024.11.28"
+release = "2025.6.26"
 libhsl = "/home/alexis/Bureau/git/hsl/libhsl/libHSL.v$release/"
 
 # HSL_SUBSET
@@ -32,6 +33,29 @@ hsl_precision = Dict{Char, String}('i' => "integer",
                                    'c' => "complex",
                                    'z' => "double_complex")
 
+julia_to_c = Dict{String, String}("Cvoid" => "void ",
+                                  "UInt8" => "char ",
+                                  "Ref{UInt8}" => "char *",
+                                  "Ptr{UInt8}" => "char *",
+                                  "Float32" => "float ",
+                                  "Ref{Float32}" => "float *",
+                                  "Ptr{Float32}" => "float *",
+                                  "Float64" => "double ",
+                                  "Ref{Float64}" => "double *",
+                                  "Ptr{Float64}" => "double *",
+                                  "Cint" => "int ",
+                                  "Ref{Cint}" => "int *",
+                                  "Ptr{Cint}" => "int *",
+                                  "Int64" => "long ",
+                                  "Ref{Int64}" => "long *",
+                                  "Ptr{Int64}" => "long *",
+                                  "ComplexF32" => "float _Complex ",
+                                  "Ref{ComplexF32}" => "float _Complex *",
+                                  "Ptr{ComplexF32}" => "float _Complex *",
+                                  "ComplexF64" => "double _Complex ",
+                                  "Ref{ComplexF64}" => "double _Complex *",
+                                  "Ptr{ComplexF64}" => "double _Complex *",)
+
 # Function to easily get the extension of a file
 function file_extension(file::String)
   pos_dot = findlast(==('.'), file)
@@ -444,6 +468,9 @@ function main(name::String="all"; verbose::Bool=false, f90::Bool=false, print_in
   # We use it to update src/wrappers.jl
   list_include = String[]
 
+  # path_header = joinpath(@__DIR__, "headers", "libhsl", "libhsl.h")
+  # file_header = open(path_header, "w")
+
   for (root, dirs, files) in walkdir(libhsl)
 
     # We don't want to go inside "examples", metis" and "libhsl" folders
@@ -463,6 +490,11 @@ function main(name::String="all"; verbose::Bool=false, f90::Bool=false, print_in
 
         path_wrapper = joinpath(@__DIR__, "..", "src", "Fortran", name_hsl_package, "$(package).jl") |> normpath
         file_wrapper = open(path_wrapper, "w")
+        path_header = joinpath(@__DIR__, "headers", "libhsl", "$(package).h")
+        file_header = open(path_header, "w")
+
+        # write(file_header, '\n')
+        write(file_header, "/* C interface to Fortran routines from $(uppercase(package)) */\n")
 
         # Debug mode (also replace `package == name` by `'/' ∉ package`)
         # path_wrapper = joinpath("..", "src", "Fortran", "debug.jl")
@@ -503,8 +535,14 @@ function main(name::String="all"; verbose::Bool=false, f90::Bool=false, print_in
 
             if name_hsl_package == "libhsl"
               (fname ∉ symbols) && @warn "Unable to find the symbol $fname in the shared library libhsl."
-              write(file_wrapper, "function $signature\n")
-              write(file_wrapper, "  @ccall libhsl.$fname(")
+              if output_type ∉ ("Cvoid", "Cint", "Float32", "Float64")
+                signature = replace(signature, "$(fname[1:end-1])(" => "$(fname[1:end-1])(result,")
+                write(file_wrapper, "function $signature\n")
+                write(file_wrapper, "  @ccall libhsl.$fname(result::Ref{$(output_type)}, ")
+              else
+                write(file_wrapper, "function $signature\n")
+                write(file_wrapper, "  @ccall libhsl.$fname(")
+              end
               for k = 1:narguments
                 if types[k] == ""
                   format = false
@@ -526,9 +564,27 @@ function main(name::String="all"; verbose::Bool=false, f90::Bool=false, print_in
                 format = false
                 @warn "Unable to determine the output type"
               end
-              write(file_wrapper, ")::$(output_type)\n")
+              if output_type ∉ ("Cvoid", "Cint", "Float32", "Float64")
+                write(file_wrapper, ")::Cvoid\n")
+              else
+                write(file_wrapper, ")::$(output_type)\n")
+              end
               write(file_wrapper, "end\n")
 
+              # C Headers
+              if !("Ptr{Ptr{UInt8}}" in types)
+                if output_type ∉ ("Cvoid", "Cint", "Float32", "Float64")
+                  write(file_header, "void $fname($(julia_to_c[output_type])*result, ")
+                else
+                  write(file_header, "$(julia_to_c[output_type])$fname(")
+                end
+                for k = 1:narguments
+                  write(file_header, "$(julia_to_c[types[k]])$(arguments[k])")
+                  (k < narguments) && write(file_header, ", ")
+                end
+                write(file_header, ");\n")
+              end
+
             elseif name_hsl_package == "hsl_subset"
               if endswith(fname, "r_")
                 variant = 0
@@ -586,12 +642,15 @@ function main(name::String="all"; verbose::Bool=false, f90::Bool=false, print_in
             index < num_fnames && write(file_wrapper, "\n")
           end
         end
+        close(file_header)
         close(file_wrapper)
         format && format_file(path_wrapper, YASStyle())
       end
     end
   end
 
+  # close(file_header)
+
   if print_include
     for str in list_include
       println(str)

diff --git a/src/Fortran/libhsl/pa16.jl b/src/Fortran/libhsl/pa16.jl
@@ -13,8 +13,9 @@ function pa16cd(a, n, root, e, w, f, ig, cr)
                         w::Ptr{ComplexF64}, f::Ptr{Float64}, ig::Ptr{Cint}, cr::Ptr{Float64})::Cvoid
 end
 
-function pa16dd(z, n, a)
-  @ccall libhsl.pa16dd_(z::Ref{ComplexF64}, n::Ref{Cint}, a::Ptr{ComplexF64})::ComplexF64
+function pa16dd(result, z, n, a)
+  @ccall libhsl.pa16dd_(result::Ref{ComplexF64}, z::Ref{ComplexF64}, n::Ref{Cint},
+                        a::Ptr{ComplexF64})::Cvoid
 end
 
 function pa16a(a, n, root, e, w, f, ig)
@@ -32,6 +33,7 @@ function pa16c(a, n, root, e, w, f, ig, cr)
                        w::Ptr{ComplexF32}, f::Ptr{Float32}, ig::Ptr{Cint}, cr::Ptr{Float32})::Cvoid
 end
 
-function pa16d(z, n, a)
-  @ccall libhsl.pa16d_(z::Ref{ComplexF32}, n::Ref{Cint}, a::Ptr{ComplexF32})::ComplexF32
+function pa16d(result, z, n, a)
+  @ccall libhsl.pa16d_(result::Ref{ComplexF32}, z::Ref{ComplexF32}, n::Ref{Cint},
+                       a::Ptr{ComplexF32})::Cvoid
 end
diff --git a/src/Fortran/libhsl/pa17.jl b/src/Fortran/libhsl/pa17.jl
@@ -13,6 +13,7 @@ function pa17cd(a, n, root, e, w, ig, cr)
                         w::Ptr{Float64}, ig::Ptr{Cint}, cr::Ptr{Float64})::Cvoid
 end
 
-function pa17dd(z, n, a)
-  @ccall libhsl.pa17dd_(z::Ref{ComplexF64}, n::Ref{Cint}, a::Ptr{Float64})::ComplexF64
+function pa17dd(result, z, n, a)
+  @ccall libhsl.pa17dd_(result::Ref{ComplexF64}, z::Ref{ComplexF64}, n::Ref{Cint},
+                        a::Ptr{Float64})::Cvoid
 end
diff --git a/src/HSL.jl b/src/HSL.jl
@@ -46,6 +46,7 @@ include("hsl_mc64.jl")
 include("kb07.jl")
 include("mc19.jl")
 include("mc21.jl")
+include("mc26.jl")
 include("mc29.jl")
 include("mc30.jl")
 include("mc77.jl")

diff --git a/src/mc26.jl b/src/mc26.jl
@@ -0,0 +1,60 @@
+export mc26
+
+"""
+    jptr, iptr, ata = mc26(A::SparseMatrixCSC{Float32,Cint})
+    jptr, iptr, ata = mc26(A::SparseMatrixCSC{Float64,Cint})
+
+Compute the structure and optionally the values of `AᵀA`.
+
+Returns:
+- `jptr`: column pointer array for the compressed sparse column (CSC) format of `AᵀA`.
+- `iptr`: row indices corresponding to each nonzero entry in `AᵀA`.
+- `ata`: numerical values of `AᵀA` corresponding to the `(iptr, jptr)` structure.
+"""
+function mc26 end
+
+for (mc26ir, mc26ar, T) in ((:mc26i , :mc26a , :Float32),
+                            (:mc26id, :mc26ad, :Float64))
+  @eval begin
+    function mc26(A::SparseMatrixCSC{$T,Cint})
+      icntl = zeros(Cint, 5)
+      info = zeros(Cint, 5)
+      $mc26ir(icntl)  # Set default control parameters
+
+      m, n = size(A)
+      itype = Cint(1)  # CSC format
+      nz = nnz(A)
+      a = A.nzval
+      ind = A.rowval
+
+      lata = max(1, nz + n)  # initial guess
+      iata = Vector{Cint}(undef, lata)
+      ata  = Vector{$T}(undef, lata)
+      jptr = Vector{Cint}(undef, n + 1)
+      iptr = Vector{Cint}(undef, lata)
+      iw = zeros(Cint, max(1, 2 * n))
+      nzata = Ref{Cint}(0)
+      success = false
+
+      while !success
+        $mc26ar(m, n, itype, nz, a, ind, lata, iata, ata, jptr, iptr, nzata, iw, icntl, info)
+
+        if info[1] == -4
+          lata = max(lata + n, iata[1])
+          resize!(iata, lata)
+          resize!(ata,  lata)
+          resize!(iptr, lata)
+        elseif info[1] == 0
+          success = true
+        else
+          error("mc26a failed with info[1] = $(info[1])")
+        end
+      end
+
+      resize!(ata, nzata[])
+      resize!(iptr, nzata[])
+
+      return jptr, iptr, ata
+    end
+  end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -13,6 +13,7 @@ if LIBHSL_isfunctional()
   include("test_kb07.jl")
   include("test_mc19.jl")
   include("test_mc21.jl")
+  # include("test_mc26.jl")
   include("test_mc29.jl")
   include("test_mc30.jl")
   include("test_mc77.jl")

diff --git a/test/test_mc26.jl b/test/test_mc26.jl
@@ -0,0 +1,19 @@
+@testset "mc26" begin
+    # Define the matrix A:
+    # [1.0 4.0]
+    # [0.0 3.0]
+    # [5.0 0.0]
+    rowval = Cint[1, 2, 3, 1]  # rows of nonzeros (Fortran-style indexing)
+    colptr = Cint[1, 3, 5]     # column pointer (CSC)
+    nzval  = [1.0, 3.0, 5.0, 4.0]  # nonzero values
+
+    @testset "Precision $T" for T in (Float32, Float64)
+        A = SparseMatrixCSC{T,Cint}(3, 2, colptr, rowval, T.(nzval))
+
+        jptr, iptr, ata = mc26(A)
+
+        @test jptr == Cint[1, 3, 4]   # two columns in AᵀA, jptr[3] = 4 means 3 nnz
+        @test iptr == Cint[1, 2, 2]   # row indices for columns
+        @test ata ≈ T[26.0, 4.0, 25.0] # AᵀA = [26 4; 4 25], stored by columns
+    end
+end