From cd71ea5760b17ca6cead35dd73de7ebfbed0616f Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 27 May 2026 11:35:20 -0400 Subject: [PATCH 1/2] Add SciML test setup and correctness tests Restructure test/ around a GROUP env split (Core/QA, plus per-factorization subgroups), add Aqua + ExplicitImports QA, SafeTestsets isolation, and a GitHub Actions CI matrix. Substantially expands correctness coverage for cs_qr, cs_lu, cs_cholesky: known-answer matrices, identity / diagonal / tridiagonal cases, n=1 trivial case, input-matrix immutability, rhs immutability, multiple solves on one factorization, idempotent explicit finalize, cross-factorization consistency (QR vs LU on non-singular, QR vs LU vs Cholesky on SPD), LU pivoting on a zero-diagonal matrix, Hermitian (non-symmetric) Cholesky with complex off-diagonals, and overdetermined least-squares coefficient match against dense backslash. Two QA-driven source/compat fixes: - src/CXSparse.jl: `using LinearAlgebra` -> explicit `using LinearAlgebra: LinearAlgebra, ldiv!` (ExplicitImports). - Project.toml: add compat for Random and Test (Aqua). Bump julia compat 1.10 -> 1.11 to match the actual CXSparse_jll v4 SuiteSparse_jll requirement (resolver fails on 1.10). 626 tests pass on Julia 1.11. Co-Authored-By: Chris Rackauckas --- .github/workflows/CI.yml | 51 ++++++++++++ Project.toml | 12 ++- src/CXSparse.jl | 2 +- test/cholesky_tests.jl | 163 ++++++++++++++++++++++++++++++++++++ test/cross_tests.jl | 76 +++++++++++++++++ test/lu_tests.jl | 154 ++++++++++++++++++++++++++++++++++ test/qa.jl | 17 ++++ test/qr_tests.jl | 168 +++++++++++++++++++++++++++++++++++++ test/runtests.jl | 175 +++------------------------------------ test/shared.jl | 35 ++++++++ 10 files changed, 688 insertions(+), 165 deletions(-) create mode 100644 .github/workflows/CI.yml create mode 100644 test/cholesky_tests.jl create mode 100644 test/cross_tests.jl create mode 100644 test/lu_tests.jl create mode 100644 test/qa.jl create mode 100644 test/qr_tests.jl create mode 100644 test/shared.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..af1f534 --- /dev/null +++ b/.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 diff --git a/Project.toml b/Project.toml index 7c2d694..e946ce1 100644 --- a/Project.toml +++ b/Project.toml @@ -9,14 +9,22 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] +Aqua = "0.8" CXSparse_jll = "4" +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"] diff --git a/src/CXSparse.jl b/src/CXSparse.jl index 86ef1e6..373bd67 100644 --- a/src/CXSparse.jl +++ b/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 diff --git a/test/cholesky_tests.jl b/test/cholesky_tests.jl new file mode 100644 index 0000000..a7a2582 --- /dev/null +++ b/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 diff --git a/test/cross_tests.jl b/test/cross_tests.jl new file mode 100644 index 0000000..3ba0051 --- /dev/null +++ b/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 diff --git a/test/lu_tests.jl b/test/lu_tests.jl new file mode 100644 index 0000000..13eb4a3 --- /dev/null +++ b/test/lu_tests.jl @@ -0,0 +1,154 @@ +include("shared.jl") + +@testset "CXSparse cs_lu" begin + Random.seed!(0xCC55_F00D) + + @testset "square well-conditioned ($Tv, $Ti, n=$n)" for Tv in ELTYPES, + Ti in IDXTYPES, n in (5, 25, 100) + Adense = _randmat(Tv, n, n) + Tv(n) * I + A = _convert(Adense, Tv, Ti) + b = _randvec(Tv, n) + + F = cs_lu(A) + @test size(F) == (n, n) + x = F \ b + @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_lu(A) + A2 = _convert(_randmat(Tv, 5, 3), Tv, Ti) + @test_throws ErrorException cs_lu(A2) + end + + @testset "singular matrix is rejected ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + # A rank-deficient square matrix: row 2 == 2 × row 1. + A = _convert(Tv[1 2; 2 4], Tv, Ti) + @test_throws ErrorException cs_lu(A) + end + + @testset "dimension checks ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert(Tv[2 0; 0 3], Tv, Ti) + F = cs_lu(A) + @test_throws DimensionMismatch (F \ Tv[1, 2, 3]) + @test_throws DimensionMismatch ldiv!(zeros(Tv, 3), F, Tv[1, 2]) + @test_throws DimensionMismatch ldiv!(zeros(Tv, 2), F, Tv[1]) + end + + @testset "n=1 trivial case ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert(reshape(Tv[5.0], 1, 1), Tv, Ti) + F = cs_lu(A) + x = F \ Tv[10.0] + @test x ≈ 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_lu(A) + @test F \ b ≈ b + end + + @testset "diagonal matrix ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + d = Tv <: Complex ? Tv[2+0im, -3+1im, 4-2im, 1+0im] : Tv[2, -3, 4, 1] + A = _convert(Diagonal(d), Tv, Ti) + b = _randvec(Tv, 4) + F = cs_lu(A) + @test F \ b ≈ b ./ d + end + + @testset "requires pivoting (zero on diagonal) ($Tv, $Ti)" for + Tv in ELTYPES, Ti in IDXTYPES + # Without pivoting, LU breaks at A[1,1] = 0. Partial pivoting must swap + # rows. Exact solution exists: x = [1, 1]. + A = _convert(Tv[0 2; 1 0], Tv, Ti) + b = Tv[2, 1] + F = cs_lu(A) + @test F \ b ≈ Tv[1, 1] + end + + @testset "upper-triangular ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + # A = [2 1 0; 0 3 1; 0 0 4]; b = [5, 7, 8] ⇒ back-sub: x3=2, x2=5/3, x1=5/3 + A = _convert(Tv[2 1 0; 0 3 1; 0 0 4], Tv, Ti) + b = Tv[5, 7, 8] + F = cs_lu(A) + @test F \ b ≈ Matrix(A) \ b + end + + @testset "tridiagonal ($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_lu(A) + @test F \ b ≈ Matrix(Adense) \ b rtol=1e-10 + end + + @testset "input matrix is not mutated by factorization or solve ($Tv, $Ti)" for + Tv in ELTYPES, Ti in IDXTYPES + A = _convert(_randmat(Tv, 8, 8) + Tv(8) * I, Tv, Ti) + snap = _snapshot(A) + F = cs_lu(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 + A = _convert(_randmat(Tv, 6, 6) + Tv(6) * I, Tv, Ti) + b = _randvec(Tv, 6) + b_orig = copy(b) + F = cs_lu(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 + A = _convert(_randmat(Tv, 10, 10) + Tv(10) * I, Tv, Ti) + F = cs_lu(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 + A = _convert(_randmat(Tv, 5, 5) + Tv(5) * I, Tv, Ti) + F = cs_lu(A) + b = _randvec(Tv, 5) + x = Vector{Tv}(undef, 5) + ret = ldiv!(x, F, b) + @test ret === x + end + + @testset "small known-answer ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + # A = [2 1; 1 3]; b = [3, 4] ⇒ det=5, x = [(3*3 - 1*4)/5, (2*4 - 1*3)/5] = [1, 1] + A = _convert(Tv[2 1; 1 3], Tv, Ti) + b = Tv[3, 4] + F = cs_lu(A) + @test F \ b ≈ Tv[1, 1] + end + + @testset "explicit finalize is safe and idempotent" begin + A = sparse(randn(5, 5) + 5 * I) + F = cs_lu(A) + _ = F \ randn(5) + finalize(F) + finalize(F) + @test F.S == C_NULL + @test F.N == C_NULL + end +end diff --git a/test/qa.jl b/test/qa.jl new file mode 100644 index 0000000..4896984 --- /dev/null +++ b/test/qa.jl @@ -0,0 +1,17 @@ +using CXSparse +using Aqua +using ExplicitImports +using Test + +@testset "Aqua" begin + Aqua.test_all(CXSparse; ambiguities = false, piracies = false) + # Ambiguities are checked package-locally to avoid noise from deps like + # SparseArrays exposing ambiguous fallbacks not under our control. + Aqua.test_ambiguities(CXSparse; recursive = false) +end + +@testset "ExplicitImports" begin + @test check_no_implicit_imports(CXSparse) === nothing + @test check_no_stale_explicit_imports(CXSparse) === nothing + @test check_all_qualified_accesses_via_owners(CXSparse) === nothing +end diff --git a/test/qr_tests.jl b/test/qr_tests.jl new file mode 100644 index 0000000..63b07fa --- /dev/null +++ b/test/qr_tests.jl @@ -0,0 +1,168 @@ +include("shared.jl") + +@testset "CXSparse cs_qr" begin + Random.seed!(0xCC55_F00D) + + @testset "well-conditioned square ($Tv, $Ti, n=$n)" for Tv in ELTYPES, + Ti in IDXTYPES, n in (5, 25, 100) + + Adense = _randmat(Tv, n, n) + Tv(n) * I + A = _convert(Adense, Tv, Ti) + b = _randvec(Tv, n) + + F = cs_qr(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 "overdetermined m > n: matches dense LS ($Tv, $Ti, m=$m, n=$n)" for + Tv in ELTYPES, Ti in IDXTYPES, (m, n) in ((5, 3), (10, 4), (50, 10)) + + Adense = _randmat(Tv, m, n) + A = _convert(Adense, Tv, Ti) + b = _randvec(Tv, m) + + F = cs_qr(A) + @test size(F) == (m, n) + x_cs = F \ b + x_ref = Adense \ b + # Coefficient-wise match (not just residual norm) — both algorithms + # should arrive at the same least-squares minimizer. + @test x_cs ≈ x_ref rtol=1e-8 + @test norm(A * x_cs - b) ≈ norm(Adense * x_ref - b) rtol=1e-8 + end + + @testset "rejects underdetermined m < n ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert([1.0 2.0 3.0; 4.0 5.0 6.0], Tv, Ti) + @test_throws ErrorException cs_qr(A) + end + + @testset "dimension checks ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert(Tv[1 0; 1 1], Tv, Ti) + F = cs_qr(A) + @test_throws DimensionMismatch (F \ Tv[1, 2, 3]) + @test_throws DimensionMismatch ldiv!(zeros(Tv, 3), F, Tv[1, 2]) + @test_throws DimensionMismatch ldiv!(zeros(Tv, 2), F, Tv[1]) + end + + @testset "rank-deficient: solve does not throw" begin + # cs_qr is NOT rank-revealing — entries may be non-finite. Assert only + # that the call sequence runs to completion without error. + A = sparse([1.0 2.0; 0.5 1.0; 2.0 4.0]) # rank 1 + b = [1.0, 1.0, 1.0] + F = cs_qr(A) + x = F \ b + @test x isa Vector + @test length(x) == size(F, 2) + end + + @testset "n=1 trivial case ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert(reshape(Tv[3.0], 1, 1), Tv, Ti) + F = cs_qr(A) + x = F \ Tv[6.0] + @test x ≈ 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_qr(A) + x = F \ b + @test x ≈ b + end + + @testset "diagonal matrix ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + d = Tv <: Complex ? Tv[2+0im, -3+1im, 4-2im, 1+0im] : Tv[2, -3, 4, 1] + A = _convert(Diagonal(d), Tv, Ti) + b = _randvec(Tv, 4) + F = cs_qr(A) + x = F \ b + @test x ≈ b ./ d + end + + @testset "tridiagonal SPD-like ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + n = 20 + 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_qr(A) + x = F \ b + @test x ≈ Matrix(Adense) \ b rtol=1e-10 + end + + @testset "input matrix is not mutated by factorization or solve ($Tv, $Ti)" for + Tv in ELTYPES, Ti in IDXTYPES + A = _convert(_randmat(Tv, 8, 8) + Tv(8) * I, Tv, Ti) + snap = _snapshot(A) + F = cs_qr(A) + @test _unchanged(A, snap) + b = _randvec(Tv, 8) + _ = F \ b + @test _unchanged(A, snap) + end + + @testset "rhs `b` is not mutated by solve ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + A = _convert(_randmat(Tv, 6, 6) + Tv(6) * I, Tv, Ti) + b = _randvec(Tv, 6) + b_orig = copy(b) + F = cs_qr(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 + A = _convert(_randmat(Tv, 10, 10) + Tv(10) * I, Tv, Ti) + F = cs_qr(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 "small known-answer matrix ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES + # A upper-triangular: back-substitution by hand. + # [2 1; 0 3] x = [4; 6] ⇒ x2 = 2, x1 = (4 - 1*2)/2 = 1 + A = _convert(Tv[2 1; 0 3], Tv, Ti) + b = Tv[4, 6] + F = cs_qr(A) + @test F \ b ≈ Tv[1.0, 2.0] + end + + @testset "ldiv! returns the destination vector ($Tv, $Ti)" for + Tv in ELTYPES, Ti in IDXTYPES + A = _convert(_randmat(Tv, 5, 5) + Tv(5) * I, Tv, Ti) + F = cs_qr(A) + b = _randvec(Tv, 5) + x = Vector{Tv}(undef, 5) + ret = ldiv!(x, F, b) + @test ret === x + end + + @testset "explicit finalize is safe and idempotent" begin + A = sparse(randn(5, 5) + 5 * I) + F = cs_qr(A) + b = randn(5) + _ = F \ b + finalize(F) + # A second finalize call must not crash (pointers should be NULL'd). + finalize(F) + @test F.S == C_NULL + @test F.N == C_NULL + end +end diff --git a/test/runtests.jl b/test/runtests.jl index fa38cc4..42c65be 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,174 +1,25 @@ -using CXSparse -using LinearAlgebra -using SparseArrays -using Test -using Random +using SafeTestsets -const ELTYPES = (Float64, ComplexF64) -const IDXTYPES = (Int32, Int64) +const GROUP = get(ENV, "GROUP", "All") -# Convert a generic sparse matrix to the (Tv, Ti) variant used by a particular -# CXSparse entry point. CXSparse_jll only ships Float64/ComplexF64 × Int32/Int64. -_convert(A::SparseMatrixCSC, ::Type{Tv}, ::Type{Ti}) where {Tv,Ti} = - SparseMatrixCSC{Tv,Ti}(A) - -@testset "CXSparse cs_qr" begin - Random.seed!(0xCC55_F00D) - - @testset "well-conditioned square ($Tv, $Ti, n=$n)" for Tv in ELTYPES, - Ti in IDXTYPES, n in (5, 25, 100) - - # Sparse-ish, diagonally dominant → non-singular. - Adense = (Tv <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n)) + - Tv(n) * I - A = _convert(sparse(Adense), Tv, Ti) - b = (Tv <: Complex ? complex.(randn(n), randn(n)) : randn(n)) - b = Vector{Tv}(b) - - F = cs_qr(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) - - x_pre = zeros(Tv, n) - @test ldiv!(x_pre, F, b) === x_pre - @test x_pre ≈ x - end - - @testset "overdetermined m > n ($Tv, $Ti)" for Tv in ELTYPES, Ti in IDXTYPES - # 5×3 overdetermined system; check that the residual matches the - # dense least-squares solution to high precision. - m, n = 5, 3 - Adense = Tv <: Complex ? complex.(randn(m, n), randn(m, n)) : randn(m, n) - b_ = Tv <: Complex ? complex.(randn(m), randn(m)) : randn(m) - b = Vector{Tv}(b_) - A = _convert(sparse(Adense), Tv, Ti) - - F = cs_qr(A) - @test size(F) == (m, n) - x_cs = F \ b - x_ref = Adense \ b - @test norm(A * x_cs - b) ≈ norm(Adense * x_ref - b) rtol=1e-8 - end - - @testset "rejects underdetermined m < n" for Tv in ELTYPES, Ti in IDXTYPES - A = _convert(sparse([1.0 2.0 3.0; 4.0 5.0 6.0]), Tv, Ti) - @test_throws ErrorException cs_qr(A) - end - - @testset "dimension checks" begin - A = sparse([1.0 0.0; 1.0 1.0]) - F = cs_qr(A) - @test_throws DimensionMismatch (F \ [1.0, 2.0, 3.0]) - @test_throws DimensionMismatch ldiv!(zeros(3), F, [1.0, 2.0]) - @test_throws DimensionMismatch ldiv!(zeros(2), F, [1.0]) +@time begin + if GROUP == "All" || GROUP == "Core" || GROUP == "QR" + @safetestset "cs_qr" include("qr_tests.jl") end - @testset "rank-deficient: solve does not throw" begin - # cs_qr is NOT rank-revealing — `x` may contain non-finite entries on - # truly rank-deficient matrices. We only assert that calling factor + - # solve does not throw. Use SPQR if you need a rank-revealing result. - A = sparse([1.0 2.0; 0.5 1.0; 2.0 4.0]) # rank 1 - b = [1.0, 1.0, 1.0] - F = cs_qr(A) - x = F \ b - @test x isa Vector - @test length(x) == size(F, 2) + if GROUP == "All" || GROUP == "Core" || GROUP == "LU" + @safetestset "cs_lu" include("lu_tests.jl") end -end - -@testset "CXSparse cs_lu" begin - Random.seed!(0xCC55_F00D) - - @testset "square well-conditioned ($Tv, $Ti, n=$n)" for Tv in ELTYPES, - Ti in IDXTYPES, n in (5, 25, 100) - Adense = (Tv <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n)) + - Tv(n) * I - A = _convert(sparse(Adense), Tv, Ti) - b_ = Tv <: Complex ? complex.(randn(n), randn(n)) : randn(n) - b = Vector{Tv}(b_) - F = cs_lu(A) - @test size(F) == (n, n) - x = F \ b - @test norm(A * x - b) < 1e-8 * norm(b) - x_pre = zeros(Tv, n) - @test ldiv!(x_pre, F, b) === x_pre - @test x_pre ≈ x + if GROUP == "All" || GROUP == "Core" || GROUP == "Cholesky" + @safetestset "cs_cholesky" include("cholesky_tests.jl") end - @testset "rejects non-square" for Tv in ELTYPES, Ti in IDXTYPES - A = _convert(sparse(randn(3, 5)), Tv, Ti) - @test_throws ErrorException cs_lu(A) + if GROUP == "All" || GROUP == "Core" || GROUP == "Cross" + @safetestset "cross-factorization" include("cross_tests.jl") end -end - -@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) - - # SPD (real) / Hermitian PD (complex) via M*M' + n*I. - Mraw = Tv <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n) - Adense = Mraw * Mraw' + Tv(n) * I - # Adense is dense Float64 / ComplexF64 either way; sparsify and convert. - A = _convert(sparse(Adense), Tv, Ti) - b_ = Tv <: Complex ? complex.(randn(n), randn(n)) : randn(n) - b = Vector{Tv}(b_) - - 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) - - 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(sparse(randn(3, 5)), Tv, Ti) - @test_throws ErrorException cs_cholesky(A) - end - - @testset "rejects non-positive-definite" begin - # Negative-definite — cs_chol returns NULL on the first non-positive - # pivot, which we surface as an error. - A = sparse(Float64[-1.0 0.0; 0.0 -1.0]) - @test_throws ErrorException cs_cholesky(A) - end - - @testset "dimension checks" begin - A = sparse(Float64[2.0 0.0; 0.0 2.0]) - F = cs_cholesky(A) - @test_throws DimensionMismatch (F \ [1.0, 2.0, 3.0]) - @test_throws DimensionMismatch ldiv!(zeros(3), F, [1.0, 2.0]) - end -end -@testset "Finalizer doesn't crash" begin - # Force GC of factorizations and ensure we don't double-free. - for _ in 1:50 - A = sparse(randn(10, 10) + 10 * I) - Asym = Matrix(A) * Matrix(A)' + 10 * I - Asym = sparse(Asym) - F = cs_qr(A) - F2 = cs_lu(A) - F3 = cs_cholesky(Asym) - _ = F \ randn(10) - _ = F2 \ randn(10) - _ = F3 \ randn(10) + if GROUP == "All" || GROUP == "QA" + @safetestset "Quality Assurance" include("qa.jl") end - GC.gc() - GC.gc() - @test true end diff --git a/test/shared.jl b/test/shared.jl new file mode 100644 index 0000000..a84321d --- /dev/null +++ b/test/shared.jl @@ -0,0 +1,35 @@ +using CXSparse +using LinearAlgebra +using SparseArrays +using Test +using Random + +const ELTYPES = (Float64, ComplexF64) +const IDXTYPES = (Int32, Int64) + +# Convert a generic dense/sparse matrix to the (Tv, Ti) variant. CXSparse_jll +# ships only Float64/ComplexF64 × Int32/Int64. +_convert(A::AbstractMatrix, ::Type{Tv}, ::Type{Ti}) where {Tv,Ti} = + SparseMatrixCSC{Tv,Ti}(sparse(A)) +_convert(A::SparseMatrixCSC, ::Type{Tv}, ::Type{Ti}) where {Tv,Ti} = + SparseMatrixCSC{Tv,Ti}(A) + +# Random vector of the right element type. +_randvec(::Type{T}, n::Integer) where {T<:Real} = Vector{T}(randn(n)) +_randvec(::Type{T}, n::Integer) where {T<:Complex} = + Vector{T}(complex.(randn(n), randn(n))) + +# Random dense matrix of the right element type. +_randmat(::Type{T}, m::Integer, n::Integer) where {T<:Real} = randn(m, n) +_randmat(::Type{T}, m::Integer, n::Integer) where {T<:Complex} = + complex.(randn(m, n), randn(m, n)) + +# Snapshot the storage of a SparseMatrixCSC so we can verify it is not mutated +# by factorization. Returns deep copies of (colptr, rowval, nzval). +_snapshot(A::SparseMatrixCSC) = (copy(A.colptr), copy(A.rowval), copy(A.nzval)) + +# Compare a snapshot with the current state of A. +function _unchanged(A::SparseMatrixCSC, snap) + cp, rv, nz = snap + return A.colptr == cp && A.rowval == rv && A.nzval == nz +end From b0e9fc6d12e692b4d198c883f3402601d03895b1 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 27 May 2026 12:17:43 -0400 Subject: [PATCH 2/2] Loosen CXSparse_jll compat to "4, 400" MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CXSparse_jll has two distinct version series in the General registry: - 4.4.x — uses the upstream CXSparse 4.4 release version directly - 400.400.200 — uses the SuiteSparse-bundle version scheme (4.4.0 → 400.400.200) with looser SuiteSparse_jll compat that admits 7.7 and 7.8 The previous compat `"4"` selected only the 4.4.x series, which pins SuiteSparse_jll ~7.7. Julia 1.12 ships SuiteSparse_jll 7.8.3 in its stdlib, so resolution fails on 1.12 (CI was red on every Julia 1 / 1.12 job). Allowing both series ("4, 400") lets the resolver pick the appropriate one per Julia version. Tests pass on Julia 1.11 (v4.4.0) and Julia 1.12 (v400.400.200), 626 tests each. Co-Authored-By: Chris Rackauckas --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e946ce1..af3fc84 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] Aqua = "0.8" -CXSparse_jll = "4" +CXSparse_jll = "4, 400" ExplicitImports = "1" LinearAlgebra = "1" Random = "1"