From 9f7129df1489460571d9e7fbdd7cfd029eb07bc1 Mon Sep 17 00:00:00 2001
From: Jash <jashambaliya1@gmail.com>
Date: Mon, 27 Apr 2026 10:05:57 +0530
Subject: [PATCH 1/5] Add HSL wrapper

---
 .github/workflows/HSL.yml |  18 +++++
 Project.toml              |   3 +
 ext/LinearSolveHSLExt.jl  | 164 ++++++++++++++++++++++++++++++++++++++
 src/LinearSolve.jl        |   4 +-
 src/extension_algs.jl     |  72 +++++++++++++++++
 test/hsl/Project.toml     |   7 ++
 test/hsl/hsl.jl           |  63 +++++++++++++++
 test/runtests.jl          |   5 ++
 8 files changed, 335 insertions(+), 1 deletion(-)
 create mode 100644 .github/workflows/HSL.yml
 create mode 100644 ext/LinearSolveHSLExt.jl
 create mode 100644 test/hsl/Project.toml
 create mode 100644 test/hsl/hsl.jl

diff --git a/.github/workflows/HSL.yml b/.github/workflows/HSL.yml
new file mode 100644
index 000000000..08c500e94
--- /dev/null
+++ b/.github/workflows/HSL.yml
@@ -0,0 +1,18 @@
+name: "HSL Tests"
+
+on:
+  workflow_dispatch:
+
+jobs:
+  hsl-tests:
+    name: "LinearSolveHSL"
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v3
+        with:
+          version: "1"
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          GROUP: "LinearSolveHSL"
diff --git a/Project.toml b/Project.toml
index e1a28ee3d..8f510fcd8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -46,6 +46,7 @@ FastLapackInterface = "29a986be-02c6-4525-aec4-84b980013641"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ginkgo = "4c8bd3c9-ead9-4b5e-a625-08f1338ba0ec"
 HYPRE = "b5ffcf37-a2bd-41ab-a3da-4bd9bc8ad771"
+HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
@@ -84,6 +85,7 @@ LinearSolveFastLapackInterfaceExt = "FastLapackInterface"
 LinearSolveForwardDiffExt = "ForwardDiff"
 LinearSolveGinkgoExt = ["Ginkgo", "SparseArrays"]
 LinearSolveHYPREExt = "HYPRE"
+LinearSolveHSLExt = ["HSL", "SparseArrays"]
 LinearSolveIterativeSolversExt = "IterativeSolvers"
 LinearSolveKernelAbstractionsExt = "KernelAbstractions"
 LinearSolveKrylovKitExt = "KrylovKit"
@@ -126,6 +128,7 @@ ForwardDiff = "0.10.38, 1"
 GPUArraysCore = "0.2"
 Ginkgo = "1"
 HYPRE = "1.7"
+HSL = "0.5"
 InteractiveUtils = "1.10"
 IterativeSolvers = "0.9.4"
 KernelAbstractions = "0.9.30"
diff --git a/ext/LinearSolveHSLExt.jl b/ext/LinearSolveHSLExt.jl
new file mode 100644
index 000000000..abc39c2df
--- /dev/null
+++ b/ext/LinearSolveHSLExt.jl
@@ -0,0 +1,164 @@
+module LinearSolveHSLExt
+
+using HSL
+using LinearSolve
+using LinearSolve: LinearVerbosity, OperatorAssumptions
+using SciMLBase: SciMLBase, ReturnCode
+using SciMLLogging: @SciMLMessage
+using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC
+using LinearAlgebra: ishermitian, issymmetric
+
+mutable struct HSLMA57Cache{M}
+    ma57::M
+end
+
+mutable struct HSLMA97Cache{M}
+    ma97::M
+end
+
+function _is_symmetric_like(A)
+    T = eltype(A)
+    if T <: Real
+        return issymmetric(A)
+    else
+        return ishermitian(A)
+    end
+end
+
+function _as_int_csc(A::SparseMatrixCSC{T, Int}) where {T}
+    return A
+end
+
+function _as_int_csc(A::AbstractSparseMatrixCSC{T}) where {T}
+    return SparseMatrixCSC{T, Int}(A)
+end
+
+function LinearSolve.init_cacheval(
+        alg::LinearSolve.HSLMA57Factorization,
+        A::AbstractSparseMatrixCSC{<:Union{Float32, Float64}}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return HSLMA57Cache(HSL.Ma57(A; alg.kwargs...))
+end
+
+function LinearSolve.init_cacheval(
+        alg::LinearSolve.HSLMA57Factorization,
+        A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return nothing
+end
+
+function LinearSolve.init_cacheval(
+        alg::LinearSolve.HSLMA97Factorization,
+        A::AbstractSparseMatrixCSC{<:Union{Float32, Float64, ComplexF32, ComplexF64}}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    A_int = _as_int_csc(A)
+    return HSLMA97Cache(HSL.Ma97(A_int; alg.kwargs...))
+end
+
+function LinearSolve.init_cacheval(
+        alg::LinearSolve.HSLMA97Factorization,
+        A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    return nothing
+end
+
+function SciMLBase.solve!(
+        cache::LinearSolve.LinearCache,
+        alg::LinearSolve.HSLMA57Factorization;
+        kwargs...
+    )
+    A = convert(AbstractMatrix, cache.A)
+    A isa AbstractSparseMatrixCSC ||
+        error("HSLMA57Factorization currently supports only sparse CSC matrices")
+    size(A, 1) == size(A, 2) ||
+        error("HSLMA57Factorization requires a square matrix")
+    _is_symmetric_like(A) ||
+        error("HSLMA57Factorization requires a symmetric/Hermitian matrix")
+
+    hcache = LinearSolve.@get_cacheval(cache, :HSLMA57Factorization)
+    hcache === nothing && error(
+        "HSLMA57Factorization supports `AbstractSparseMatrixCSC{<:Union{Float32, Float64}}`"
+    )
+
+    try
+        if cache.isfresh
+            hcache.ma57 = HSL.Ma57(A; alg.kwargs...)
+            HSL.ma57_factorize!(hcache.ma57)
+            cache.isfresh = false
+        end
+
+        cache.u .= HSL.ma57_solve(hcache.ma57, cache.b)
+
+        return SciMLBase.build_linear_solution(
+            alg, cache.u, nothing, cache; retcode = ReturnCode.Success
+        )
+    catch err
+        if err isa HSL.Ma57Exception
+            @SciMLMessage(
+                "MA57 failed: $(err.msg) (flag $(err.flag))",
+                cache.verbose,
+                :solver_failure
+            )
+            return SciMLBase.build_linear_solution(
+                alg, cache.u, nothing, cache; retcode = ReturnCode.Failure
+            )
+        end
+        rethrow(err)
+    end
+end
+
+function SciMLBase.solve!(
+        cache::LinearSolve.LinearCache,
+        alg::LinearSolve.HSLMA97Factorization;
+        kwargs...
+    )
+    A = convert(AbstractMatrix, cache.A)
+    A isa AbstractSparseMatrixCSC ||
+        error("HSLMA97Factorization currently supports only sparse CSC matrices")
+    size(A, 1) == size(A, 2) ||
+        error("HSLMA97Factorization requires a square matrix")
+    _is_symmetric_like(A) ||
+        error("HSLMA97Factorization requires a symmetric/Hermitian matrix")
+
+    hcache = LinearSolve.@get_cacheval(cache, :HSLMA97Factorization)
+    hcache === nothing && error(
+        "HSLMA97Factorization supports `AbstractSparseMatrixCSC{<:Union{Float32, Float64, ComplexF32, ComplexF64}}`"
+    )
+
+    try
+        if cache.isfresh
+            A_int = _as_int_csc(A)
+            hcache.ma97 = HSL.Ma97(A_int; alg.kwargs...)
+            HSL.ma97_factorize!(hcache.ma97, matrix_type = alg.matrix_type)
+            cache.isfresh = false
+        end
+
+        cache.u .= HSL.ma97_solve(hcache.ma97, cache.b)
+
+        return SciMLBase.build_linear_solution(
+            alg, cache.u, nothing, cache; retcode = ReturnCode.Success
+        )
+    catch err
+        if err isa HSL.Ma97Exception
+            @SciMLMessage(
+                "MA97 failed: $(err.msg) (flag $(err.flag))",
+                cache.verbose,
+                :solver_failure
+            )
+            return SciMLBase.build_linear_solution(
+                alg, cache.u, nothing, cache; retcode = ReturnCode.Failure
+            )
+        end
+        rethrow(err)
+    end
+end
+
+end
diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl
index ae40d1eee..00a2d6268 100644
--- a/src/LinearSolve.jl
+++ b/src/LinearSolve.jl
@@ -483,6 +483,7 @@ for alg in (
         :DiagonalFactorization, :CholeskyFactorization, :BunchKaufmanFactorization,
         :CHOLMODFactorization, :LDLtFactorization, :AppleAccelerateLUFactorization,
         :MKLLUFactorization, :MetalLUFactorization, :CUSOLVERRFFactorization, :ParUFactorization,
+        :HSLMA57Factorization, :HSLMA97Factorization,
         :STRUMPACKFactorization,
     )
     @eval needs_square_A(::$(alg)) = true
@@ -526,7 +527,8 @@ export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
     BunchKaufmanFactorization, CHOLMODFactorization, LDLtFactorization,
     CUSOLVERRFFactorization, CliqueTreesFactorization, ParUFactorization,
     STRUMPACKFactorization, MUMPSFactorization,
-    SpecializedLUFactorization, SpecializedQRFactorization
+    SpecializedLUFactorization, SpecializedQRFactorization,
+    HSLMA57Factorization, HSLMA97Factorization
 
 export LinearSolveFunction, DirectLdiv!, show_algorithm_choices
 
diff --git a/src/extension_algs.jl b/src/extension_algs.jl
index a064372aa..d080da6b7 100644
--- a/src/extension_algs.jl
+++ b/src/extension_algs.jl
@@ -1461,6 +1461,78 @@ struct MUMPSFactorization <: AbstractSparseFactorization
     end
 end
 
+"""
+`HSLMA57Factorization(; kwargs...)`
+
+A sparse symmetric direct solver powered by
+[HSL.jl](https://github.com/JuliaSmoothOptimizers/HSL.jl)'s MA57 backend.
+
+Keyword arguments are forwarded to `HSL.Ma57(...)`.
+
+!!! note
+
+    Using this solver requires loading `HSL.jl` and `SparseArrays`:
+    ```julia
+    using HSL, SparseArrays
+    ```
+
+    `HSL.jl` requires a manual installation of `HSL_jll.jl` due proprietary
+    licensing. See `HSL.jl` installation instructions:
+    https://github.com/JuliaSmoothOptimizers/HSL.jl
+"""
+struct HSLMA57Factorization{K} <: AbstractSparseFactorization
+    kwargs::K
+
+    function HSLMA57Factorization(; kwargs...)
+        ext = Base.get_extension(@__MODULE__, :LinearSolveHSLExt)
+        if ext === nothing
+            error("HSLMA57Factorization requires `using HSL, SparseArrays`. Note: `HSL.jl` requires manual installation of `HSL_jll.jl`.")
+        end
+        return new{typeof(kwargs)}(kwargs)
+    end
+end
+
+"""
+`HSLMA97Factorization(; matrix_type = :real_indef, kwargs...)`
+
+A sparse symmetric/Hermitian direct solver powered by
+[HSL.jl](https://github.com/JuliaSmoothOptimizers/HSL.jl)'s MA97 backend.
+
+## Keyword Arguments
+
+  - `matrix_type`: Passed to `HSL.ma97_factorize!`. Supported values include
+    `:real_spd`, `:real_indef`, `:herm_pd`, `:herm_indef`, `:cmpl_indef`.
+  - `kwargs...`: Forwarded to `HSL.Ma97(...)`.
+
+!!! note
+
+    Using this solver requires loading `HSL.jl` and `SparseArrays`:
+    ```julia
+    using HSL, SparseArrays
+    ```
+
+    `HSL.jl` requires a manual installation of `HSL_jll.jl` due proprietary
+    licensing. See `HSL.jl` installation instructions:
+    https://github.com/JuliaSmoothOptimizers/HSL.jl
+"""
+struct HSLMA97Factorization{K} <: AbstractSparseFactorization
+    matrix_type::Symbol
+    kwargs::K
+
+    function HSLMA97Factorization(; matrix_type::Symbol = :real_indef, kwargs...)
+        ext = Base.get_extension(@__MODULE__, :LinearSolveHSLExt)
+        if ext === nothing
+            error("HSLMA97Factorization requires `using HSL, SparseArrays`. Note: `HSL.jl` requires manual installation of `HSL_jll.jl`.")
+        end
+
+        matrix_type ∈ (:real_spd, :real_indef, :herm_pd, :herm_indef, :cmpl_indef) || error(
+            "Unsupported matrix_type: $(matrix_type). Expected one of :real_spd, :real_indef, :herm_pd, :herm_indef, :cmpl_indef."
+        )
+
+        return new{typeof(kwargs)}(matrix_type, kwargs)
+    end
+end
+
 """
     ElementalJL(; method = :LU)
 
diff --git a/test/hsl/Project.toml b/test/hsl/Project.toml
new file mode 100644
index 000000000..07a8c90fa
--- /dev/null
+++ b/test/hsl/Project.toml
@@ -0,0 +1,7 @@
+[deps]
+HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/hsl/hsl.jl b/test/hsl/hsl.jl
new file mode 100644
index 000000000..86b94bee0
--- /dev/null
+++ b/test/hsl/hsl.jl
@@ -0,0 +1,63 @@
+using HSL
+using LinearSolve
+using LinearAlgebra
+using Random
+using SparseArrays
+using Test
+
+Random.seed!(1234)
+const REQUIRE_FUNCTIONAL_HSL =
+    lowercase(get(ENV, "LINEARSOLVE_HSL_REQUIRE_FUNCTIONAL", "false")) == "true"
+
+if !HSL.LIBHSL_isfunctional()
+    if REQUIRE_FUNCTIONAL_HSL
+        error(
+            "LINEARSOLVE_HSL_REQUIRE_FUNCTIONAL=true but HSL_jll is not functional in this environment."
+        )
+    else
+        @info "Skipping LinearSolveHSL tests: HSL_jll is not functional in this environment"
+        @test true
+    end
+else
+    @testset "HSL MA57 wrapper" begin
+        n = 40
+        A = sprand(Float64, n, n, 0.12)
+        A = A + A' + 40.0I
+        b = rand(Float64, n)
+        prob = LinearProblem(A, b)
+
+        sol = solve(prob, HSLMA57Factorization())
+        @test sol.retcode == ReturnCode.Success
+        @test A * sol.u ≈ b rtol = 1.0e-9 atol = 1.0e-11
+
+        cache = init(prob, HSLMA57Factorization())
+        sol1 = solve!(cache)
+        cache.b = rand(Float64, n)
+        b2 = copy(cache.b)
+        sol2 = solve!(cache)
+        @test sol1.retcode == ReturnCode.Success
+        @test sol2.retcode == ReturnCode.Success
+        @test A * sol2.u ≈ b2 rtol = 1.0e-9 atol = 1.0e-11
+    end
+
+    @testset "HSL MA97 wrapper" begin
+        n = 35
+        A = sprand(Float64, n, n, 0.1)
+        A = A + A' + 30.0I
+        b = rand(Float64, n)
+        prob = LinearProblem(A, b)
+
+        sol = solve(prob, HSLMA97Factorization(matrix_type = :real_spd))
+        @test sol.retcode == ReturnCode.Success
+        @test A * sol.u ≈ b rtol = 1.0e-9 atol = 1.0e-11
+
+        cache = init(prob, HSLMA97Factorization(matrix_type = :real_spd))
+        sol1 = solve!(cache)
+        cache.b = rand(Float64, n)
+        b2 = copy(cache.b)
+        sol2 = solve!(cache)
+        @test sol1.retcode == ReturnCode.Success
+        @test sol2.retcode == ReturnCode.Success
+        @test A * sol2.u ≈ b2 rtol = 1.0e-9 atol = 1.0e-11
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 8b91949db..88c8e0af7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -161,6 +161,11 @@ else
         @time @safetestset "Pardiso" include("pardiso/pardiso.jl")
     end
 
+    if GROUP == "LinearSolveHSL"
+        activate_group_env("hsl")
+        @time @safetestset "HSL" include("hsl/hsl.jl")
+    end
+
     if Base.Sys.islinux() && GROUP == "LinearSolveMUMPS"
         activate_group_env("mumps")
         @time @safetestset "MUMPS" include("mumps/mumps.jl")

From 766903e1c37c83cbe3346b37ef552c65aa3fc833 Mon Sep 17 00:00:00 2001
From: Jash <jashambaliya1@gmail.com>
Date: Mon, 27 Apr 2026 10:46:31 +0530
Subject: [PATCH 2/5] added guards to HSL in resolve

---
 test/core/resolve.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/test/core/resolve.jl b/test/core/resolve.jl
index 4079cdac8..2b2e54b4b 100644
--- a/test/core/resolve.jl
+++ b/test/core/resolve.jl
@@ -64,6 +64,14 @@ for alg in vcat(
         ) &&
             (
             !(alg == MUMPSFactorization) || HAS_MUMPS
+        ) &&
+            (
+            !(alg == HSLMA57Factorization) ||
+                Base.get_extension(LinearSolve, :LinearSolveHSLExt) !== nothing
+        ) &&
+            (
+            !(alg == HSLMA97Factorization) ||
+                Base.get_extension(LinearSolve, :LinearSolveHSLExt) !== nothing
         )
         A = [1.0 2.0; 3.0 4.0]
         alg in [

From b674d2e5d2bd1fac30a6506fb66715fb02d184cf Mon Sep 17 00:00:00 2001
From: Jash <jashambaliya1@gmail.com>
Date: Tue, 5 May 2026 20:31:06 +0530
Subject: [PATCH 3/5] patched MA57 and MA97  to accept Symmetric/Hermitian
 sparse wrappers and to solve in-place through the cache instead of allocating
 via HSL.ma57_solve / HSL.ma97_solve.

---
 ext/LinearSolveHSLExt.jl | 80 +++++++++++++++++++++++++---------------
 test/hsl/hsl.jl          | 11 ++++++
 2 files changed, 62 insertions(+), 29 deletions(-)

diff --git a/ext/LinearSolveHSLExt.jl b/ext/LinearSolveHSLExt.jl
index abc39c2df..03b3a7e40 100644
--- a/ext/LinearSolveHSLExt.jl
+++ b/ext/LinearSolveHSLExt.jl
@@ -6,10 +6,11 @@ using LinearSolve: LinearVerbosity, OperatorAssumptions
 using SciMLBase: SciMLBase, ReturnCode
 using SciMLLogging: @SciMLMessage
 using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC
-using LinearAlgebra: ishermitian, issymmetric
+using LinearAlgebra: Hermitian, Symmetric, ishermitian, issymmetric, parent
 
-mutable struct HSLMA57Cache{M}
+mutable struct HSLMA57Cache{M, W}
     ma57::M
+    work::W
 end
 
 mutable struct HSLMA97Cache{M}
@@ -25,6 +26,17 @@ function _is_symmetric_like(A)
     end
 end
 
+function _sparse_csc_matrix(A::AbstractSparseMatrixCSC)
+    return A
+end
+
+function _sparse_csc_matrix(A::Union{
+        Symmetric{T, <:AbstractSparseMatrixCSC{T}},
+        Hermitian{T, <:AbstractSparseMatrixCSC{T}},
+    }) where {T}
+    return parent(A)
+end
+
 function _as_int_csc(A::SparseMatrixCSC{T, Int}) where {T}
     return A
 end
@@ -33,13 +45,13 @@ function _as_int_csc(A::AbstractSparseMatrixCSC{T}) where {T}
     return SparseMatrixCSC{T, Int}(A)
 end
 
-function LinearSolve.init_cacheval(
-        alg::LinearSolve.HSLMA57Factorization,
-        A::AbstractSparseMatrixCSC{<:Union{Float32, Float64}}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol,
-        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
-    )
-    return HSLMA57Cache(HSL.Ma57(A; alg.kwargs...))
+_hsl_rhs_ncols(b::AbstractVector) = 1
+_hsl_rhs_ncols(b::AbstractMatrix) = size(b, 2)
+
+function _resize_ma57_work!(hcache::HSLMA57Cache, b)
+    lwork = hcache.ma57.n * _hsl_rhs_ncols(b)
+    length(hcache.work) == lwork || resize!(hcache.work, lwork)
+    return hcache.work
 end
 
 function LinearSolve.init_cacheval(
@@ -48,17 +60,13 @@ function LinearSolve.init_cacheval(
         maxiters::Int, abstol, reltol,
         verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
     )
-    return nothing
-end
-
-function LinearSolve.init_cacheval(
-        alg::LinearSolve.HSLMA97Factorization,
-        A::AbstractSparseMatrixCSC{<:Union{Float32, Float64, ComplexF32, ComplexF64}}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol,
-        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
-    )
-    A_int = _as_int_csc(A)
-    return HSLMA97Cache(HSL.Ma97(A_int; alg.kwargs...))
+    A_sparse = _sparse_csc_matrix(A)
+    if !(A_sparse isa AbstractSparseMatrixCSC{<:Union{Float32, Float64}})
+        return nothing
+    end
+    ma57 = HSL.Ma57(A_sparse; alg.kwargs...)
+    work = Vector{eltype(A_sparse)}(undef, ma57.n * _hsl_rhs_ncols(b))
+    return HSLMA57Cache(ma57, work)
 end
 
 function LinearSolve.init_cacheval(
@@ -67,7 +75,14 @@ function LinearSolve.init_cacheval(
         maxiters::Int, abstol, reltol,
         verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
     )
-    return nothing
+    A_sparse = _sparse_csc_matrix(A)
+    if !(A_sparse isa AbstractSparseMatrixCSC{
+            <:Union{Float32, Float64, ComplexF32, ComplexF64},
+        })
+        return nothing
+    end
+    A_int = _as_int_csc(A_sparse)
+    return HSLMA97Cache(HSL.Ma97(A_int; alg.kwargs...))
 end
 
 function SciMLBase.solve!(
@@ -76,9 +91,10 @@ function SciMLBase.solve!(
         kwargs...
     )
     A = convert(AbstractMatrix, cache.A)
-    A isa AbstractSparseMatrixCSC ||
+    A_sparse = _sparse_csc_matrix(A)
+    A_sparse isa AbstractSparseMatrixCSC ||
         error("HSLMA57Factorization currently supports only sparse CSC matrices")
-    size(A, 1) == size(A, 2) ||
+    size(A_sparse, 1) == size(A_sparse, 2) ||
         error("HSLMA57Factorization requires a square matrix")
     _is_symmetric_like(A) ||
         error("HSLMA57Factorization requires a symmetric/Hermitian matrix")
@@ -90,12 +106,14 @@ function SciMLBase.solve!(
 
     try
         if cache.isfresh
-            hcache.ma57 = HSL.Ma57(A; alg.kwargs...)
+            hcache.ma57 = HSL.Ma57(A_sparse; alg.kwargs...)
             HSL.ma57_factorize!(hcache.ma57)
             cache.isfresh = false
         end
 
-        cache.u .= HSL.ma57_solve(hcache.ma57, cache.b)
+        _resize_ma57_work!(hcache, cache.b)
+        copyto!(cache.u, cache.b)
+        HSL.ma57_solve!(hcache.ma57, cache.u, hcache.work)
 
         return SciMLBase.build_linear_solution(
             alg, cache.u, nothing, cache; retcode = ReturnCode.Success
@@ -107,6 +125,7 @@ function SciMLBase.solve!(
                 cache.verbose,
                 :solver_failure
             )
+            cache.isfresh = true
             return SciMLBase.build_linear_solution(
                 alg, cache.u, nothing, cache; retcode = ReturnCode.Failure
             )
@@ -121,9 +140,10 @@ function SciMLBase.solve!(
         kwargs...
     )
     A = convert(AbstractMatrix, cache.A)
-    A isa AbstractSparseMatrixCSC ||
+    A_sparse = _sparse_csc_matrix(A)
+    A_sparse isa AbstractSparseMatrixCSC ||
         error("HSLMA97Factorization currently supports only sparse CSC matrices")
-    size(A, 1) == size(A, 2) ||
+    size(A_sparse, 1) == size(A_sparse, 2) ||
         error("HSLMA97Factorization requires a square matrix")
     _is_symmetric_like(A) ||
         error("HSLMA97Factorization requires a symmetric/Hermitian matrix")
@@ -135,13 +155,14 @@ function SciMLBase.solve!(
 
     try
         if cache.isfresh
-            A_int = _as_int_csc(A)
+            A_int = _as_int_csc(A_sparse)
             hcache.ma97 = HSL.Ma97(A_int; alg.kwargs...)
             HSL.ma97_factorize!(hcache.ma97, matrix_type = alg.matrix_type)
             cache.isfresh = false
         end
 
-        cache.u .= HSL.ma97_solve(hcache.ma97, cache.b)
+        copyto!(cache.u, cache.b)
+        HSL.ma97_solve!(hcache.ma97, cache.u)
 
         return SciMLBase.build_linear_solution(
             alg, cache.u, nothing, cache; retcode = ReturnCode.Success
@@ -153,6 +174,7 @@ function SciMLBase.solve!(
                 cache.verbose,
                 :solver_failure
             )
+            cache.isfresh = true
             return SciMLBase.build_linear_solution(
                 alg, cache.u, nothing, cache; retcode = ReturnCode.Failure
             )
diff --git a/test/hsl/hsl.jl b/test/hsl/hsl.jl
index 86b94bee0..b954d36f8 100644
--- a/test/hsl/hsl.jl
+++ b/test/hsl/hsl.jl
@@ -25,10 +25,12 @@ else
         A = A + A' + 40.0I
         b = rand(Float64, n)
         prob = LinearProblem(A, b)
+        prob_sym = LinearProblem(Symmetric(A), b)
 
         sol = solve(prob, HSLMA57Factorization())
         @test sol.retcode == ReturnCode.Success
         @test A * sol.u ≈ b rtol = 1.0e-9 atol = 1.0e-11
+        @test solve(prob_sym, HSLMA57Factorization()).u ≈ sol.u
 
         cache = init(prob, HSLMA57Factorization())
         sol1 = solve!(cache)
@@ -46,10 +48,12 @@ else
         A = A + A' + 30.0I
         b = rand(Float64, n)
         prob = LinearProblem(A, b)
+        prob_sym = LinearProblem(Symmetric(A), b)
 
         sol = solve(prob, HSLMA97Factorization(matrix_type = :real_spd))
         @test sol.retcode == ReturnCode.Success
         @test A * sol.u ≈ b rtol = 1.0e-9 atol = 1.0e-11
+        @test solve(prob_sym, HSLMA97Factorization(matrix_type = :real_spd)).u ≈ sol.u
 
         cache = init(prob, HSLMA97Factorization(matrix_type = :real_spd))
         sol1 = solve!(cache)
@@ -59,5 +63,12 @@ else
         @test sol1.retcode == ReturnCode.Success
         @test sol2.retcode == ReturnCode.Success
         @test A * sol2.u ≈ b2 rtol = 1.0e-9 atol = 1.0e-11
+
+        Ac = ComplexF64.(A)
+        bc = rand(ComplexF64, n)
+        prob_herm = LinearProblem(Hermitian(Ac), bc)
+        solc = solve(prob_herm, HSLMA97Factorization(matrix_type = :herm_pd))
+        @test solc.retcode == ReturnCode.Success
+        @test Hermitian(Ac) * solc.u ≈ bc rtol = 1.0e-9 atol = 1.0e-11
     end
 end

From cbd0f07f74216f3caa0205a3998167e1f88051ac Mon Sep 17 00:00:00 2001
From: Jash <jashambaliya1@gmail.com>
Date: Tue, 5 May 2026 20:34:44 +0530
Subject: [PATCH 4/5] ran runic

---
 ext/LinearSolveHSLExt.jl | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/ext/LinearSolveHSLExt.jl b/ext/LinearSolveHSLExt.jl
index 03b3a7e40..7779acba4 100644
--- a/ext/LinearSolveHSLExt.jl
+++ b/ext/LinearSolveHSLExt.jl
@@ -30,10 +30,12 @@ function _sparse_csc_matrix(A::AbstractSparseMatrixCSC)
     return A
 end
 
-function _sparse_csc_matrix(A::Union{
-        Symmetric{T, <:AbstractSparseMatrixCSC{T}},
-        Hermitian{T, <:AbstractSparseMatrixCSC{T}},
-    }) where {T}
+function _sparse_csc_matrix(
+        A::Union{
+            Symmetric{T, <:AbstractSparseMatrixCSC{T}},
+            Hermitian{T, <:AbstractSparseMatrixCSC{T}},
+        }
+    ) where {T}
     return parent(A)
 end
 
@@ -76,9 +78,11 @@ function LinearSolve.init_cacheval(
         verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
     )
     A_sparse = _sparse_csc_matrix(A)
-    if !(A_sparse isa AbstractSparseMatrixCSC{
-            <:Union{Float32, Float64, ComplexF32, ComplexF64},
-        })
+    if !(
+            A_sparse isa AbstractSparseMatrixCSC{
+                <:Union{Float32, Float64, ComplexF32, ComplexF64},
+            }
+        )
         return nothing
     end
     A_int = _as_int_csc(A_sparse)

From 07f8251fb40d493a193af0d56dbe39fb93fa2e24 Mon Sep 17 00:00:00 2001
From: Jash <jashambaliya1@gmail.com>
Date: Thu, 11 Jun 2026 14:19:49 +0530
Subject: [PATCH 5/5] use the new test_groups in test

---
 .github/workflows/HSL.yml | 18 ------------------
 test/test_groups.toml     |  5 +++++
 2 files changed, 5 insertions(+), 18 deletions(-)
 delete mode 100644 .github/workflows/HSL.yml

diff --git a/.github/workflows/HSL.yml b/.github/workflows/HSL.yml
deleted file mode 100644
index 08c500e94..000000000
--- a/.github/workflows/HSL.yml
+++ /dev/null
@@ -1,18 +0,0 @@
-name: "HSL Tests"
-
-on:
-  workflow_dispatch:
-
-jobs:
-  hsl-tests:
-    name: "LinearSolveHSL"
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v6
-      - uses: julia-actions/setup-julia@v3
-        with:
-          version: "1"
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
-        env:
-          GROUP: "LinearSolveHSL"
diff --git a/test/test_groups.toml b/test/test_groups.toml
index fd0457df4..abe61bd46 100644
--- a/test/test_groups.toml
+++ b/test/test_groups.toml
@@ -73,6 +73,11 @@ versions = ["lts", "1"]
 runner = "self-hosted"
 num_threads = 2
 
+[LinearSolveHSL]
+versions = ["lts", "1"]
+runner = "self-hosted"
+num_threads = 2
+
 [LinearSolveGinkgo]
 versions = ["lts", "1"]
 runner = "self-hosted"