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..7779acba4
--- /dev/null
+++ b/ext/LinearSolveHSLExt.jl
@@ -0,0 +1,190 @@
+module LinearSolveHSLExt
+
+using HSL
+using LinearSolve
+using LinearSolve: LinearVerbosity, OperatorAssumptions
+using SciMLBase: SciMLBase, ReturnCode
+using SciMLLogging: @SciMLMessage
+using SparseArrays: AbstractSparseMatrixCSC, SparseMatrixCSC
+using LinearAlgebra: Hermitian, Symmetric, ishermitian, issymmetric, parent
+
+mutable struct HSLMA57Cache{M, W}
+    ma57::M
+    work::W
+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 _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
+
+function _as_int_csc(A::AbstractSparseMatrixCSC{T}) where {T}
+    return SparseMatrixCSC{T, Int}(A)
+end
+
+_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(
+        alg::LinearSolve.HSLMA57Factorization,
+        A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    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(
+        alg::LinearSolve.HSLMA97Factorization,
+        A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions
+    )
+    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!(
+        cache::LinearSolve.LinearCache,
+        alg::LinearSolve.HSLMA57Factorization;
+        kwargs...
+    )
+    A = convert(AbstractMatrix, cache.A)
+    A_sparse = _sparse_csc_matrix(A)
+    A_sparse isa AbstractSparseMatrixCSC ||
+        error("HSLMA57Factorization currently supports only sparse CSC matrices")
+    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")
+
+    hcache = LinearSolve.@get_cacheval(cache, :HSLMA57Factorization)
+    hcache === nothing && error(
+        "HSLMA57Factorization supports `AbstractSparseMatrixCSC{<:Union{Float32, Float64}}`"
+    )
+
+    try
+        if cache.isfresh
+            hcache.ma57 = HSL.Ma57(A_sparse; alg.kwargs...)
+            HSL.ma57_factorize!(hcache.ma57)
+            cache.isfresh = false
+        end
+
+        _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
+        )
+    catch err
+        if err isa HSL.Ma57Exception
+            @SciMLMessage(
+                "MA57 failed: $(err.msg) (flag $(err.flag))",
+                cache.verbose,
+                :solver_failure
+            )
+            cache.isfresh = true
+            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_sparse = _sparse_csc_matrix(A)
+    A_sparse isa AbstractSparseMatrixCSC ||
+        error("HSLMA97Factorization currently supports only sparse CSC matrices")
+    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")
+
+    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_sparse)
+            hcache.ma97 = HSL.Ma97(A_int; alg.kwargs...)
+            HSL.ma97_factorize!(hcache.ma97, matrix_type = alg.matrix_type)
+            cache.isfresh = false
+        end
+
+        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
+        )
+    catch err
+        if err isa HSL.Ma97Exception
+            @SciMLMessage(
+                "MA97 failed: $(err.msg) (flag $(err.flag))",
+                cache.verbose,
+                :solver_failure
+            )
+            cache.isfresh = true
+            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/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 [
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..b954d36f8
--- /dev/null
+++ b/test/hsl/hsl.jl
@@ -0,0 +1,74 @@
+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)
+        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)
+        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)
+        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)
+        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
+
+        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
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")
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"