From 004cf138781eca59b9374c0b90035947634c08ff Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Wed, 12 Jan 2022 18:28:30 -0500
Subject: [PATCH 01/16] Use concrete references

---
 src/hsl_ma97.jl | 60 +++++++++++++++++++++++++------------------------
 1 file changed, 31 insertions(+), 29 deletions(-)

diff --git a/src/hsl_ma97.jl b/src/hsl_ma97.jl
index adf91ef..2c160bf 100644
--- a/src/hsl_ma97.jl
+++ b/src/hsl_ma97.jl
@@ -82,13 +82,13 @@ mutable struct Ma97_Control{T <: Ma97Real}
                   zeros(Cint, 5), zeros(T, 10))
 
     if T == Float32
-      ccall((:ma97_default_control_s, libhsl_ma97), Nothing, (Ref{Ma97_Control},), control)
+      ccall((:ma97_default_control_s, libhsl_ma97), Nothing, (Ref{Ma97_Control{Float32}},), control)
     elseif T == Float64
-      ccall((:ma97_default_control_d, libhsl_ma97), Nothing, (Ref{Ma97_Control},), control)
+      ccall((:ma97_default_control_d, libhsl_ma97), Nothing, (Ref{Ma97_Control{Float64}},), control)
     elseif T == ComplexF32
-      ccall((:ma97_default_control_c, libhsl_ma97), Nothing, (Ref{Ma97_Control},), control)
+      ccall((:ma97_default_control_c, libhsl_ma97), Nothing, (Ref{Ma97_Control{Float32}},), control)
     elseif T == ComplexF64
-      ccall((:ma97_default_control_z, libhsl_ma97), Nothing, (Ref{Ma97_Control},), control)
+      ccall((:ma97_default_control_z, libhsl_ma97), Nothing, (Ref{Ma97_Control{Float64}},), control)
     end
     control.f_arrays = 1  # Use 1-based indexing for arrays, avoiding copies.
     control.print_level = print_level
@@ -250,10 +250,10 @@ for (fname, typ) in ((:ma97_finalise_s, Float32),
                      (:ma97_finalise_d, Float64),
                      (:ma97_finalise_c, ComplexF32),
                      (:ma97_finalise_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
-    function ma97_finalize(ma97 :: Ma97{$typ, $(data_map[typ])})
+    function ma97_finalize(ma97 :: Ma97{$typ, $S})
       ccall(($(string(fname)), libhsl_ma97), Nothing,
             (Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}),
              ma97.__akeep,   ma97.__fkeep)
@@ -267,17 +267,17 @@ for (fname, freename, typ) in ((:ma97_analyse_s, :ma97_free_akeep_s, Float32),
                                (:ma97_analyse_d, :ma97_free_akeep_d, Float64),
                                (:ma97_analyse_c, :ma97_free_akeep_c, ComplexF32),
                                (:ma97_analyse_z, :ma97_free_akeep_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
     function ma97_csc(n :: Int, colptr :: Vector{Ti}, rowval :: Vector{Ti}, nzval :: Vector{$typ}; kwargs...) where {Ti <: Integer}
-      control = Ma97_Control{$(data_map[typ])}(; kwargs...)
-      info = Ma97_Info{$(data_map[typ])}()
-      M = Ma97{$typ, $(data_map[typ])}([convert(Ptr{Nothing}, C_NULL)], [convert(Ptr{Nothing}, C_NULL)], n, colptr, rowval, nzval, control, info)
+      control = Ma97_Control{$S}(; kwargs...)
+      info = Ma97_Info{$S}()
+      M = Ma97{$typ, $S}([convert(Ptr{Nothing}, C_NULL)], [convert(Ptr{Nothing}, C_NULL)], n, colptr, rowval, nzval, control, info)
 
       # Perform symbolic analysis.
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{Cint}),
+            (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{Cint}),
              1,    M.n,  M.colptr,  M.rowval,  C_NULL,    M.__akeep,      M.control,         M.info,         C_NULL)
 
       if M.info.flag < 0
@@ -322,18 +322,18 @@ for (fname, freename, typ) in ((:ma97_analyse_coord_s, :ma97_free_akeep_s, Float
                                (:ma97_analyse_coord_d, :ma97_free_akeep_d, Float64),
                                (:ma97_analyse_coord_c, :ma97_free_akeep_c, ComplexF32),
                                (:ma97_analyse_coord_z, :ma97_free_akeep_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
     function ma97_coord(n :: Int, cols :: Vector{Ti}, rows :: Vector{Ti}, nzval :: Vector{$typ}; kwargs...) where {Ti <: Integer}
-      control = Ma97_Control{$(data_map[typ])}(; kwargs...)
-      info = Ma97_Info{$(data_map[typ])}()
-      M = Ma97{$typ, $(data_map[typ])}([convert(Ptr{Nothing}, C_NULL)], [convert(Ptr{Nothing}, C_NULL)], n, convert(Vector{Cint}, cols), convert(Vector{Cint}, rows), nzval, control, info)
+      control = Ma97_Control{$S}(; kwargs...)
+      info = Ma97_Info{$S}()
+      M = Ma97{$typ, $S}([convert(Ptr{Nothing}, C_NULL)], [convert(Ptr{Nothing}, C_NULL)], n, convert(Vector{Cint}, cols), convert(Vector{Cint}, rows), nzval, control, info)
       nz = length(cols)
 
       # Perform symbolic analysis.
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{Cint}),
+            (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{Cint}),
              M.n,  nz,   M.rowval,  M.colptr,  C_NULL,    M.__akeep,      M.control,         M.info,         C_NULL)
 
       if M.info.flag < 0
@@ -353,14 +353,14 @@ for (fname, typ) in ((:ma97_factor_s, Float32),
                      (:ma97_factor_d, Float64),
                      (:ma97_factor_c, ComplexF32),
                      (:ma97_factor_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
-    function ma97_factorize!(ma97 :: Ma97{$typ, $(data_map[typ])}; matrix_type :: Symbol=:real_indef)
+    function ma97_factorize!(ma97 :: Ma97{$typ, $S}; matrix_type :: Symbol=:real_indef)
       t = matrix_types97[matrix_type]
 
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ},  Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{$(data_map[typ])}),
+            (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ},  Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{Cdouble}),
              t,    C_NULL,    C_NULL,    ma97.nzval, ma97.__akeep,   ma97.__fkeep,   ma97.control,      ma97.info,   C_NULL)
 
       if ma97.info.flag < 0
@@ -395,16 +395,16 @@ for (fname, typ) in ((:ma97_solve_s, Float32),
                      (:ma97_solve_d, Float64),
                      (:ma97_solve_c, ComplexF32),
                      (:ma97_solve_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
-    function ma97_solve!(ma97 :: Ma97{$typ, $(data_map[typ])}, b :: Array{$typ}; job :: Symbol=:A)
+    function ma97_solve!(ma97 :: Ma97{$typ, $S}, b :: Array{$typ}; job :: Symbol=:A)
       size(b, 1) == ma97.n || throw(Ma97Exception("Ma97: rhs size mismatch", 0))
       nrhs = size(b, 2)
 
       j = jobs97[job]
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Cint, Cint, Ptr{$typ}, Cint,   Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}),
+            (Cint, Cint, Ptr{$typ}, Cint,   Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}),
              j,    nrhs, b,         ma97.n, ma97.__akeep,   ma97.__fkeep,   ma97.control,      ma97.info)
 
       if ma97.info.flag < 0
@@ -437,6 +437,7 @@ for (fname, typ) in ((:ma97_factor_solve_s, Float32),
                      (:ma97_factor_solve_c, ComplexF32),
                      (:ma97_factor_solve_z, ComplexF64))
 
+  S = data_map[typ]
   @eval begin
 
     function ma97_solve!(A :: SparseMatrixCSC{$typ,Int}, b :: Array{$typ}; matrix_type :: Symbol=:real_indef)
@@ -445,7 +446,7 @@ for (fname, typ) in ((:ma97_factor_solve_s, Float32),
       size(b, 1) == M.n || throw(Ma97Exception("Ma97: rhs size mismatch", 0))
       nrhs = size(b, 2)
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Cint, Ptr{$typ}, Cint, Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{$(data_map[typ])}),
+            (Cint, Ptr{Cint}, Ptr{Cint}, Ptr{$typ}, Cint, Ptr{$typ}, Cint, Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{Cdouble}),
              t,    M.colptr,  M.rowval,  M.nzval,   nrhs, b,         M.n,  M.__akeep,      M.__fkeep,      M.control,      M.info,      C_NULL)
 
       if M.info.flag < 0
@@ -466,9 +467,10 @@ for (indef, posdef, typ) in ((:ma97_enquire_indef_s, :ma97_enquire_posdef_s, Flo
                              (:ma97_enquire_indef_c, :ma97_enquire_posdef_c, ComplexF32),
                              (:ma97_enquire_indef_z, :ma97_enquire_posdef_z, ComplexF64))
 
+  S = data_map[typ]
   @eval begin
 
-    function ma97_inquire(ma97 :: Ma97{$typ, $(data_map[typ])}; matrix_type :: Symbol=:real_indef)
+    function ma97_inquire(ma97 :: Ma97{$typ, $S}; matrix_type :: Symbol=:real_indef)
       if matrix_type in [:real_indef, :herm_indef, :cmpl_indef]
         piv_order = zeros(Cint, ma97.n)
         # AMBUSH ALERT: although Julia will call the C interface of the library
@@ -476,13 +478,13 @@ for (indef, posdef, typ) in ((:ma97_enquire_indef_s, :ma97_enquire_posdef_s, Flo
         # documentation says d should be n x 2, we must declare 2 x n.
         d = zeros($typ, 2, ma97.n)
         ccall(($(string(indef)), libhsl_ma97), Nothing,
-              (Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{Cint}, Ptr{$typ}),
+              (Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{Cint}, Ptr{$typ}),
                ma97.__akeep,   ma97.__fkeep,   ma97.control,      ma97.info,      piv_order, d)
         ret = (piv_order, d)
       else
         d = zeros($typ, ma97.n)
         ccall(($(string(posdef)), libhsl_ma97), Nothing,
-              (Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}, Ptr{$typ}),
+              (Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}, Ptr{$typ}),
                ma97.__akeep,   ma97.__fkeep,   ma97.control,      ma97.info,      d)
         ret = d
       end
@@ -505,14 +507,14 @@ for (fname, typ) in ((:ma97_alter_s, Float32),
                      (:ma97_alter_d, Float64),
                      (:ma97_alter_c, ComplexF32),
                      (:ma97_alter_z, ComplexF64))
-
+  S = data_map[typ]
   @eval begin
 
-    function ma97_alter!(ma97 :: Ma97{$typ, $(data_map[typ])}, d :: Array{$typ, 2})
+    function ma97_alter!(ma97 :: Ma97{$typ, $S}, d :: Array{$typ, 2})
       n, m = size(d)
       (m == ma97.n && n == 2) || throw(Ma97Exception("Ma97: input array d must be n x 2", 0))
       ccall(($(string(fname)), libhsl_ma97), Nothing,
-            (Ptr{$typ}, Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control}, Ref{Ma97_Info}),
+            (Ptr{$typ}, Ptr{Ptr{Nothing}}, Ptr{Ptr{Nothing}}, Ref{Ma97_Control{$S}}, Ref{Ma97_Info{$S}}),
              d,         ma97.__akeep,   ma97.__fkeep,   ma97.control,      ma97.info)
 
       if ma97.info.flag < 0

From 38c502d26c77327d75ea0ca04bf9848684c3d2d4 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Thu, 13 Jan 2022 10:15:47 -0500
Subject: [PATCH 02/16] Test allocations

---
 test/test_ma97.jl | 16 +++++++++++++++-
 1 file changed, 15 insertions(+), 1 deletion(-)

diff --git a/test/test_ma97.jl b/test/test_ma97.jl
index e36b03e..1c9d5a7 100644
--- a/test/test_ma97.jl
+++ b/test/test_ma97.jl
@@ -1,4 +1,15 @@
-for T in (Float32, Float64, ComplexF32, ComplexF64)
+function getallocs(A)
+    T = eltype(A)
+    n = size(A,1)
+    ma97 = Ma97(A)
+    matrix_type = T in (ComplexF32, ComplexF64) ? :herm_indef : :real_indef
+    allocs = @allocated ma97_factorize!(ma97, matrix_type=matrix_type)
+    b = rand(T, n)
+    x = copy(b)
+    allocs += @allocated ma97_solve!(ma97, x)
+end
+
+for T in (Float32, Float64, ComplexF32, ComplexF64 )
 
   @info "Testing hsl_ma97 with $T data"
   matrix_type = T in (ComplexF32, ComplexF64) ? :herm_indef : :real_indef
@@ -43,6 +54,9 @@ for T in (Float32, Float64, ComplexF32, ComplexF64)
   x = ma97_solve(A, b, matrix_type=matrix_type)
   @test norm(x - x_exact) ≤ ϵ * norm(x_exact)
 
+  # Test allocations
+  @test getallocs(A) == 0
+
   # Test symmetric positive definite A.
   A = A * A';
   matrix_type = T in (ComplexF32, ComplexF64) ? :herm_pd : :real_spd

From 210d1c9638fc6216d4c97b6681f94d9b19a06099 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Fri, 14 Jan 2022 13:40:07 -0500
Subject: [PATCH 03/16] Add ma86 to build

---
 deps/build.jl          | 14 ++++++++++++--
 deps/build_hsl_ma86.jl | 11 +++++++++++
 src/HSL.jl             |  4 ++++
 3 files changed, 27 insertions(+), 2 deletions(-)
 create mode 100644 deps/build_hsl_ma86.jl

diff --git a/deps/build.jl b/deps/build.jl
index 29cabe3..5d0bfc9 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -16,15 +16,19 @@ const hsl_ma57_patch = joinpath(hsl_ma57_path, "get_factors.patch")
 
 const hsl_ma97_version = "2.6.0"
 const hsl_ma97_sha256 = "be5fe822674be93e3d2e1a7d7ed6c5ad831b91cf8ca5150beb473f67af5fcb66"
-
 const hsl_ma97_path = haskey(ENV, "HSL_MA97_PATH") ? ENV["HSL_MA97_PATH"] : joinpath(@__DIR__, "downloads")
 const hsl_ma97_archive = joinpath(hsl_ma97_path, "hsl_ma97-$hsl_ma97_version.tar.gz")
 
+const hsl_ma86_version = "1.7.0"
+const hsl_ma86_sha256 = "845c65bee0bf31507d7b99c87773997012b7b16434da349ad5c361ceb257191e"
+const hsl_ma86_path = haskey(ENV, "HSL_MA86_PATH") ? ENV["HSL_MA86_PATH"] : joinpath(@__DIR__, "downloads")
+const hsl_ma86_archive = joinpath(hsl_ma86_path, "hsl_ma86-$hsl_ma86_version.tar.gz")
+
 const so         = Sys.isapple() ? "dylib" : "so"
 const all_load   = Sys.isapple() ? "-all_load" : "--whole-archive"
 const noall_load = Sys.isapple() ? "-noall_load" : "--no-whole-archive"
 
-hsl_archives = [hsl_ma57_archive, hsl_ma97_archive]
+hsl_archives = [hsl_ma57_archive, hsl_ma97_archive, hsl_ma86_archive]
 
 if any(isfile.(hsl_archives))
   products = Product[]
@@ -45,6 +49,12 @@ if any(isfile.(hsl_archives))
     include("build_hsl_ma57.jl")
   end
 
+  if isfile(hsl_ma86_archive)
+    @info "building ma86"
+    push!(products, FileProduct(prefix, "lib/libhsl_ma86.$so", :libhsl_ma86))
+    include("build_hsl_ma86.jl")
+  end
+
   if isfile(hsl_ma97_archive)
     @info "building ma97"
     push!(products, FileProduct(prefix, "lib/libhsl_ma97.$so", :libhsl_ma97))
diff --git a/deps/build_hsl_ma86.jl b/deps/build_hsl_ma86.jl
new file mode 100644
index 0000000..0ad7224
--- /dev/null
+++ b/deps/build_hsl_ma86.jl
@@ -0,0 +1,11 @@
+open(hsl_ma86_archive) do f
+  if bytes2hex(sha256(f)) != hsl_ma86_sha256
+    error("SHA256 of HSL MA86 doesn't match")
+  end
+end
+run(`tar -zxf $hsl_ma86_archive -C $builddir`)
+cd("$builddir/hsl_ma86-$hsl_ma86_version")
+run(`./configure --prefix=$usrdir F77=gfortran CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --with-blas="-L$libopenblas_dir -lopenblas" --with-metis="-L$libmetis_dir -lmetis"`)
+run(`make install`)
+run(`gfortran -fPIC -shared -Wl,$all_load $libdir/libhsl_ma86.a -L$libopenblas_dir -lopenblas -L$libmetis_dir -lmetis -lgomp -Wl,$noall_load -o $libdir/libhsl_ma86.$so`)
+cd(@__DIR__)
\ No newline at end of file
diff --git a/src/HSL.jl b/src/HSL.jl
index 8063ba7..b0b9e1e 100644
--- a/src/HSL.jl
+++ b/src/HSL.jl
@@ -21,6 +21,7 @@ const data_map = Dict{Type, Type}(Float32 => Cfloat,
                                   ComplexF32 => Cfloat,
                                   ComplexF64 => Cdouble)
 
+
 # package-specific definitions
 if (@isdefined libhsl_ma57) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma57.jl")
@@ -31,5 +32,8 @@ end
 if (@isdefined libhsl_ma97) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma97.jl")
 end
+if (@isdefined libhsl_ma86) || haskey(ENV, "DOCUMENTER_KEY")
+  include("hsl_ma86.jl")
+end
 
 end

From 9809d127952eb74872b87e043e40175f4be1d812 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Mon, 24 Jan 2022 18:01:30 -0500
Subject: [PATCH 04/16] Basic ma86 wrapper

---
 src/hsl_ma86.jl | 196 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 196 insertions(+)
 create mode 100644 src/hsl_ma86.jl

diff --git a/src/hsl_ma86.jl b/src/hsl_ma86.jl
new file mode 100644
index 0000000..125b8ad
--- /dev/null
+++ b/src/hsl_ma86.jl
@@ -0,0 +1,196 @@
+
+const Ma86Data = Union{Float32, Float64, ComplexF32, ComplexF64}
+const Ma86Real = Union{Cfloat, Cdouble}
+
+function appendtype(fname, T)
+  typesuffix = Dict(Float32=>"s", Float64=>"d", ComplexF32=>"c", ComplexF64=>"z")
+  return string(fname) * "_" * typesuffix[T]
+end
+
+mutable struct Ma86_Control{S <: Ma86Real}
+  f_arrays::Cint
+  diagnostics_level::Cint
+  unit_diagnostics::Cint
+  unit_error::Cint
+  unit_warning::Cint
+  nemin::Cint
+  nb::Cint
+  action::Cint
+  nbi::Cint
+  bool_size::Cint
+  small_::S
+  static_::S
+  u::S
+  umin::S
+  scaling::Cint
+  function Ma86_Control(::Type{T}) where T
+    S = data_map[T]
+    control = new{S}(0,0,0,0,0,0,0,0,0,0,zero(S),zero(S),zero(S),zero(S),0)
+    ma86_default_control(T,control)
+    control.f_arrays = 1  # Use 1-based indexing for arrays, avoiding copies.
+    return control
+  end
+end
+Ma86_Control{S}() where S <: Ma86Real = Ma86_Control(S)
+
+mutable struct Ma86_Info{S <: Ma86Real}
+  detlog::S
+  detsign::Cint
+  flag::Cint
+  matrix_rank::Cint
+  maxdepth::Cint
+  num_delay::Cint
+  num_factor::Clong
+  num_flops::Clong
+  num_neg::Cint
+  num_nodes::Cint
+  num_nothresh::Cint
+  num_perturbed::Cint
+  num_two::Cint
+  pool_size::Cint
+  stat::Cint
+  usmall::S
+  function Ma86_Info(::Type{T}) where T <: Ma86Data
+    S = data_map[T]
+    new{S}(
+      zero(S), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Clong), 
+      zero(Clong), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), 
+      zero(Cint), zero(S)
+    )
+  end
+end
+Ma86_Info{S}() where S <: Ma86Real = Ma86_Info(S)
+
+struct Keep
+  ptr::Vector{Ptr{Cvoid}}
+end
+@inline Base.cconvert(::Type{Ref{Ptr{Cvoid}}}, akeep::Keep) = akeep.ptr
+@inline Base.unsafe_convert(::Type{Ptr{Ptr{Cvoid}}}, akeep::Keep) = 
+  Base.unsafe_convert(Ptr{Ptr{Cvoid}}, akeep.ptr)
+Keep() = Keep([C_NULL])
+isnull(akeep::Keep) = akeep.ptr[1] == C_NULL
+
+for T in (Float32, Float64, ComplexF32, ComplexF64) 
+  S = data_map[T]
+
+  @eval function ma86_default_control(::Type{$T}, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_default_control, T)), libhsl_ma86), Cvoid,
+      (Ref{Ma86_Control{$S}},),
+      control
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_analyse(::Type{$T}, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control, 
+                              info::Ma86_Info)
+    @assert length(ptr) == n+1
+    @assert ptr[end] == length(row) + 1
+    @assert length(order) == n
+    ccall(($(appendtype(:ma86_analyse, T)), libhsl_ma86), Cvoid,
+      (Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}, 
+      Ref{Ma86_Info{$S}}),
+      n, ptr, row, order, keep, control, info
+    )
+  end
+
+  if T <: Ma86Real
+    @eval function ma86_factor(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                              info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                                    val::Vector{$T}, order::Vector{Cint}, keep::Keep, 
+                                    control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                    nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                    scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  else
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                               row::Vector{Cint}, val::Vector{$T}, 
+                               order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                               info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                                     row::Vector{Cint}, val::Vector{$T}, 
+                                     order::Vector{Cint}, keep::Keep, 
+                                     control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                     nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                     scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  end
+
+  @eval function ma86_solve(job::Cint, nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                          order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                          info::Ma86_Info{$S})
+    n = size(x,1)
+    @assert job ∈ 0:4
+    @assert nrhs >= 1
+    @assert ldx >= n
+    @assert size(x,2) == nrhs
+    ccall(($(appendtype(:ma86_solve, T)), libhsl_ma86), Cvoid,
+      (Cint, Cint, Cint, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}},
+      Ref{Ma86_Info{$S}}, Ptr{$S}),
+      job, nrhs, ldx, x, order, keep, control, info, C_NULL 
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_finalise(::Type{$T}, keep::Keep, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_finalise, T)), libhsl_ma86), Cvoid,
+      (Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}),
+      keep, control
+    )
+  end
+
+end

From 83f5c21c6a409fcfa7a7a70a21934b0316e755ea Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:47:17 -0500
Subject: [PATCH 05/16] Add ma86 wrapper

---
 src/hsl_ma86.jl | 239 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 237 insertions(+), 2 deletions(-)

diff --git a/src/hsl_ma86.jl b/src/hsl_ma86.jl
index 125b8ad..f54bbad 100644
--- a/src/hsl_ma86.jl
+++ b/src/hsl_ma86.jl
@@ -1,3 +1,4 @@
+export Ma86, Ma86_Control, Ma86_Info
 
 const Ma86Data = Union{Float32, Float64, ComplexF32, ComplexF64}
 const Ma86Real = Union{Cfloat, Cdouble}
@@ -7,6 +8,12 @@ function appendtype(fname, T)
   return string(fname) * "_" * typesuffix[T]
 end
 
+"""
+    Ma86_Control
+
+A simple wrapper around the C control structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
 mutable struct Ma86_Control{S <: Ma86Real}
   f_arrays::Cint
   diagnostics_level::Cint
@@ -33,6 +40,12 @@ mutable struct Ma86_Control{S <: Ma86Real}
 end
 Ma86_Control{S}() where S <: Ma86Real = Ma86_Control(S)
 
+"""
+    Ma86_Info
+
+A simple wrapper around the C info structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
 mutable struct Ma86_Info{S <: Ma86Real}
   detlog::S
   detsign::Cint
@@ -61,6 +74,9 @@ mutable struct Ma86_Info{S <: Ma86Real}
 end
 Ma86_Info{S}() where S <: Ma86Real = Ma86_Info(S)
 
+##############################
+# Convenience Wrapper Types 
+##############################
 struct Keep
   ptr::Vector{Ptr{Cvoid}}
 end
@@ -70,6 +86,223 @@ end
 Keep() = Keep([C_NULL])
 isnull(akeep::Keep) = akeep.ptr[1] == C_NULL
 
+const matrix_types86 = (herm_indef = Cint(-4), cmpl_indef = Cint(-5))
+const jobs86 = (
+  A    = Cint(0),  # Solve Ax = b
+  PL   = Cint(1),  # Solve PLx = b
+  D    = Cint(2),  # Solve Dx = b
+  PLt  = Cint(3),  # Solve (PL)ᴴ x = b
+  DPLt = Cint(4),  # Solve D(PL)ᴴ x = b
+)
+
+Base.@kwdef mutable struct Ma86_Flags
+  isanalysisdone::Bool = false
+  isfactordone::Bool = false
+end
+
+"""
+    Ma86{T,S}
+
+Type for the HSL Ma86 solver, a multithreaded solver for symmetric indefinite or 
+symmetric positive definite systems.
+
+# Constructors
+
+    Ma86(colptr, rowval, nzval; [control, info, order, scale, doanalyse])
+
+Creates a new Ma86 solver with the lower-triangular matrix data stored in CSC format. 
+Both `colptr` and `rowval` will be converted to vectors of `Int32`. The element type of 
+the nonzero values `nzval` will determine the data type used by the solver. The data types
+of `control`, `info`, and `scale` should match the corresponding real data type.
+
+    Ma86(A::SparseMatrix; kwargs...)
+
+Extracts the CSC data from the sparse matrix `A`. If `A` is not lower triangular, the lower
+triangular data will be pulled out automatically. The keyword arguments are identical those
+of the previous constructor.
+
+# Usage
+With a constructed `Ma86` object, you can use any of the following methods to solve the 
+linear system `Ax = b` where `b` can have any number of columns. 
+
+## 1. Call HSL steps independently
+This mirrors the C interface of Ma86. Call the following methods in the specified order:
+
+    HSL.factor!(ma86)
+    HSL.solve!(ma86, b)
+
+Or
+
+    HSL.factorsolve!(ma86, b)
+
+Note the methods above will overwrite the data in `b` with the solution `x`. Alternatively,
+you can use [`solve`](@ref) and ][`factorsolve`](@ref) which will return the solution 
+vector `x` and leave `b` unchanged.
+
+By default, the constructor will run the [`analyse!`](@ref) method on the provided sparsity 
+structure. This step can be skipped by passing `doanalyse=false` to the constructor. The 
+analysis step must then be called manually before calling either `factor!` or 
+`factorsolve!`.
+
+## 2. Use backslash or `ldiv!`
+To return the solution matrix, use `x = ma86 \\ b`, or to overwrite the `b` matrix use 
+`ldiv!(ma86, b)`.
+"""
+mutable struct Ma86{T<:Ma86Data, S<:Ma86Real}
+  __keep::Keep
+  n::Cint
+  colptr::Vector{Cint}
+  rowval::Vector{Cint}
+  nzval::Vector{T}
+  control::Ma86_Control{S}
+  info::Ma86_Info{S}
+  order::Vector{Cint}
+  scale::Vector{S}
+  flags::Ma86_Flags
+  function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}, 
+                control::Ma86_Control{S}, info::Ma86_Info{S}, order::Vector{<:Integer}, 
+                scale::Vector{<:Real}; doanalyse::Bool=true) where {T <: Ma86Data, S}
+    S == data_map[T] || throw(TypeError(:Ma86, "MA86{$T, $S}\n", data_map[T], T))
+    keep = Keep()
+    n = length(colptr) - 1
+    @assert length(order) == n
+    @assert length(scale) == n
+    flags = Ma86_Flags()
+    ma86 = new{T,S}(keep, n, colptr, rowval, nzval, control, info, order, scale, flags)
+    if doanalyse
+      analyse!(ma86)
+    end
+    finalizer(_finalize, ma86)
+  end
+end
+function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}; 
+              control::Ma86_Control=Ma86_Control(T), 
+              info::Ma86_Info=Ma86_Info(T), 
+              order::Vector{<:Integer}=collect(1:length(colptr)-1),
+              scale::Vector{<:Real}=fill(one(data_map[T]), length(colptr)-1),
+              kwargs...) where T
+  # Assumes the data is already lower triangular
+  Ma86(colptr, rowval, nzval, control, info, order, scale; kwargs...)
+end
+function Ma86(A::SparseMatrixCSC; kwargs...)
+  if !istril(A)
+    A = tril(A)
+  end
+  Ma86(A.colptr, A.rowval, A.nzval; kwargs...)
+end
+datatype(::Ma86{T}) where T = T
+realtype(::Ma86{<:Any,S}) where S = S
+
+function _finalize(m::Ma86{T}) where T
+  if !isnull(m.__keep)
+    ma86_finalise(T, m.__keep, m.control)
+  end
+  m.flags.isanalysisdone = false
+  m.flags.isfactordone = false
+  m.__keep.ptr[1] = C_NULL
+  return m
+end
+
+##############################
+# Exception Handling
+##############################
+"""
+   Ma86Exception
+
+Exception thrown by HSL MA86 algorithm. The `flag` field contains the HSL flag 
+for the error. Consult the documentation for details of the exception.
+"""
+mutable struct Ma86Exception <: Exception
+  msg::AbstractString
+  flag::Int
+end
+
+function checkerror(m::Ma86, msg)
+  if m.info.flag != 0
+    _finalize(m)  # free memory in keep before throwing an error
+    throw(Ma86Exception(msg, m.info.flag))
+  end
+end
+
+##############################
+# Convenience Wrapper
+##############################
+
+function analyse!(m::Ma86{T}) where T
+  ma86_analyse(T, m.n, m.colptr, m.rowval, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during symbol analysis.")
+  m.flags.isanalysisdone = true
+  return nothing
+end
+
+function factor!(m::Ma86{<:Ma86Data}; matrix_type::Integer=matrix_types86.herm_indef)
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factor! before analyse!")
+    return nothing
+  end
+  ma86_factor(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    m.scale
+  )
+  checkerror(m, "MA86: Error during factorization.")
+  m.flags.isfactordone = true
+  return nothing
+end
+
+function solve!(m::Ma86{T}, x::VecOrMat{T}; job::Integer=jobs86.A) where T
+  if !m.flags.isanalysisdone || !m.flags.isfactordone
+    error("Cannot call solve! before factor!")
+    return x
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_solve(Cint(job), nrhs, ldx, x, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during solve.")
+  return x
+end
+solve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T = solve!(m, copy(b); kwargs...)
+
+function factorsolve!(m::Ma86{T}, x::VecOrMat{T}; 
+  matrix_type::Integer=matrix_types86.herm_indef
+) where T <: Ma86Data
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factorsolve! before analyse!")
+    return x 
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_factor_solve(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    nrhs, ldx, x, m.scale
+  )
+  checkerror(m, "MA86: Error during factorsolve.")
+  m.flags.isfactordone = true
+  return x
+end
+
+# Out-of-place version
+factorsolve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T <: Ma86Data = 
+  factorsolve!(m, copy(b); kwargs...)
+
+# LinearAlgebra wrappers
+function Base.:\(m::Ma86{T}, b::VecOrMat{T}) where T <: Ma86Data 
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve(m, b)
+  factorsolve(m, b)
+end
+
+function LinearAlgebra.ldiv!(m::Ma86{T}, x::VecOrMat{T}) where T <: Ma86Data
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve!(m, x)
+  factorsolve!(m, x)
+end
+
+##############################
+# C Wrapper
+##############################
+
 for T in (Float32, Float64, ComplexF32, ComplexF64) 
   S = data_map[T]
 
@@ -95,7 +328,9 @@ for T in (Float32, Float64, ComplexF32, ComplexF64)
   end
 
   if T <: Ma86Real
-    @eval function ma86_factor(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
+    # Note that matrix_type argument is unused in these methods
+    # Included to make the signature match for all data types
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
                               order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
                               info::Ma86_Info{$S}, scale::Vector{$S})
       @assert length(ptr) == n+1
@@ -109,7 +344,7 @@ for T in (Float32, Float64, ComplexF32, ComplexF64)
       )
     end
 
-    @eval function ma86_factor_solve(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
                                     val::Vector{$T}, order::Vector{Cint}, keep::Keep, 
                                     control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
                                     nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 

From 956b5281a8dc8fd873104529479a9cdc75ccb45b Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:47:58 -0500
Subject: [PATCH 06/16] Don't include deps in test

---
 test/runtests.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 1fd4aee..dcb23a2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,14 +4,14 @@ using LinearAlgebra
 using SparseArrays
 using Test
 
-include("../deps/deps.jl")
+# include("../deps/deps.jl")
 
-if @isdefined libhsl_ma57
+if isdefined(HSL, :libhsl_ma57)
   include("test_ma57.jl")
-  if @isdefined libhsl_ma57_patch
+  if isdefined(HSL, :libhsl_ma57_patch)
     include("test_ma57_patch.jl")
   end
 end
-if @isdefined libhsl_ma97
+if isdefined(HSL, :libhsl_ma97)
   include("test_ma97.jl")
 end

From 949acfac584ced17c9ca0471e1d2be87b9c78e60 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:48:07 -0500
Subject: [PATCH 07/16] Add test deps

---
 test/Project.toml | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 test/Project.toml

diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..59b801f
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From 51b1395aa4cef242444f00d562e541ada049218f Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:48:37 -0500
Subject: [PATCH 08/16] Add ma86 tests

---
 test/runtests.jl  |   3 +
 test/test_ma86.jl | 178 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 181 insertions(+)
 create mode 100644 test/test_ma86.jl

diff --git a/test/runtests.jl b/test/runtests.jl
index dcb23a2..667edea 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,3 +15,6 @@ end
 if isdefined(HSL, :libhsl_ma97)
   include("test_ma97.jl")
 end
+if isdefined(HSL, :libhsl_ma86)
+  include("test_ma86.jl")
+end
diff --git a/test/test_ma86.jl b/test/test_ma86.jl
new file mode 100644
index 0000000..a4ac338
--- /dev/null
+++ b/test/test_ma86.jl
@@ -0,0 +1,178 @@
+using HSL
+using SparseArrays
+using LinearAlgebra
+using Test
+using Random
+
+function getallocs(A,b)
+    T = eltype(A)
+    n = size(A,1)
+    ma86 = Ma86(A)
+    allocs = @allocated HSL.analyse!(ma86)
+    allocs += @allocated HSL.factor!(ma86)
+    x = copy(b)
+    allocs += @allocated HSL.solve!(ma86, x)
+end
+
+T = Float64
+n = 10
+A = convert(T, 3) * sprand(T, n, n, 0.5)
+Aindef = A + A' + one(T)*I
+ma86 = Ma86(Aindef)
+ma86.flags.isanalysisdone
+HSL.factor!(ma86)
+b = rand(T,n)
+HSL.solve!(ma86, b)
+
+@testset "MA86 ($T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+    # @info "Testing hsl_ma86 with $T data"
+    n = 10
+    A = convert(T, 3) * sprand(T, n, n, 0.5)
+    Aindef = A + A' + one(T)*I
+
+    Apsd = A'A + one(T)*I
+    @test isposdef(Matrix(Apsd))
+
+
+    # Test symmetric indefinite A
+    @testset "$label" for (label,A) in (("Indefinite", Aindef), ("Positive Definite",Apsd))
+        n = size(A,1)
+        T = eltype(A)
+
+        # Test constructors
+        L = SparseMatrixCSC{T,Cint}(tril(A))
+        let ptr = L.colptr, row = L.rowval, val = L.nzval
+            # Check that new arrays aren't allocated if the input in Cint
+            let ma86 = Ma86(ptr, row, val)
+                @test ma86.colptr === ptr
+                @test ma86.rowval === row
+                @test ma86.nzval === val 
+                @test ma86.flags.isanalysisdone == true
+                @test !HSL.isnull(ma86.__keep)
+            end
+
+            # Check the keyword constructor            
+            control = Ma86_Control(T)
+            info = Ma86_Info(T)
+            S = HSL.data_map[T]
+            order = randperm(Cint(n))
+            scale = fill(one(S), n)
+            let ma86 = HSL.Ma86(ptr, row, val, 
+                    control=control, info=info, order=order, scale=scale)
+                @test ma86.order === order
+                @test ma86.scale === scale 
+                @test ma86.control === control
+                @test ma86.info === info 
+            end
+
+            # Pass in Int64 data
+            let ma86 = HSL.Ma86(Int64.(ptr), Int64.(row), val, 
+                    order=randperm(Int64(n)), scale=Float64.(scale))
+                @test ma86.colptr !== ptr
+                @test ma86.rowval !== row
+                @test ma86.order !== order
+                @test ma86.scale !== scale 
+            end
+
+            # Pass in lower trianglar matrix
+            # Make sure it doesn't reallocate new data
+            let ma86 = HSL.Ma86(L) 
+                @test ma86.colptr === L.colptr
+                @test ma86.rowval === L.rowval
+                @test ma86.nzval === L.nzval
+            end
+
+            # Check analyse skipping
+            let ma86 = HSL.Ma86(L, doanalyse=false)
+                @test ma86.flags.isanalysisdone == false
+            end
+        end
+
+        @testset "with $(size(b,2)) rhs vectors" for b in (rand(T, n), rand(T, n, 3))
+            n = size(A,1)
+            T = eltype(A)
+
+            # factor, then solve
+            ma86_0 = Ma86(A)
+            let ma86 = ma86_0
+                x = copy(b)
+                HSL.factor!(ma86)
+                @test ma86.info.flag == 0
+                HSL.solve!(ma86, x)
+                @test ma86.info.flag == 0
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+                @test HSL.solve(ma86, b) ≈ x  # out-of-place solve
+                @test ma86 \ b ≈ x            # backslash
+                x2 = copy(b)
+                ldiv!(ma86, x2)
+                @test x2 ≈ x
+            end
+
+            # in-place
+            ma86_1 = Ma86(A)
+            let ma86 = ma86_1
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # out-of-place
+            ma86_2 = Ma86(A)
+            let ma86 = ma86_2
+                x = HSL.factorsolve(ma86, b)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # Test linear algebra
+            ma86_3 = Ma86(A)
+            x = ma86_3 \ b
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            ma86_4 = Ma86(A)
+            x = copy(b)
+            ldiv!(ma86_4, x)
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            # Test errors
+            ma86_err = Ma86(Apsd, doanalyse=false)
+            let ma86 = ma86_err
+                @test_throws ErrorException HSL.factor!(ma86)
+                b = rand(T, n)
+                @test_throws ErrorException HSL.solve!(ma86, b)
+                @test_throws ErrorException HSL.factorsolve!(ma86, b)
+
+                HSL.analyse!(ma86)
+                @test_throws ErrorException HSL.solve!(ma86, b) 
+
+                b = rand(T, n)
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                x = copy(b)
+                HSL.solve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                b_bad = rand(T, 2n)
+                @test_throws DimensionMismatch HSL.solve!(ma86, b_bad)
+                @test_throws DimensionMismatch HSL.factorsolve!(ma86, b_bad)
+
+                # Test allocations
+                @test getallocs(A, b) == 0
+            end
+            
+            # Test HSL Errors
+            # Prints out annoying warning message from MC78.
+            # Not sure how to avoid this.
+            # Should we even allow for non-square matrices to be passed to the constructor?
+            B = sprand(T, n, n, 0.1)
+            ma86_hslerr = Ma86(B, doanalyse=false)
+            ma86_hslerr.control.diagnostics_level = -1
+            ma86_hslerr.control.unit_diagnostics = -1
+            ma86_hslerr.control.unit_error = -1
+            ma86_hslerr.control.unit_warning = -1
+            @test_throws HSL.Ma86Exception HSL.analyse!(ma86_hslerr)
+        end
+    end
+end
+

From 8859ca6e05193d17c162f6b60e89e8e37790bf07 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Fri, 14 Jan 2022 13:40:07 -0500
Subject: [PATCH 09/16] Add ma86 to build

---
 deps/build.jl          | 17 ++++++++++++++---
 deps/build_hsl_ma86.jl | 11 +++++++++++
 src/HSL.jl             |  4 ++++
 3 files changed, 29 insertions(+), 3 deletions(-)
 create mode 100644 deps/build_hsl_ma86.jl

diff --git a/deps/build.jl b/deps/build.jl
index cc11739..3f142e5 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -56,24 +56,29 @@ const hsl_ma57_patch = joinpath(hsl_ma57_path, "get_factors.patch")
 ##############################
 # MA97
 ##############################
-const hsl_ma97_verions = [
+const hsl_ma97_versions = [
     HSLVersion("hsl_ma97", "2.7.0", "ac3a081d3a28e9ecb8871ce769f4ced2a5ffa5a9c36defbd2c844ae3493ccb37", ".tar.gz"),
     HSLVersion("hsl_ma97", "2.7.0", "8221b607d96554d7a57cc60483c7305ef43a8785dc4171ac2e8da087900a1100", ".zip"),
     HSLVersion("hsl_ma97", "2.6.0", "be5fe822674be93e3d2e1a7d7ed6c5ad831b91cf8ca5150beb473f67af5fcb66", ".tar.gz"),
 ]
 const hsl_ma97_path = haskey(ENV, "HSL_MA97_PATH") ? ENV["HSL_MA97_PATH"] : joinpath(@__DIR__, "downloads")
-hsl_ma97_version = findversion(hsl_ma97_verions, hsl_ma97_path)
+hsl_ma97_version = findversion(hsl_ma97_versions, hsl_ma97_path)
 const hsl_ma97_archive = isnothing(hsl_ma97_version) ? "" : joinpath(hsl_ma97_path, getname(hsl_ma97_version))
 
 ##############################
 # Build
 ##############################
-const hsl_archives = [hsl_ma57_archive, hsl_ma97_archive]
+const hsl_archives = [hsl_ma57_archive, hsl_ma97_archive, hsl_ma86_archive]
 
 const HSL_FC = haskey(ENV, "HSL_FC") ? ENV["HSL_FC"] : "gfortran"
 const HSL_F77 = haskey(ENV, "HSL_F77") ? ENV["HSL_F77"] : HSL_FC
 const HSL_CC = haskey(ENV, "HSL_CC") ? ENV["HSL_CC"] : "gcc"
 
+const hsl_ma86_version = "1.7.0"
+const hsl_ma86_sha256 = "845c65bee0bf31507d7b99c87773997012b7b16434da349ad5c361ceb257191e"
+const hsl_ma86_path = haskey(ENV, "HSL_MA86_PATH") ? ENV["HSL_MA86_PATH"] : joinpath(@__DIR__, "downloads")
+const hsl_ma86_archive = joinpath(hsl_ma86_path, "hsl_ma86-$hsl_ma86_version.tar.gz")
+
 const so         = Sys.isapple() ? "dylib" : "so"
 const all_load   = Sys.isapple() ? "-all_load" : "--whole-archive"
 const noall_load = Sys.isapple() ? "-noall_load" : "--no-whole-archive"
@@ -99,6 +104,12 @@ if any(isfile.(hsl_archives))
     include("build_hsl_ma57.jl")
   end
 
+  if isfile(hsl_ma86_archive)
+    @info "building ma86"
+    push!(products, FileProduct(prefix, "lib/libhsl_ma86.$so", :libhsl_ma86))
+    include("build_hsl_ma86.jl")
+  end
+
   if isfile(hsl_ma97_archive)
     @info "building ma97"
     push!(products, FileProduct(prefix, "lib/libhsl_ma97.$so", :libhsl_ma97))
diff --git a/deps/build_hsl_ma86.jl b/deps/build_hsl_ma86.jl
new file mode 100644
index 0000000..0ad7224
--- /dev/null
+++ b/deps/build_hsl_ma86.jl
@@ -0,0 +1,11 @@
+open(hsl_ma86_archive) do f
+  if bytes2hex(sha256(f)) != hsl_ma86_sha256
+    error("SHA256 of HSL MA86 doesn't match")
+  end
+end
+run(`tar -zxf $hsl_ma86_archive -C $builddir`)
+cd("$builddir/hsl_ma86-$hsl_ma86_version")
+run(`./configure --prefix=$usrdir F77=gfortran CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --with-blas="-L$libopenblas_dir -lopenblas" --with-metis="-L$libmetis_dir -lmetis"`)
+run(`make install`)
+run(`gfortran -fPIC -shared -Wl,$all_load $libdir/libhsl_ma86.a -L$libopenblas_dir -lopenblas -L$libmetis_dir -lmetis -lgomp -Wl,$noall_load -o $libdir/libhsl_ma86.$so`)
+cd(@__DIR__)
\ No newline at end of file
diff --git a/src/HSL.jl b/src/HSL.jl
index bd25aa7..7249f89 100644
--- a/src/HSL.jl
+++ b/src/HSL.jl
@@ -24,6 +24,7 @@ const data_map = Dict{Type, Type}(Float32 => Cfloat,
                                   ComplexF32 => Cfloat,
                                   ComplexF64 => Cdouble)
 
+
 # package-specific definitions
 if (@isdefined libhsl_ma57) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma57.jl")
@@ -34,5 +35,8 @@ end
 if (@isdefined libhsl_ma97) || haskey(ENV, "DOCUMENTER_KEY")
   include("hsl_ma97.jl")
 end
+if (@isdefined libhsl_ma86) || haskey(ENV, "DOCUMENTER_KEY")
+  include("hsl_ma86.jl")
+end
 
 end

From a71a59f79764dd49fa435deedd8948903700836f Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Mon, 24 Jan 2022 18:01:30 -0500
Subject: [PATCH 10/16] Basic ma86 wrapper

---
 src/hsl_ma86.jl | 196 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 196 insertions(+)
 create mode 100644 src/hsl_ma86.jl

diff --git a/src/hsl_ma86.jl b/src/hsl_ma86.jl
new file mode 100644
index 0000000..125b8ad
--- /dev/null
+++ b/src/hsl_ma86.jl
@@ -0,0 +1,196 @@
+
+const Ma86Data = Union{Float32, Float64, ComplexF32, ComplexF64}
+const Ma86Real = Union{Cfloat, Cdouble}
+
+function appendtype(fname, T)
+  typesuffix = Dict(Float32=>"s", Float64=>"d", ComplexF32=>"c", ComplexF64=>"z")
+  return string(fname) * "_" * typesuffix[T]
+end
+
+mutable struct Ma86_Control{S <: Ma86Real}
+  f_arrays::Cint
+  diagnostics_level::Cint
+  unit_diagnostics::Cint
+  unit_error::Cint
+  unit_warning::Cint
+  nemin::Cint
+  nb::Cint
+  action::Cint
+  nbi::Cint
+  bool_size::Cint
+  small_::S
+  static_::S
+  u::S
+  umin::S
+  scaling::Cint
+  function Ma86_Control(::Type{T}) where T
+    S = data_map[T]
+    control = new{S}(0,0,0,0,0,0,0,0,0,0,zero(S),zero(S),zero(S),zero(S),0)
+    ma86_default_control(T,control)
+    control.f_arrays = 1  # Use 1-based indexing for arrays, avoiding copies.
+    return control
+  end
+end
+Ma86_Control{S}() where S <: Ma86Real = Ma86_Control(S)
+
+mutable struct Ma86_Info{S <: Ma86Real}
+  detlog::S
+  detsign::Cint
+  flag::Cint
+  matrix_rank::Cint
+  maxdepth::Cint
+  num_delay::Cint
+  num_factor::Clong
+  num_flops::Clong
+  num_neg::Cint
+  num_nodes::Cint
+  num_nothresh::Cint
+  num_perturbed::Cint
+  num_two::Cint
+  pool_size::Cint
+  stat::Cint
+  usmall::S
+  function Ma86_Info(::Type{T}) where T <: Ma86Data
+    S = data_map[T]
+    new{S}(
+      zero(S), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Clong), 
+      zero(Clong), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), zero(Cint), 
+      zero(Cint), zero(S)
+    )
+  end
+end
+Ma86_Info{S}() where S <: Ma86Real = Ma86_Info(S)
+
+struct Keep
+  ptr::Vector{Ptr{Cvoid}}
+end
+@inline Base.cconvert(::Type{Ref{Ptr{Cvoid}}}, akeep::Keep) = akeep.ptr
+@inline Base.unsafe_convert(::Type{Ptr{Ptr{Cvoid}}}, akeep::Keep) = 
+  Base.unsafe_convert(Ptr{Ptr{Cvoid}}, akeep.ptr)
+Keep() = Keep([C_NULL])
+isnull(akeep::Keep) = akeep.ptr[1] == C_NULL
+
+for T in (Float32, Float64, ComplexF32, ComplexF64) 
+  S = data_map[T]
+
+  @eval function ma86_default_control(::Type{$T}, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_default_control, T)), libhsl_ma86), Cvoid,
+      (Ref{Ma86_Control{$S}},),
+      control
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_analyse(::Type{$T}, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control, 
+                              info::Ma86_Info)
+    @assert length(ptr) == n+1
+    @assert ptr[end] == length(row) + 1
+    @assert length(order) == n
+    ccall(($(appendtype(:ma86_analyse, T)), libhsl_ma86), Cvoid,
+      (Cint, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}, 
+      Ref{Ma86_Info{$S}}),
+      n, ptr, row, order, keep, control, info
+    )
+  end
+
+  if T <: Ma86Real
+    @eval function ma86_factor(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
+                              order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                              info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+                                    val::Vector{$T}, order::Vector{Cint}, keep::Keep, 
+                                    control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                    nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                    scale::Vector{$S})
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  else
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                               row::Vector{Cint}, val::Vector{$T}, 
+                               order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                               info::Ma86_Info{$S}, scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      ccall(($(appendtype(:ma86_factor, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, scale
+      )
+    end
+
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, 
+                                     row::Vector{Cint}, val::Vector{$T}, 
+                                     order::Vector{Cint}, keep::Keep, 
+                                     control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
+                                     nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                                     scale::Vector{$S})
+      @assert matrix_type ∈ (-4,-5)
+      @assert length(ptr) == n+1
+      @assert length(row) + 1 == ptr[end]
+      @assert length(val) + 1 == ptr[end]
+      @assert length(order) == n
+      @assert nrhs >= 1
+      @assert ldx >= n
+      @assert length(scale) == n
+      @assert size(x,1) == n
+      @assert size(x,2) == nrhs
+      ccall(($(appendtype(:ma86_factor_solve, T)), libhsl_ma86), Cvoid,
+        (Cint, Cint, Ref{Cint}, Ref{Cint}, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, 
+        Ref{Ma86_Control{$S}}, Ref{Ma86_Info{$S}}, Cint, Cint, Ref{$T}, Ref{$S}),
+        matrix_type, n, ptr, row, val, order, keep, control, info, nrhs, ldx, x, scale
+      )
+    end
+  end
+
+  @eval function ma86_solve(job::Cint, nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 
+                          order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
+                          info::Ma86_Info{$S})
+    n = size(x,1)
+    @assert job ∈ 0:4
+    @assert nrhs >= 1
+    @assert ldx >= n
+    @assert size(x,2) == nrhs
+    ccall(($(appendtype(:ma86_solve, T)), libhsl_ma86), Cvoid,
+      (Cint, Cint, Cint, Ref{$T}, Ref{Cint}, Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}},
+      Ref{Ma86_Info{$S}}, Ptr{$S}),
+      job, nrhs, ldx, x, order, keep, control, info, C_NULL 
+    )
+  end
+
+  # Add type as argument to avoid ambiguity
+  @eval function ma86_finalise(::Type{$T}, keep::Keep, control::Ma86_Control{$S})
+    ccall(($(appendtype(:ma86_finalise, T)), libhsl_ma86), Cvoid,
+      (Ptr{Ptr{Cvoid}}, Ref{Ma86_Control{$S}}),
+      keep, control
+    )
+  end
+
+end

From 4374c725a9a7a3db2a46135b6d9f0e999f086e89 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:47:17 -0500
Subject: [PATCH 11/16] Add ma86 wrapper

---
 src/hsl_ma86.jl | 239 +++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 237 insertions(+), 2 deletions(-)

diff --git a/src/hsl_ma86.jl b/src/hsl_ma86.jl
index 125b8ad..f54bbad 100644
--- a/src/hsl_ma86.jl
+++ b/src/hsl_ma86.jl
@@ -1,3 +1,4 @@
+export Ma86, Ma86_Control, Ma86_Info
 
 const Ma86Data = Union{Float32, Float64, ComplexF32, ComplexF64}
 const Ma86Real = Union{Cfloat, Cdouble}
@@ -7,6 +8,12 @@ function appendtype(fname, T)
   return string(fname) * "_" * typesuffix[T]
 end
 
+"""
+    Ma86_Control
+
+A simple wrapper around the C control structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
 mutable struct Ma86_Control{S <: Ma86Real}
   f_arrays::Cint
   diagnostics_level::Cint
@@ -33,6 +40,12 @@ mutable struct Ma86_Control{S <: Ma86Real}
 end
 Ma86_Control{S}() where S <: Ma86Real = Ma86_Control(S)
 
+"""
+    Ma86_Info
+
+A simple wrapper around the C info structure for C HSL Ma86 interface. Consult HSL 
+  documentation for more details.
+"""
 mutable struct Ma86_Info{S <: Ma86Real}
   detlog::S
   detsign::Cint
@@ -61,6 +74,9 @@ mutable struct Ma86_Info{S <: Ma86Real}
 end
 Ma86_Info{S}() where S <: Ma86Real = Ma86_Info(S)
 
+##############################
+# Convenience Wrapper Types 
+##############################
 struct Keep
   ptr::Vector{Ptr{Cvoid}}
 end
@@ -70,6 +86,223 @@ end
 Keep() = Keep([C_NULL])
 isnull(akeep::Keep) = akeep.ptr[1] == C_NULL
 
+const matrix_types86 = (herm_indef = Cint(-4), cmpl_indef = Cint(-5))
+const jobs86 = (
+  A    = Cint(0),  # Solve Ax = b
+  PL   = Cint(1),  # Solve PLx = b
+  D    = Cint(2),  # Solve Dx = b
+  PLt  = Cint(3),  # Solve (PL)ᴴ x = b
+  DPLt = Cint(4),  # Solve D(PL)ᴴ x = b
+)
+
+Base.@kwdef mutable struct Ma86_Flags
+  isanalysisdone::Bool = false
+  isfactordone::Bool = false
+end
+
+"""
+    Ma86{T,S}
+
+Type for the HSL Ma86 solver, a multithreaded solver for symmetric indefinite or 
+symmetric positive definite systems.
+
+# Constructors
+
+    Ma86(colptr, rowval, nzval; [control, info, order, scale, doanalyse])
+
+Creates a new Ma86 solver with the lower-triangular matrix data stored in CSC format. 
+Both `colptr` and `rowval` will be converted to vectors of `Int32`. The element type of 
+the nonzero values `nzval` will determine the data type used by the solver. The data types
+of `control`, `info`, and `scale` should match the corresponding real data type.
+
+    Ma86(A::SparseMatrix; kwargs...)
+
+Extracts the CSC data from the sparse matrix `A`. If `A` is not lower triangular, the lower
+triangular data will be pulled out automatically. The keyword arguments are identical those
+of the previous constructor.
+
+# Usage
+With a constructed `Ma86` object, you can use any of the following methods to solve the 
+linear system `Ax = b` where `b` can have any number of columns. 
+
+## 1. Call HSL steps independently
+This mirrors the C interface of Ma86. Call the following methods in the specified order:
+
+    HSL.factor!(ma86)
+    HSL.solve!(ma86, b)
+
+Or
+
+    HSL.factorsolve!(ma86, b)
+
+Note the methods above will overwrite the data in `b` with the solution `x`. Alternatively,
+you can use [`solve`](@ref) and ][`factorsolve`](@ref) which will return the solution 
+vector `x` and leave `b` unchanged.
+
+By default, the constructor will run the [`analyse!`](@ref) method on the provided sparsity 
+structure. This step can be skipped by passing `doanalyse=false` to the constructor. The 
+analysis step must then be called manually before calling either `factor!` or 
+`factorsolve!`.
+
+## 2. Use backslash or `ldiv!`
+To return the solution matrix, use `x = ma86 \\ b`, or to overwrite the `b` matrix use 
+`ldiv!(ma86, b)`.
+"""
+mutable struct Ma86{T<:Ma86Data, S<:Ma86Real}
+  __keep::Keep
+  n::Cint
+  colptr::Vector{Cint}
+  rowval::Vector{Cint}
+  nzval::Vector{T}
+  control::Ma86_Control{S}
+  info::Ma86_Info{S}
+  order::Vector{Cint}
+  scale::Vector{S}
+  flags::Ma86_Flags
+  function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}, 
+                control::Ma86_Control{S}, info::Ma86_Info{S}, order::Vector{<:Integer}, 
+                scale::Vector{<:Real}; doanalyse::Bool=true) where {T <: Ma86Data, S}
+    S == data_map[T] || throw(TypeError(:Ma86, "MA86{$T, $S}\n", data_map[T], T))
+    keep = Keep()
+    n = length(colptr) - 1
+    @assert length(order) == n
+    @assert length(scale) == n
+    flags = Ma86_Flags()
+    ma86 = new{T,S}(keep, n, colptr, rowval, nzval, control, info, order, scale, flags)
+    if doanalyse
+      analyse!(ma86)
+    end
+    finalizer(_finalize, ma86)
+  end
+end
+function Ma86(colptr::Vector{<:Integer}, rowval::Vector{<:Integer}, nzval::Vector{T}; 
+              control::Ma86_Control=Ma86_Control(T), 
+              info::Ma86_Info=Ma86_Info(T), 
+              order::Vector{<:Integer}=collect(1:length(colptr)-1),
+              scale::Vector{<:Real}=fill(one(data_map[T]), length(colptr)-1),
+              kwargs...) where T
+  # Assumes the data is already lower triangular
+  Ma86(colptr, rowval, nzval, control, info, order, scale; kwargs...)
+end
+function Ma86(A::SparseMatrixCSC; kwargs...)
+  if !istril(A)
+    A = tril(A)
+  end
+  Ma86(A.colptr, A.rowval, A.nzval; kwargs...)
+end
+datatype(::Ma86{T}) where T = T
+realtype(::Ma86{<:Any,S}) where S = S
+
+function _finalize(m::Ma86{T}) where T
+  if !isnull(m.__keep)
+    ma86_finalise(T, m.__keep, m.control)
+  end
+  m.flags.isanalysisdone = false
+  m.flags.isfactordone = false
+  m.__keep.ptr[1] = C_NULL
+  return m
+end
+
+##############################
+# Exception Handling
+##############################
+"""
+   Ma86Exception
+
+Exception thrown by HSL MA86 algorithm. The `flag` field contains the HSL flag 
+for the error. Consult the documentation for details of the exception.
+"""
+mutable struct Ma86Exception <: Exception
+  msg::AbstractString
+  flag::Int
+end
+
+function checkerror(m::Ma86, msg)
+  if m.info.flag != 0
+    _finalize(m)  # free memory in keep before throwing an error
+    throw(Ma86Exception(msg, m.info.flag))
+  end
+end
+
+##############################
+# Convenience Wrapper
+##############################
+
+function analyse!(m::Ma86{T}) where T
+  ma86_analyse(T, m.n, m.colptr, m.rowval, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during symbol analysis.")
+  m.flags.isanalysisdone = true
+  return nothing
+end
+
+function factor!(m::Ma86{<:Ma86Data}; matrix_type::Integer=matrix_types86.herm_indef)
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factor! before analyse!")
+    return nothing
+  end
+  ma86_factor(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    m.scale
+  )
+  checkerror(m, "MA86: Error during factorization.")
+  m.flags.isfactordone = true
+  return nothing
+end
+
+function solve!(m::Ma86{T}, x::VecOrMat{T}; job::Integer=jobs86.A) where T
+  if !m.flags.isanalysisdone || !m.flags.isfactordone
+    error("Cannot call solve! before factor!")
+    return x
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_solve(Cint(job), nrhs, ldx, x, m.order, m.__keep, m.control, m.info)
+  checkerror(m, "MA86: Error during solve.")
+  return x
+end
+solve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T = solve!(m, copy(b); kwargs...)
+
+function factorsolve!(m::Ma86{T}, x::VecOrMat{T}; 
+  matrix_type::Integer=matrix_types86.herm_indef
+) where T <: Ma86Data
+  if isnull(m.__keep) || !m.flags.isanalysisdone
+    error("Cannot call factorsolve! before analyse!")
+    return x 
+  end
+  size(x,1) != m.n && throw(DimensionMismatch("x must have $(m.n) rows, got $(size(x,1))."))
+  ldx = Cint(stride(x,2))
+  nrhs = Cint(size(x,2))
+  ma86_factor_solve(
+    Cint(matrix_type), m.n, m.colptr, m.rowval, m.nzval, m.order, m.__keep, m.control, m.info, 
+    nrhs, ldx, x, m.scale
+  )
+  checkerror(m, "MA86: Error during factorsolve.")
+  m.flags.isfactordone = true
+  return x
+end
+
+# Out-of-place version
+factorsolve(m::Ma86{T}, b::VecOrMat{T}; kwargs...) where T <: Ma86Data = 
+  factorsolve!(m, copy(b); kwargs...)
+
+# LinearAlgebra wrappers
+function Base.:\(m::Ma86{T}, b::VecOrMat{T}) where T <: Ma86Data 
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve(m, b)
+  factorsolve(m, b)
+end
+
+function LinearAlgebra.ldiv!(m::Ma86{T}, x::VecOrMat{T}) where T <: Ma86Data
+  !m.flags.isanalysisdone && analyse!(m)
+  m.flags.isfactordone && return solve!(m, x)
+  factorsolve!(m, x)
+end
+
+##############################
+# C Wrapper
+##############################
+
 for T in (Float32, Float64, ComplexF32, ComplexF64) 
   S = data_map[T]
 
@@ -95,7 +328,9 @@ for T in (Float32, Float64, ComplexF32, ComplexF64)
   end
 
   if T <: Ma86Real
-    @eval function ma86_factor(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
+    # Note that matrix_type argument is unused in these methods
+    # Included to make the signature match for all data types
+    @eval function ma86_factor(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, val::Vector{$T}, 
                               order::Vector{Cint}, keep::Keep, control::Ma86_Control{$S}, 
                               info::Ma86_Info{$S}, scale::Vector{$S})
       @assert length(ptr) == n+1
@@ -109,7 +344,7 @@ for T in (Float32, Float64, ComplexF32, ComplexF64)
       )
     end
 
-    @eval function ma86_factor_solve(n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
+    @eval function ma86_factor_solve(matrix_type::Cint, n::Cint, ptr::Vector{Cint}, row::Vector{Cint}, 
                                     val::Vector{$T}, order::Vector{Cint}, keep::Keep, 
                                     control::Ma86_Control{$S}, info::Ma86_Info{$S}, 
                                     nrhs::Cint, ldx::Cint, x::VecOrMat{$T}, 

From 20ec82901f65d98b7455c995b0db019584a870c8 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:47:58 -0500
Subject: [PATCH 12/16] Don't include deps in test

---
 test/runtests.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 1fd4aee..dcb23a2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,14 +4,14 @@ using LinearAlgebra
 using SparseArrays
 using Test
 
-include("../deps/deps.jl")
+# include("../deps/deps.jl")
 
-if @isdefined libhsl_ma57
+if isdefined(HSL, :libhsl_ma57)
   include("test_ma57.jl")
-  if @isdefined libhsl_ma57_patch
+  if isdefined(HSL, :libhsl_ma57_patch)
     include("test_ma57_patch.jl")
   end
 end
-if @isdefined libhsl_ma97
+if isdefined(HSL, :libhsl_ma97)
   include("test_ma97.jl")
 end

From 8b90b1bdb7e579015818d749e7c4a4dcc1c96907 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:48:07 -0500
Subject: [PATCH 13/16] Add test deps

---
 test/Project.toml | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 test/Project.toml

diff --git a/test/Project.toml b/test/Project.toml
new file mode 100644
index 0000000..59b801f
--- /dev/null
+++ b/test/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From d921e47cba1be0b1c745bd61a41d807513979156 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 17:48:37 -0500
Subject: [PATCH 14/16] Add ma86 tests

---
 test/runtests.jl  |   3 +
 test/test_ma86.jl | 178 ++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 181 insertions(+)
 create mode 100644 test/test_ma86.jl

diff --git a/test/runtests.jl b/test/runtests.jl
index dcb23a2..667edea 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,3 +15,6 @@ end
 if isdefined(HSL, :libhsl_ma97)
   include("test_ma97.jl")
 end
+if isdefined(HSL, :libhsl_ma86)
+  include("test_ma86.jl")
+end
diff --git a/test/test_ma86.jl b/test/test_ma86.jl
new file mode 100644
index 0000000..a4ac338
--- /dev/null
+++ b/test/test_ma86.jl
@@ -0,0 +1,178 @@
+using HSL
+using SparseArrays
+using LinearAlgebra
+using Test
+using Random
+
+function getallocs(A,b)
+    T = eltype(A)
+    n = size(A,1)
+    ma86 = Ma86(A)
+    allocs = @allocated HSL.analyse!(ma86)
+    allocs += @allocated HSL.factor!(ma86)
+    x = copy(b)
+    allocs += @allocated HSL.solve!(ma86, x)
+end
+
+T = Float64
+n = 10
+A = convert(T, 3) * sprand(T, n, n, 0.5)
+Aindef = A + A' + one(T)*I
+ma86 = Ma86(Aindef)
+ma86.flags.isanalysisdone
+HSL.factor!(ma86)
+b = rand(T,n)
+HSL.solve!(ma86, b)
+
+@testset "MA86 ($T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+    # @info "Testing hsl_ma86 with $T data"
+    n = 10
+    A = convert(T, 3) * sprand(T, n, n, 0.5)
+    Aindef = A + A' + one(T)*I
+
+    Apsd = A'A + one(T)*I
+    @test isposdef(Matrix(Apsd))
+
+
+    # Test symmetric indefinite A
+    @testset "$label" for (label,A) in (("Indefinite", Aindef), ("Positive Definite",Apsd))
+        n = size(A,1)
+        T = eltype(A)
+
+        # Test constructors
+        L = SparseMatrixCSC{T,Cint}(tril(A))
+        let ptr = L.colptr, row = L.rowval, val = L.nzval
+            # Check that new arrays aren't allocated if the input in Cint
+            let ma86 = Ma86(ptr, row, val)
+                @test ma86.colptr === ptr
+                @test ma86.rowval === row
+                @test ma86.nzval === val 
+                @test ma86.flags.isanalysisdone == true
+                @test !HSL.isnull(ma86.__keep)
+            end
+
+            # Check the keyword constructor            
+            control = Ma86_Control(T)
+            info = Ma86_Info(T)
+            S = HSL.data_map[T]
+            order = randperm(Cint(n))
+            scale = fill(one(S), n)
+            let ma86 = HSL.Ma86(ptr, row, val, 
+                    control=control, info=info, order=order, scale=scale)
+                @test ma86.order === order
+                @test ma86.scale === scale 
+                @test ma86.control === control
+                @test ma86.info === info 
+            end
+
+            # Pass in Int64 data
+            let ma86 = HSL.Ma86(Int64.(ptr), Int64.(row), val, 
+                    order=randperm(Int64(n)), scale=Float64.(scale))
+                @test ma86.colptr !== ptr
+                @test ma86.rowval !== row
+                @test ma86.order !== order
+                @test ma86.scale !== scale 
+            end
+
+            # Pass in lower trianglar matrix
+            # Make sure it doesn't reallocate new data
+            let ma86 = HSL.Ma86(L) 
+                @test ma86.colptr === L.colptr
+                @test ma86.rowval === L.rowval
+                @test ma86.nzval === L.nzval
+            end
+
+            # Check analyse skipping
+            let ma86 = HSL.Ma86(L, doanalyse=false)
+                @test ma86.flags.isanalysisdone == false
+            end
+        end
+
+        @testset "with $(size(b,2)) rhs vectors" for b in (rand(T, n), rand(T, n, 3))
+            n = size(A,1)
+            T = eltype(A)
+
+            # factor, then solve
+            ma86_0 = Ma86(A)
+            let ma86 = ma86_0
+                x = copy(b)
+                HSL.factor!(ma86)
+                @test ma86.info.flag == 0
+                HSL.solve!(ma86, x)
+                @test ma86.info.flag == 0
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+                @test HSL.solve(ma86, b) ≈ x  # out-of-place solve
+                @test ma86 \ b ≈ x            # backslash
+                x2 = copy(b)
+                ldiv!(ma86, x2)
+                @test x2 ≈ x
+            end
+
+            # in-place
+            ma86_1 = Ma86(A)
+            let ma86 = ma86_1
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # out-of-place
+            ma86_2 = Ma86(A)
+            let ma86 = ma86_2
+                x = HSL.factorsolve(ma86, b)
+                @test A*x ≈ b atol=sqrt(eps(norm(x)))
+            end
+
+            # Test linear algebra
+            ma86_3 = Ma86(A)
+            x = ma86_3 \ b
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            ma86_4 = Ma86(A)
+            x = copy(b)
+            ldiv!(ma86_4, x)
+            @test A*x ≈ b atol=sqrt(eps(norm(x)))
+
+            # Test errors
+            ma86_err = Ma86(Apsd, doanalyse=false)
+            let ma86 = ma86_err
+                @test_throws ErrorException HSL.factor!(ma86)
+                b = rand(T, n)
+                @test_throws ErrorException HSL.solve!(ma86, b)
+                @test_throws ErrorException HSL.factorsolve!(ma86, b)
+
+                HSL.analyse!(ma86)
+                @test_throws ErrorException HSL.solve!(ma86, b) 
+
+                b = rand(T, n)
+                x = copy(b)
+                HSL.factorsolve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                x = copy(b)
+                HSL.solve!(ma86, x)
+                @test Apsd*x ≈ b atol=sqrt(eps(norm(x)))
+
+                b_bad = rand(T, 2n)
+                @test_throws DimensionMismatch HSL.solve!(ma86, b_bad)
+                @test_throws DimensionMismatch HSL.factorsolve!(ma86, b_bad)
+
+                # Test allocations
+                @test getallocs(A, b) == 0
+            end
+            
+            # Test HSL Errors
+            # Prints out annoying warning message from MC78.
+            # Not sure how to avoid this.
+            # Should we even allow for non-square matrices to be passed to the constructor?
+            B = sprand(T, n, n, 0.1)
+            ma86_hslerr = Ma86(B, doanalyse=false)
+            ma86_hslerr.control.diagnostics_level = -1
+            ma86_hslerr.control.unit_diagnostics = -1
+            ma86_hslerr.control.unit_error = -1
+            ma86_hslerr.control.unit_warning = -1
+            @test_throws HSL.Ma86Exception HSL.analyse!(ma86_hslerr)
+        end
+    end
+end
+

From ba2c8d01fbb3cc4108e358eb944c34fac0e319cd Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 18:03:21 -0500
Subject: [PATCH 15/16] Modify tests a bit

---
 test/test_ma86.jl | 14 +++++++-------
 test/test_ma97.jl |  3 +++
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/test/test_ma86.jl b/test/test_ma86.jl
index a4ac338..52b8e97 100644
--- a/test/test_ma86.jl
+++ b/test/test_ma86.jl
@@ -165,13 +165,13 @@ HSL.solve!(ma86, b)
             # Prints out annoying warning message from MC78.
             # Not sure how to avoid this.
             # Should we even allow for non-square matrices to be passed to the constructor?
-            B = sprand(T, n, n, 0.1)
-            ma86_hslerr = Ma86(B, doanalyse=false)
-            ma86_hslerr.control.diagnostics_level = -1
-            ma86_hslerr.control.unit_diagnostics = -1
-            ma86_hslerr.control.unit_error = -1
-            ma86_hslerr.control.unit_warning = -1
-            @test_throws HSL.Ma86Exception HSL.analyse!(ma86_hslerr)
+            # B = sprand(T, n, n, 0.1)
+            # ma86_hslerr = Ma86(B, doanalyse=false)
+            # ma86_hslerr.control.diagnostics_level = -1
+            # ma86_hslerr.control.unit_diagnostics = -1
+            # ma86_hslerr.control.unit_error = -1
+            # ma86_hslerr.control.unit_warning = -1
+            # @test_throws HSL.Ma86Exception HSL.analyse!(ma86_hslerr)
         end
     end
 end
diff --git a/test/test_ma97.jl b/test/test_ma97.jl
index 1c9d5a7..17f7818 100644
--- a/test/test_ma97.jl
+++ b/test/test_ma97.jl
@@ -1,3 +1,6 @@
+using Random
+Random.seed!(1)
+
 function getallocs(A)
     T = eltype(A)
     n = size(A,1)

From 6cba88a0cd213d9fb0ff652a6605b3de83608374 Mon Sep 17 00:00:00 2001
From: Brian Jackson <bjack205@gmail.com>
Date: Tue, 25 Jan 2022 18:03:34 -0500
Subject: [PATCH 16/16] Update build

---
 deps/build.jl          | 17 +++++++++++------
 deps/build_hsl_ma86.jl | 19 ++++++++++---------
 2 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/deps/build.jl b/deps/build.jl
index 3f142e5..64e376c 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -10,7 +10,7 @@ struct HSLVersion
     ext::String
 end
 
-getname(ver::HSLVersion) = ver.algname * "-" * ver.version * ver.ext
+getname(ver::HSLVersion, ext=true) = ver.algname * "-" * ver.version * (ext ? ver.ext : "")
 
 function checksha(version::HSLVersion, filepath)
     if isfile(filepath)
@@ -65,6 +65,16 @@ const hsl_ma97_path = haskey(ENV, "HSL_MA97_PATH") ? ENV["HSL_MA97_PATH"] : join
 hsl_ma97_version = findversion(hsl_ma97_versions, hsl_ma97_path)
 const hsl_ma97_archive = isnothing(hsl_ma97_version) ? "" : joinpath(hsl_ma97_path, getname(hsl_ma97_version))
 
+##############################
+# MA86
+##############################
+const hsl_ma86_versions = [
+  HSLVersion("hsl_ma86", "1.7.0",  "845c65bee0bf31507d7b99c87773997012b7b16434da349ad5c361ceb257191e", ".tar.gz"),
+]
+const hsl_ma86_path = haskey(ENV, "HSL_MA86_PATH") ? ENV["HSL_MA86_PATH"] : joinpath(@__DIR__, "downloads")
+hsl_ma86_version = findversion(hsl_ma86_versions, hsl_ma86_path)
+const hsl_ma86_archive = isnothing(hsl_ma86_version) ? "" : joinpath(hsl_ma86_path, getname(hsl_ma86_version))
+
 ##############################
 # Build
 ##############################
@@ -74,11 +84,6 @@ const HSL_FC = haskey(ENV, "HSL_FC") ? ENV["HSL_FC"] : "gfortran"
 const HSL_F77 = haskey(ENV, "HSL_F77") ? ENV["HSL_F77"] : HSL_FC
 const HSL_CC = haskey(ENV, "HSL_CC") ? ENV["HSL_CC"] : "gcc"
 
-const hsl_ma86_version = "1.7.0"
-const hsl_ma86_sha256 = "845c65bee0bf31507d7b99c87773997012b7b16434da349ad5c361ceb257191e"
-const hsl_ma86_path = haskey(ENV, "HSL_MA86_PATH") ? ENV["HSL_MA86_PATH"] : joinpath(@__DIR__, "downloads")
-const hsl_ma86_archive = joinpath(hsl_ma86_path, "hsl_ma86-$hsl_ma86_version.tar.gz")
-
 const so         = Sys.isapple() ? "dylib" : "so"
 const all_load   = Sys.isapple() ? "-all_load" : "--whole-archive"
 const noall_load = Sys.isapple() ? "-noall_load" : "--no-whole-archive"
diff --git a/deps/build_hsl_ma86.jl b/deps/build_hsl_ma86.jl
index 0ad7224..4b31b9e 100644
--- a/deps/build_hsl_ma86.jl
+++ b/deps/build_hsl_ma86.jl
@@ -1,11 +1,12 @@
-open(hsl_ma86_archive) do f
-  if bytes2hex(sha256(f)) != hsl_ma86_sha256
-    error("SHA256 of HSL MA86 doesn't match")
-  end
+version = hsl_ma86_version
+archive = hsl_ma86_archive
+if version.ext == ".tar.gz"
+    run(`tar -zxf $archive -C $builddir`)
+elseif hsl_ma97_version.ext == ".zip"
+    run(`unzip -o $archive -d $builddir`)
 end
-run(`tar -zxf $hsl_ma86_archive -C $builddir`)
-cd("$builddir/hsl_ma86-$hsl_ma86_version")
-run(`./configure --prefix=$usrdir F77=gfortran CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --with-blas="-L$libopenblas_dir -lopenblas" --with-metis="-L$libmetis_dir -lmetis"`)
+name = getname(version, false)
+cd("$builddir/$name")
+run(`./configure --prefix=$usrdir FC=$(HSL_FC) F77=$(HSL_F77) CC=$(HSL_CC) CFLAGS=-fPIC FFLAGS="-fPIC -fopenmp" FCFLAGS="-fPIC -fopenmp" --with-blas="-L$libopenblas_dir -lopenblas" --with-metis="-L$libmetis_dir -lmetis"`)
 run(`make install`)
-run(`gfortran -fPIC -shared -Wl,$all_load $libdir/libhsl_ma86.a -L$libopenblas_dir -lopenblas -L$libmetis_dir -lmetis -lgomp -Wl,$noall_load -o $libdir/libhsl_ma86.$so`)
-cd(@__DIR__)
\ No newline at end of file
+cd(@__DIR__)