Add HSL wrapper

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"

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

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
 

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)
 

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 [

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"