Add SciML test setup and correctness tests

.github/workflows/CI.yml

@@ -0,0 +1,51 @@
+name: CI
+on:
+  push:
+    branches:
+      - main
+    tags: ['*']
+  pull_request:
+  workflow_dispatch:
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}-${{ github.head_ref || github.run_id }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.group }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 60
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.11'
+          - '1'
+        os:
+          - ubuntu-latest
+        group:
+          - Core
+          - QA
+        include:
+          - version: '1'
+            os: macos-latest
+            group: Core
+          - version: '1'
+            os: windows-latest
+            group: Core
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: ${{ matrix.version }}
+          arch: x64
+      - uses: julia-actions/cache@v2
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        env:
+          GROUP: ${{ matrix.group }}
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v4
+        with:
+          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false

Project.toml

@@ -9,14 +9,22 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-CXSparse_jll = "4"
+Aqua = "0.8"
+CXSparse_jll = "4, 400"
+ExplicitImports = "1"
 LinearAlgebra = "1"
+Random = "1"
+SafeTestsets = "0.1"
 SparseArrays = "1"
-julia = "1.10"
+Test = "1"
+julia = "1.11"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Random", "Test"]
+test = ["Aqua", "ExplicitImports", "Random", "SafeTestsets", "Test"]

src/CXSparse.jl

@@ -1,7 +1,7 @@
 module CXSparse
 
 using CXSparse_jll: libcxsparse
-using LinearAlgebra
+using LinearAlgebra: LinearAlgebra, ldiv!
 using SparseArrays: SparseArrays, SparseMatrixCSC, getcolptr, rowvals, nonzeros
 
 export cs_qr, cs_lu, cs_cholesky, CSQR, CSLU, CSCholesky

test/cholesky_tests.jl

@@ -0,0 +1,163 @@
+include("shared.jl")
+
+@testset "CXSparse cs_cholesky" begin
+    Random.seed!(0xCC55_F00D)
+
+    @testset "SPD ($Tv, $Ti, n=$n)" for Tv in ELTYPES,
+            Ti in IDXTYPES, n in (5, 25, 100)
+
+        Mraw = _randmat(Tv, n, n)
+        Adense = Mraw * Mraw' + Tv(n) * I
+        A = _convert(Adense, Tv, Ti)
+        b = _randvec(Tv, n)
+
+        F = cs_cholesky(A)
+        @test size(F) == (n, n)
+        @test size(F, 1) == n
+        @test size(F, 2) == n
+
+        x = F \ b
+        @test length(x) == n
+        @test norm(A * x - b) < 1e-8 * norm(b)
+        @test x ≈ Adense \ b rtol=1e-8
+
+        x_pre = zeros(Tv, n)
+        @test ldiv!(x_pre, F, b) === x_pre
+        @test x_pre ≈ x
+    end
+
+    @testset "rejects non-square ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        A = _convert(_randmat(Tv, 3, 5), Tv, Ti)
+        @test_throws ErrorException cs_cholesky(A)
+    end
+
+    @testset "rejects non-positive-definite ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        A = _convert(Tv[-1 0; 0 -1], Tv, Ti)
+        @test_throws ErrorException cs_cholesky(A)
+    end
+
+    @testset "rejects zero diagonal ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        A = _convert(Tv[0 0; 0 1], Tv, Ti)
+        @test_throws ErrorException cs_cholesky(A)
+    end
+
+    @testset "dimension checks ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        A = _convert(Tv[2 0; 0 2], Tv, Ti)
+        F = cs_cholesky(A)
+        @test_throws DimensionMismatch (F \ Tv[1, 2, 3])
+        @test_throws DimensionMismatch ldiv!(zeros(Tv, 3), F, Tv[1, 2])
+    end
+
+    @testset "n=1 trivial case ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        A = _convert(reshape(Tv[9.0], 1, 1), Tv, Ti)
+        F = cs_cholesky(A)
+        @test F \ Tv[18.0] ≈ Tv[2.0]
+    end
+
+    @testset "identity matrix ($Tv, $Ti, n=$n)" for Tv in ELTYPES,
+            Ti in IDXTYPES, n in (3, 10)
+        A = _convert(Matrix{Tv}(I, n, n), Tv, Ti)
+        b = _randvec(Tv, n)
+        F = cs_cholesky(A)
+        @test F \ b ≈ b
+    end
+
+    @testset "diagonal positive matrix ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        d = Tv <: Complex ? Tv[2+0im, 3+0im, 4+0im, 1+0im] : Tv[2, 3, 4, 1]
+        A = _convert(Diagonal(d), Tv, Ti)
+        b = _randvec(Tv, 4)
+        F = cs_cholesky(A)
+        @test F \ b ≈ b ./ d
+    end
+
+    @testset "tridiagonal SPD ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        n = 25
+        Adense = Tridiagonal(fill(Tv(-1), n - 1), fill(Tv(4), n), fill(Tv(-1), n - 1))
+        A = _convert(Matrix(Adense), Tv, Ti)
+        b = _randvec(Tv, n)
+        F = cs_cholesky(A)
+        @test F \ b ≈ Matrix(Adense) \ b rtol=1e-10
+    end
+
+    @testset "Hermitian (not just real-symmetric) ($Ti)" for Ti in IDXTYPES
+        # Hermitian PD matrix with non-trivial complex off-diagonals.
+        Adense = ComplexF64[
+            4.0+0.0im   1.0+0.5im   0.0+0.0im
+            1.0-0.5im   3.0+0.0im   0.5-0.25im
+            0.0+0.0im   0.5+0.25im  2.0+0.0im
+        ]
+        @assert ishermitian(Adense)
+        A = _convert(Adense, ComplexF64, Ti)
+        b = ComplexF64[1+0im, 2-1im, -1+0.5im]
+        F = cs_cholesky(A)
+        x = F \ b
+        @test x ≈ Adense \ b rtol=1e-10
+        @test norm(A * x - b) < 1e-10 * norm(b)
+    end
+
+    @testset "small known-answer ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        # A = [4 2; 2 5]  (SPD: dets are 4, 16) ; b = [6, 7]
+        # det(A) = 16, x = [(5*6 - 2*7)/16, (4*7 - 2*6)/16] = [16/16, 16/16] = [1, 1]
+        A = _convert(Tv[4 2; 2 5], Tv, Ti)
+        b = Tv[6, 7]
+        F = cs_cholesky(A)
+        @test F \ b ≈ Tv[1, 1]
+    end
+
+    @testset "input matrix is not mutated by factorization or solve ($Tv, $Ti)" for
+            Tv in ELTYPES, Ti in IDXTYPES
+        M = _randmat(Tv, 8, 8)
+        A = _convert(M * M' + Tv(8) * I, Tv, Ti)
+        snap = _snapshot(A)
+        F = cs_cholesky(A)
+        @test _unchanged(A, snap)
+        _ = F \ _randvec(Tv, 8)
+        @test _unchanged(A, snap)
+    end
+
+    @testset "rhs `b` is not mutated by solve ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES
+        M = _randmat(Tv, 6, 6)
+        A = _convert(M * M' + Tv(6) * I, Tv, Ti)
+        b = _randvec(Tv, 6)
+        b_orig = copy(b)
+        F = cs_cholesky(A)
+        _ = F \ b
+        @test b == b_orig
+        _ = ldiv!(zeros(Tv, 6), F, b)
+        @test b == b_orig
+    end
+
+    @testset "multiple solves on same factorization ($Tv, $Ti)" for
+            Tv in ELTYPES, Ti in IDXTYPES
+        M = _randmat(Tv, 10, 10)
+        A = _convert(M * M' + Tv(10) * I, Tv, Ti)
+        F = cs_cholesky(A)
+        Adense = Matrix(A)
+        for _ in 1:5
+            b = _randvec(Tv, 10)
+            x = F \ b
+            @test x ≈ Adense \ b rtol=1e-10
+        end
+    end
+
+    @testset "ldiv! returns the destination vector ($Tv, $Ti)" for
+            Tv in ELTYPES, Ti in IDXTYPES
+        M = _randmat(Tv, 5, 5)
+        A = _convert(M * M' + Tv(5) * I, Tv, Ti)
+        F = cs_cholesky(A)
+        b = _randvec(Tv, 5)
+        x = Vector{Tv}(undef, 5)
+        @test ldiv!(x, F, b) === x
+    end
+
+    @testset "explicit finalize is safe and idempotent" begin
+        M = randn(5, 5)
+        A = sparse(M * M' + 5 * I)
+        F = cs_cholesky(A)
+        _ = F \ randn(5)
+        finalize(F)
+        finalize(F)
+        @test F.S == C_NULL
+        @test F.N == C_NULL
+    end
+end

test/cross_tests.jl

@@ -0,0 +1,76 @@
+include("shared.jl")
+
+@testset "Cross-factorization consistency" begin
+    Random.seed!(0xCC55_F00D)
+
+    # On a well-conditioned square non-singular matrix, cs_qr and cs_lu must
+    # produce solutions that match each other (and the dense reference) to
+    # high precision. On an SPD matrix, cs_cholesky must also agree.
+
+    @testset "QR vs LU on square non-singular ($Tv, $Ti, n=$n)" for
+            Tv in ELTYPES, Ti in IDXTYPES, n in (8, 30)
+        Adense = _randmat(Tv, n, n) + Tv(n) * I
+        A = _convert(Adense, Tv, Ti)
+        b = _randvec(Tv, n)
+
+        x_qr = cs_qr(A) \ b
+        x_lu = cs_lu(A) \ b
+        x_ref = Adense \ b
+        @test x_qr ≈ x_lu rtol=1e-9
+        @test x_qr ≈ x_ref rtol=1e-9
+        @test x_lu ≈ x_ref rtol=1e-9
+    end
+
+    @testset "QR vs LU vs Cholesky on SPD ($Tv, $Ti, n=$n)" for
+            Tv in ELTYPES, Ti in IDXTYPES, n in (8, 30)
+        M = _randmat(Tv, n, n)
+        Adense = M * M' + Tv(n) * I
+        A = _convert(Adense, Tv, Ti)
+        b = _randvec(Tv, n)
+
+        x_qr = cs_qr(A) \ b
+        x_lu = cs_lu(A) \ b
+        x_chol = cs_cholesky(A) \ b
+        x_ref = Adense \ b
+        @test x_qr ≈ x_lu rtol=1e-9
+        @test x_lu ≈ x_chol rtol=1e-9
+        @test x_chol ≈ x_ref rtol=1e-9
+    end
+end
+
+@testset "Finalizer GC stress" begin
+    # Repeatedly construct and discard factorizations; force GC to run finalizers
+    # and ensure we don't double-free or crash.
+    for _ in 1:50
+        A = sparse(randn(10, 10) + 10 * I)
+        M = randn(10, 10)
+        Asym = sparse(M * M' + 10 * I)
+        F = cs_qr(A)
+        F2 = cs_lu(A)
+        F3 = cs_cholesky(Asym)
+        _ = F \ randn(10)
+        _ = F2 \ randn(10)
+        _ = F3 \ randn(10)
+    end
+    GC.gc()
+    GC.gc()
+    @test true
+end
+
+@testset "Mixed element/index types interoperate" begin
+    # Each (Tv, Ti) variant should work standalone in sequence (the @ccall
+    # dispatch picks the correct cs_* symbol from libcxsparse based on the
+    # generic-function method we defined).
+    for Tv in (Float64, ComplexF64), Ti in (Int32, Int64)
+        Adense = Tv <: Complex ?
+            complex.(randn(6, 6), randn(6, 6)) + Tv(6) * I :
+            randn(6, 6) + Tv(6) * I
+        A = SparseMatrixCSC{Tv,Ti}(sparse(Adense))
+        b = Tv <: Complex ?
+            Vector{Tv}(complex.(randn(6), randn(6))) :
+            Vector{Tv}(randn(6))
+        F = cs_qr(A)
+        x = F \ b
+        @test norm(A * x - b) < 1e-8 * norm(b)
+    end
+end