Interface MA27

Project.toml

@@ -9,7 +9,14 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OpenBLAS32_jll = "656ef2d0-ae68-5445-9ca0-591084a874a2"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
+[weakdeps]
+Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
+
+[extensions]
+HSLQuadmathExt = "Quadmath"
+
 [compat]
+Quadmath = "0.5.10"
 OpenBLAS32_jll = "0.3.9"
 julia = "^1.6.0"
 

ext/HSLQuadmathExt.jl

@@ -0,0 +1,11 @@
+module HSLQuadmathExt
+
+using HSL
+using Quadmath
+using LinearAlgebra
+using SparseArrays
+import HSL.libhsl_subset
+
+include("Quadmath/wrappers.jl")
+
+end

ext/Quadmath/wrappers.jl

@@ -0,0 +1,161 @@
+function HSL.ma27iq(icntl, cntl)
+  @ccall libhsl_subset.ma27iq_(icntl::Ptr{Cint}, cntl::Ptr{Float128})::Cvoid
+end
+
+function HSL.ma27aq(n, nz, irn, icn, iw, liw, ikeep, iw1, nsteps, iflag, icntl, cntl, info, ops)
+  @ccall libhsl_subset.ma27aq_(n::Ref{Cint}, nz::Ref{Cint}, irn::Ptr{Cint}, icn::Ptr{Cint}, iw::Ptr{Cint},
+                        liw::Ref{Cint}, ikeep::Ptr{Cint}, iw1::Ptr{Cint}, nsteps::Ref{Cint},
+                        iflag::Ref{Cint}, icntl::Ptr{Cint}, cntl::Ptr{Float128}, info::Ptr{Cint},
+                        ops::Ref{Float128})::Cvoid
+end
+
+function HSL.ma27bq(n, nz, irn, icn, a, la, iw, liw, ikeep, nsteps, maxfrt, iw1, icntl, cntl, info)
+  @ccall libhsl_subset.ma27bq_(n::Ref{Cint}, nz::Ref{Cint}, irn::Ptr{Cint}, icn::Ptr{Cint},
+                        a::Ptr{Float128},
+                        la::Ref{Cint}, iw::Ptr{Cint}, liw::Ref{Cint}, ikeep::Ptr{Cint},
+                        nsteps::Ref{Cint}, maxfrt::Ref{Cint}, iw1::Ptr{Cint}, icntl::Ptr{Cint},
+                        cntl::Ptr{Float128}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27cq(n, a, la, iw, liw, w, maxfrt, rhs, iw1, nsteps, icntl, info)
+  @ccall libhsl_subset.ma27cq_(n::Ref{Cint}, a::Ptr{Float128}, la::Ref{Cint}, iw::Ptr{Cint},
+                        liw::Ref{Cint},
+                        w::Ptr{Float128}, maxfrt::Ref{Cint}, rhs::Ptr{Float128}, iw1::Ptr{Cint},
+                        nsteps::Ref{Cint}, icntl::Ptr{Cint}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27gq(n, nz, irn, icn, iw, lw, ipe, iq, flag, iwfr, icntl, info)
+  @ccall libhsl_subset.ma27gq_(n::Ref{Cint}, nz::Ref{Cint}, irn::Ptr{Cint}, icn::Ptr{Cint}, iw::Ptr{Cint},
+                        lw::Ref{Cint}, ipe::Ptr{Cint}, iq::Ptr{Cint}, flag::Ptr{Cint},
+                        iwfr::Ref{Cint}, icntl::Ptr{Cint}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27hq(n, ipe, iw, lw, iwfr, nv, nxt, lst, ipd, flag, iovflo, ncmpa, fratio)
+  @ccall libhsl_subset.ma27hq_(n::Ref{Cint}, ipe::Ptr{Cint}, iw::Ptr{Cint}, lw::Ref{Cint}, iwfr::Ref{Cint},
+                        nv::Ptr{Cint}, nxt::Ptr{Cint}, lst::Ptr{Cint}, ipd::Ptr{Cint},
+                        flag::Ptr{Cint}, iovflo::Ref{Cint}, ncmpa::Ref{Cint},
+                        fratio::Ref{Float128})::Cvoid
+end
+
+function HSL.ma27uq(n, ipe, iw, lw, iwfr, ncmpa)
+  @ccall libhsl_subset.ma27uq_(n::Ref{Cint}, ipe::Ptr{Cint}, iw::Ptr{Cint}, lw::Ref{Cint}, iwfr::Ref{Cint},
+                        ncmpa::Ref{Cint})::Cvoid
+end
+
+function HSL.ma27jq(n, nz, irn, icn, perm, iw, lw, ipe, iq, flag, iwfr, icntl, info)
+  @ccall libhsl_subset.ma27jq_(n::Ref{Cint}, nz::Ref{Cint}, irn::Ptr{Cint}, icn::Ptr{Cint},
+                        perm::Ptr{Cint},
+                        iw::Ptr{Cint}, lw::Ref{Cint}, ipe::Ptr{Cint}, iq::Ptr{Cint},
+                        flag::Ptr{Cint},
+                        iwfr::Ref{Cint}, icntl::Ptr{Cint}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27kq(n, ipe, iw, lw, iwfr, ips, ipv, nv, flag, ncmpa)
+  @ccall libhsl_subset.ma27kq_(n::Ref{Cint}, ipe::Ptr{Cint}, iw::Ptr{Cint}, lw::Ref{Cint}, iwfr::Ref{Cint},
+                        ips::Ptr{Cint}, ipv::Ptr{Cint}, nv::Ptr{Cint}, flag::Ptr{Cint},
+                        ncmpa::Ref{Cint})::Cvoid
+end
+
+function HSL.ma27lq(n, ipe, nv, ips, ne, na, nd, nsteps, nemin)
+  @ccall libhsl_subset.ma27lq_(n::Ref{Cint}, ipe::Ptr{Cint}, nv::Ptr{Cint}, ips::Ptr{Cint}, ne::Ptr{Cint},
+                        na::Ptr{Cint}, nd::Ptr{Cint}, nsteps::Ref{Cint}, nemin::Ref{Cint})::Cvoid
+end
+
+function HSL.ma27mq(n, nz, irn, icn, perm, na, ne, nd, nsteps, lstki, lstkr, iw, info, ops)
+  @ccall libhsl_subset.ma27mq_(n::Ref{Cint}, nz::Ref{Cint}, irn::Ptr{Cint}, icn::Ptr{Cint},
+                        perm::Ptr{Cint},
+                        na::Ptr{Cint}, ne::Ptr{Cint}, nd::Ptr{Cint}, nsteps::Ref{Cint},
+                        lstki::Ptr{Cint}, lstkr::Ptr{Cint}, iw::Ptr{Cint}, info::Ptr{Cint},
+                        ops::Ref{Float128})::Cvoid
+end
+
+function HSL.ma27nq(n, nz, nz1, a, la, irn, icn, iw, liw, perm, iw2, icntl, info)
+  @ccall libhsl_subset.ma27nq_(n::Ref{Cint}, nz::Ref{Cint}, nz1::Ref{Cint}, a::Ptr{Float128},
+                        la::Ref{Cint},
+                        irn::Ptr{Cint}, icn::Ptr{Cint}, iw::Ptr{Cint}, liw::Ref{Cint},
+                        perm::Ptr{Cint}, iw2::Ptr{Cint}, icntl::Ptr{Cint}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27oq(n, nz, a, la, iw, liw, perm, nstk, nsteps, maxfrt, nelim, iw2, icntl, cntl, info)
+  @ccall libhsl_subset.ma27oq_(n::Ref{Cint}, nz::Ref{Cint}, a::Ptr{Float128}, la::Ref{Cint}, iw::Ptr{Cint},
+                        liw::Ref{Cint}, perm::Ptr{Cint}, nstk::Ptr{Cint}, nsteps::Ref{Cint},
+                        maxfrt::Ref{Cint}, nelim::Ptr{Cint}, iw2::Ptr{Cint}, icntl::Ptr{Cint},
+                        cntl::Ptr{Float128}, info::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27pq(a, iw, j1, j2, itop, ireal, ncmpbr, ncmpbi)
+  @ccall libhsl_subset.ma27pq_(a::Ptr{Float128}, iw::Ptr{Cint}, j1::Ref{Cint}, j2::Ref{Cint},
+                        itop::Ref{Cint}, ireal::Ref{Cint}, ncmpbr::Ref{Cint},
+                        ncmpbi::Ref{Cint})::Cvoid
+end
+
+function HSL.ma27qq(n, a, la, iw, liw, w, maxfnt, rhs, iw2, nblk, latop, icntl)
+  @ccall libhsl_subset.ma27qq_(n::Ref{Cint}, a::Ptr{Float128}, la::Ref{Cint}, iw::Ptr{Cint},
+                        liw::Ref{Cint},
+                        w::Ptr{Float128}, maxfnt::Ref{Cint}, rhs::Ptr{Float128}, iw2::Ptr{Cint},
+                        nblk::Ref{Cint}, latop::Ref{Cint}, icntl::Ptr{Cint})::Cvoid
+end
+
+function HSL.ma27rq(n, a, la, iw, liw, w, maxfnt, rhs, iw2, nblk, latop, icntl)
+  @ccall libhsl_subset.ma27rq_(n::Ref{Cint}, a::Ptr{Float128}, la::Ref{Cint}, iw::Ptr{Cint},
+                        liw::Ref{Cint},
+                        w::Ptr{Float128}, maxfnt::Ref{Cint}, rhs::Ptr{Float128}, iw2::Ptr{Cint},
+                        nblk::Ref{Cint}, latop::Ref{Cint}, icntl::Ptr{Cint})::Cvoid
+end
+
+
+for (iname, aname, bname, cname, T) in ((:ma27iq, :ma27aq, :ma27bq, :ma27cq, :Float128),)
+  @eval begin
+    function HSL.Ma27(A::SparseMatrixCSC{$T})
+      n = LinearAlgebra.checksquare(A)
+      nz = nnz(A)
+      irn, icn, a = findnz(A)
+      la = 4 * nz
+      liw = 2 * (2 * nz + 3 * n + 1)
+      iw = Vector{Cint}(undef, liw)
+      w = Vector{$T}(undef, 0)
+      ikeep = Vector{Cint}(undef, 3 * n)
+      iw1 = Vector{Cint}(undef, 2 * n)
+      nsteps = Ref{Cint}(0)
+      maxfrt = Ref{Cint}(1)
+      iflag = 0
+      icntl = Vector{Cint}(undef, 30)
+      cntl = Vector{$T}(undef, 5)
+      info = Vector{Cint}(undef, 20)
+      ops = Ref{$T}(0)
+      HSL.$iname(icntl, cntl)
+      return HSL.Ma27{$T}(n, nz, Cint.(irn), Cint.(icn), a, la, iw, liw, w, ikeep, iw1,
+                     nsteps, maxfrt, iflag, icntl, cntl, info, ops)
+    end
+
+    function HSL.ma27_update!(F::HSL.Ma27{$T}, A::SparseMatrixCSC{$T})
+      copyto!(F.a, A.nzval) 
+      return F
+    end
+
+    function HSL.ma27_analyze!(F::HSL.Ma27{$T})
+      HSL.$aname(F.n, F.nz, F.irn, F.icn, F.iw, F.liw, F.ikeep, F.iw1, F.nsteps, F.iflag, F.icntl, F.cntl, F.info, F.ops)
+      return F
+    end
+
+    function HSL.ma27_factorize!(F::HSL.Ma27{$T})
+      HSL.$bname(F.n, F.nz, F.irn, F.icn, F.a, F.la, F.iw, F.liw, F.ikeep, F.nsteps, F.maxfrt, F.iw1, F.icntl, F.cntl, F.info)
+      resize!(F.w, F.maxfrt[])
+      return F
+    end
+
+    function HSL.ma27_solve!(F::HSL.Ma27{$T}, x::Vector{$T}, b::Vector{$T})
+      copyto!(x, b)
+      HSL.$cname(F.n, F.a, F.la, F.iw, F.liw, F.w, F.maxfrt, x, F.iw1, F.nsteps, F.icntl, F.info)
+      return b
+    end
+
+    function HSL.ma27(A::SparseMatrixCSC{$T}, b::Vector{$T})
+      F = HSL.Ma27(A)
+      HSL.ma27_analyze!(F)
+      HSL.ma27_factorize!(F)
+      x = Vector{$T}(undef, A.n)
+      HSL.ma27_solve!(F, x, b)
+    end
+  end
+end

gen/analyzer.jl

@@ -2,8 +2,8 @@
 using HSL_jll
 using JuliaFormatter
 
-release = "2023.11.7"
-libhsl = "/home/alexis/Bureau/git/hsl/libhsl/libHSL-$release/"
+release = "2024.11.8"
+libhsl = "/home/alexis/Bureau/git/hsl/libhsl/libHSL.v$release/"
 
 # Symbols of the shared library libhsl
 symbols_path = "symbols.txt"

gen/wrapper.jl

@@ -11,7 +11,7 @@ function wrapper(name::String, headers::Vector{String}, optimized::Bool)
   @info "Wrapping $name"
 
   cd(@__DIR__)
-  include_dir = joinpath(HSL_jll.artifact_dir, "include")
+  include_dir = joinpath(HSL_jll.artifact_dir, "include", "libhsl")
   options = load_options(joinpath(@__DIR__, "hsl.toml"))
   options["general"]["output_file_path"] = joinpath("..", "src", "C", "$(name).jl")
 
@@ -47,7 +47,7 @@ function hsl_headers(include::String, package::String, precisions::Vector{Char})
 end
 
 function main(name::String="all"; optimized::Bool=false)
-  include = joinpath(HSL_jll.artifact_dir, "include")
+  include = joinpath(HSL_jll.artifact_dir, "include", "libhsl")
 
   if name == "all" || name == "libhsl"
     wrapper("libhsl", [joinpath(include, "libhsl.h")], optimized)

src/Fortran/ma27.jl

@@ -192,3 +192,21 @@ function ma27r(n, a, la, iw, liw, w, maxfnt, rhs, iw2, nblk, latop, icntl)
                        w::Ptr{Float32}, maxfnt::Ref{Cint}, rhs::Ptr{Float32}, iw2::Ptr{Cint},
                        nblk::Ref{Cint}, latop::Ref{Cint}, icntl::Ptr{Cint})::Cvoid
 end
+
+# Routines available in the extension HSLQuadmathExt.jl
+function ma27iq end
+function ma27aq end
+function ma27bq end
+function ma27cq end
+function ma27gq end
+function ma27hq end
+function ma27uq end
+function ma27jq end
+function ma27kq end
+function ma27lq end
+function ma27mq end
+function ma27nq end
+function ma27oq end
+function ma27pq end
+function ma27qq end
+function ma27rq end

src/HSL.jl

@@ -3,13 +3,19 @@ module HSL
 using Libdl
 using LinearAlgebra
 using SparseArrays
+import LinearAlgebra: checksquare
 
 if haskey(ENV, "JULIA_HSL_LIBRARY_PATH")
   const libhsl = joinpath(ENV["JULIA_HSL_LIBRARY_PATH"], "libhsl.$dlext")
+  const libhsl_subset = joinpath(ENV["JULIA_HSL_LIBRARY_PATH"], "libhsl_subset.$dlext")
+  const libhsl_subset_64 = joinpath(ENV["JULIA_HSL_LIBRARY_PATH"], "libhsl_subset_64.$dlext")
   const HSL_INSTALLATION = "CUSTOM"
 else
   using OpenBLAS32_jll
-  using HSL_jll
+  import HSL_jll
+  const libhsl = HSL_jll.libhsl
+  const libhsl_subset = HSL_jll.libhsl_subset
+  const libhsl_subset_64 = HSL_jll.libhsl_subset_64
   const HSL_INSTALLATION = "ARTIFACT"
 end
 
@@ -37,6 +43,7 @@ include("wrappers.jl")
 include("hsl_ma57.jl")
 include("hsl_ma97.jl")
 include("kb07.jl")
+include("ma27.jl")
 include("mc21.jl")
 include("mc77.jl")
 

src/ma27.jl

@@ -0,0 +1,79 @@
+export ma27
+
+mutable struct Ma27{T}
+  n::Cint
+  nz::Cint
+  irn::Vector{Cint}
+  icn::Vector{Cint}
+  a::Vector{T}
+  la::Cint
+  iw::Vector{Cint}
+  liw::Cint
+  w::Vector{T}
+  ikeep::Vector{Cint}
+  iw1::Vector{Cint}
+  nsteps::Ref{Cint}
+  maxfrt::Ref{Cint}
+  iflag::Cint
+  icntl::Vector{Cint}
+  cntl::Vector{T}
+  info::Vector{Cint}
+  ops::Ref{T}
+end
+
+for (iname, aname, bname, cname, T) in ((:ma27i , :ma27a , :ma27b , :ma27c , :Float32),
+                                        (:ma27id, :ma27ad, :ma27bd, :ma27cd, :Float64))
+  @eval begin
+    function Ma27(A::SparseMatrixCSC{$T})
+      n = checksquare(A)
+      nz = nnz(A)
+      irn, icn, a = findnz(A)
+      la = 4 * nz
+      liw = 2 * (2 * nz + 3 * n + 1)
+      iw = Vector{Cint}(undef, liw)
+      w = Vector{$T}(undef, 0)
+      ikeep = Vector{Cint}(undef, 3 * n)
+      iw1 = Vector{Cint}(undef, 2 * n)
+      nsteps = Ref{Cint}(0)
+      maxfrt = Ref{Cint}(1)
+      iflag = 0
+      icntl = Vector{Cint}(undef, 30)
+      cntl = Vector{$T}(undef, 5)
+      info = Vector{Cint}(undef, 20)
+      ops = Ref{$T}(0)
+      $iname(icntl, cntl)
+      return HSL.Ma27{$T}(n, nz, Cint.(irn), Cint.(icn), a, la, iw, liw, w, ikeep, iw1,
+                          nsteps, maxfrt, iflag, icntl, cntl, info, ops)
+    end
+
+    function ma27_update!(F::Ma27{$T}, A::SparseMatrixCSC{$T})
+      copyto!(F.a, A.nzval) 
+      return F
+    end
+
+    function ma27_analyze!(F::Ma27{$T})
+      $aname(F.n, F.nz, F.irn, F.icn, F.iw, F.liw, F.ikeep, F.iw1, F.nsteps, F.iflag, F.icntl, F.cntl, F.info, F.ops)
+      return F
+    end
+
+    function ma27_factorize!(F::Ma27{$T})
+      $bname(F.n, F.nz, F.irn, F.icn, F.a, F.la, F.iw, F.liw, F.ikeep, F.nsteps, F.maxfrt, F.iw1, F.icntl, F.cntl, F.info)
+      resize!(F.w, F.maxfrt[])
+      return F
+    end
+
+    function ma27_solve!(F::Ma27{$T}, x::Vector{$T}, b::Vector{$T})
+      copyto!(x, b)
+      $cname(F.n, F.a, F.la, F.iw, F.liw, F.w, F.maxfrt, x, F.iw1, F.nsteps, F.icntl, F.info)
+      return b
+    end
+
+    function ma27(A::SparseMatrixCSC{$T}, b::Vector{$T})
+      F = Ma27(A)
+      ma27_analyze!(F)
+      ma27_factorize!(F)
+      x = Vector{$T}(undef, A.n)
+      ma27_solve!(F, x, b)
+    end
+  end
+end

test/runtests.jl

@@ -10,6 +10,7 @@ if LIBHSL_isfunctional()
   include("test_hsl_ma57.jl")
   include("test_hsl_ma97.jl")
   include("test_kb07.jl")
+  include("test_ma27.jl")
   include("test_mc21.jl")
   include("test_mc77.jl")
 else

test/test_ma27.jl

@@ -0,0 +1,21 @@
+@testset "ma27" begin
+  A = [2 3 0 0 0;
+       3 0 4 0 6;
+       0 4 1 5 0;
+       0 0 5 0 0;
+       0 6 0 0 1] |> sparse
+  b = [8; 45; 31; 15; 17]
+  xstar = [1; 2; 3; 4; 5]
+
+  T = Float32
+  A2 = T.(A)
+  b2 = T.(b)
+  x = ma27(A2, b2)
+  @test x ≈ xstar
+
+  T = Float64
+  A2 = T.(A)
+  b2 = T.(b)
+  x = ma27(A2, b2)
+  @test x ≈ xstar
+end