Add CUSOLVERRF.jl integration for GPU-accelerated sparse LU factorization by ChrisRackauckas · Pull Request #651 · SciML/LinearSolve.jl

ChrisRackauckas · 2025-07-28T19:42:49Z

Summary

This PR adds support for NVIDIA's cusolverRF sparse LU factorization library through a package extension, providing high-performance GPU-accelerated solving for sparse linear systems.

Motivation

CUSOLVERRF.jl provides access to NVIDIA's cusolverRF library, which offers significant performance improvements for sparse LU factorization on GPUs. This integration makes it accessible through LinearSolve.jl's unified interface.

Key Features

New CUSOLVERRFFactorization algorithm with configurable options:
- symbolic: Choose between :RF (default) or :KLU for symbolic factorization
- reuse_symbolic: Reuse symbolic factorization for matrices with same sparsity pattern
Automatic CPU-to-GPU conversion for convenience
Support for multiple right-hand sides
Adjoint solve support
Comprehensive test suite

Implementation Details

The implementation follows LinearSolve.jl's extension pattern:

Extension module in ext/LinearSolveCUSOLVERRFExt.jl
Core types and exports in src/factorization.jl and src/LinearSolve.jl
Weak dependency configuration in Project.toml
Tests in test/gpu/cusolverrf.jl

Usage Example

using LinearSolve, CUSOLVERRF, SparseArrays

# Create sparse system
A = sprand(1000, 1000, 0.01) + 5I
b = rand(1000)

# Solve with default options
prob = LinearProblem(A, b)
sol = solve(prob, CUSOLVERRFFactorization())

# Use KLU for symbolic factorization
sol = solve(prob, CUSOLVERRFFactorization(symbolic = :KLU))

Limitations

Only supports Float64 element types with Int32 indices (CUSOLVERRF limitation)
Requires CUDA-capable GPU

Testing

Tests have been added to the GPU test suite and can be run with appropriate hardware.

🤖 Generated with Claude Code

github-actions

Remaining comments which cannot be posted as a review comment to avoid GitHub Rate Limit

JuliaFormatter

[JuliaFormatter] _{reported by reviewdog 🐶}

LinearSolve.jl/test/gpu/cusolverrf.jl

Line 84 in 35073d8

[JuliaFormatter] _{reported by reviewdog 🐶}

LinearSolve.jl/test/gpu/cusolverrf.jl

Line 88 in 35073d8

end

github-actions · 2025-07-28T19:46:07Z

+        @info "CUDA not available, skipping CUSOLVERRF tests"
+        return
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:07Z

+    n = 100
+    A = sprand(n, n, 0.1) + I
+    b = rand(n)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:07Z

+    # Test with CPU sparse matrix (should auto-convert to GPU)
+    @testset "CPU Sparse Matrix" begin
+        prob = LinearProblem(A, b)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:07Z

+        # Test with default symbolic (:RF)
+        sol = solve(prob, CUSOLVERRFFactorization())
+        @test norm(A * sol.u - b) / norm(b) < 1e-10
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:07Z

+        sol_klu = solve(prob, CUSOLVERRFFactorization(symbolic = :KLU))
+        @test norm(A * sol_klu.u - b) / norm(b) < 1e-10
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:08Z

+    @testset "GPU Sparse Matrix" begin
+        A_gpu = CUDA.CUSPARSE.CuSparseMatrixCSR(A)
+        b_gpu = CuArray(b)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:08Z

+
+        prob_gpu = LinearProblem(A_gpu, b_gpu)
+        sol_gpu = solve(prob_gpu, CUSOLVERRFFactorization())
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:08Z

+        res_gpu = A_gpu * sol_gpu.u - b_gpu
+        @test norm(res_gpu) / norm(b_gpu) < 1e-10
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:08Z

+        # Create a new matrix with same pattern but different values
+        A2 = A + 0.1 * sprand(n, n, 0.01)
+        b2 = rand(n)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-07-28T19:46:08Z

+        sol2 = solve(prob2, CUSOLVERRFFactorization(reuse_symbolic = true))
+        @test norm(A2 * sol2.u - b2) / norm(b2) < 1e-10
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

oscardssmith · 2025-07-28T19:46:29Z

symbolic: Choose between :RF (default) or :KLU for symbolic factorization

Should these be types rather than symbols for type stability?

ChrisRackauckas · 2025-07-28T19:54:31Z

no it's an upstream thing https://github.com/exanauts/CUSOLVERRF.jl/blob/88b5e242e04829845cc8e06f62cefe8575cac368/README.md?plain=1#L62-L65

…tion This PR adds support for NVIDIA's cusolverRF sparse LU factorization library through a package extension. CUSOLVERRF provides high-performance GPU-accelerated factorization for sparse matrices. Key features: - New `CUSOLVERRFFactorization` algorithm with configurable symbolic factorization (RF or KLU) - Automatic CPU-to-GPU conversion for convenience - Support for multiple right-hand sides - Reusable symbolic factorization for matrices with same sparsity pattern - Adjoint solve support - Comprehensive test suite The implementation follows LinearSolve.jl's extension pattern, similar to the existing CUDSS integration. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <noreply@anthropic.com>

Include CUSOLVERRF tests in the GPU test suite when the package is available. The tests are conditionally included to avoid failures when CUSOLVERRF.jl is not installed. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <noreply@anthropic.com>

- Added CUSOLVERRF to recommended methods for sparse matrices - Added CUSOLVERRF section in the full list of solvers - Added CUSOLVERRF examples in GPU tutorial documentation - Documented supported options and limitations 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <noreply@anthropic.com>

- Updated sparse matrices recommendation to include both CUDSS.jl and CUSOLVERRF.jl - Clarified that CUDSS provides interface to NVIDIA's cuDSS library - Maintained that both offer high performance for GPU-accelerated sparse LU factorization 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <noreply@anthropic.com>

- Clarified that CUDSS works through LUFactorization() when CUDSS.jl is loaded - Explained that it automatically uses cuDSS for CuSparseMatrixCSR arrays - Removed incorrect reference to a separate CUDSS factorization type 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <noreply@anthropic.com>

github-actions · 2025-08-05T18:49:38Z

+CUSOLVERRF provides access to NVIDIA's cusolverRF library, which offers significant 
+performance improvements for sparse LU factorization on GPUs. It supports both 
+`:RF` (default) and `:KLU` symbolic factorization methods, and can reuse symbolic 


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

CUSOLVERRF provides access to NVIDIA's cusolverRF library, which offers significant

performance improvements for sparse LU factorization on GPUs. It supports both

`:RF` (default) and `:KLU` symbolic factorization methods, and can reuse symbolic

CUSOLVERRF provides access to NVIDIA's cusolverRF library, which offers significant

performance improvements for sparse LU factorization on GPUs. It supports both

`:RF` (default) and `:KLU` symbolic factorization methods, and can reuse symbolic

github-actions · 2025-08-05T23:25:45Z

+end
+
+function LinearSolve.init_cacheval(alg::LinearSolve.CUSOLVERRFFactorization,
+        A::Union{CuSparseMatrixCSR{Float64, Int32}, SparseMatrixCSC{Float64, <:Integer}}, 


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

A::Union{CuSparseMatrixCSR{Float64, Int32}, SparseMatrixCSC{Float64, <:Integer}},

A::Union{CuSparseMatrixCSR{Float64, Int32}, SparseMatrixCSC{Float64, <:Integer}},

github-actions · 2025-08-05T23:25:46Z

+    symbolic = alg.symbolic
+    # Convert to CuSparseMatrixCSR if needed
+    A_gpu = A isa CuSparseMatrixCSR ? A : CuSparseMatrixCSR(A)
+    RFLU(A_gpu; nrhs=nrhs, symbolic=symbolic)


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

RFLU(A_gpu; nrhs=nrhs, symbolic=symbolic)

RFLU(A_gpu; nrhs = nrhs, symbolic = symbolic)

github-actions · 2025-08-05T23:25:46Z

+
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::LinearSolve.CUSOLVERRFFactorization; kwargs...)
+    A = cache.A
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-05T23:25:46Z

+    else
+        error("CUSOLVERRFFactorization only supports SparseMatrixCSC or CuSparseMatrixCSR matrices")
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-05T23:25:46Z

+        if cacheval === nothing
+            # Create new factorization
+            nrhs = cache.b isa AbstractMatrix ? size(cache.b, 2) : 1
+            fact = RFLU(A_gpu; nrhs=nrhs, symbolic=alg.symbolic)


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

fact = RFLU(A_gpu; nrhs=nrhs, symbolic=alg.symbolic)

fact = RFLU(A_gpu; nrhs = nrhs, symbolic = alg.symbolic)

github-actions · 2025-08-05T23:25:46Z

+            else
+                # Create new factorization if pattern changed
+                nrhs = cache.b isa AbstractMatrix ? size(cache.b, 2) : 1
+                fact = RFLU(A_gpu; nrhs=nrhs, symbolic=alg.symbolic)


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

fact = RFLU(A_gpu; nrhs=nrhs, symbolic=alg.symbolic)

fact = RFLU(A_gpu; nrhs = nrhs, symbolic = alg.symbolic)

github-actions · 2025-08-05T23:25:46Z

+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-05T23:25:47Z

+    end
+
+    F = @get_cacheval(cache, :CUSOLVERRFFactorization)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-05T23:25:47Z

+    # Ensure b and u are on GPU
+    b_gpu = cache.b isa CUDA.CuArray ? cache.b : CUDA.CuArray(cache.b)
+    u_gpu = cache.u isa CUDA.CuArray ? cache.u : CUDA.CuArray(cache.u)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-06T00:12:11Z

+    # Solve
+    copyto!(u_gpu, b_gpu)
+    ldiv!(F, u_gpu)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-06T00:12:11Z

+
+    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)

SciMLBase.build_linear_solution(

alg, cache.u, nothing, cache; retcode = ReturnCode.Success)

github-actions · 2025-08-06T00:12:11Z

+
+end


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

end

end

github-actions · 2025-08-06T02:07:52Z

+        A_f32 = Float32.(A)
+        b_f32 = Float32.(b)
+        prob_f32 = LinearProblem(A_f32, b_f32)
+


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

github-actions · 2025-08-06T02:07:52Z

+        # This should error since CUSOLVERRF only supports Float64
+        @test_throws Exception solve(prob_f32, CUSOLVERRFFactorization())
+    end
+end


[JuliaFormatter] _{reported by reviewdog 🐶}

Suggested change

end

end

ChrisRackauckas commented Jul 28, 2025

View reviewed changes

Comment thread Project.toml Outdated

github-actions Bot reviewed Jul 28, 2025

View reviewed changes

ChrisRackauckas-Claude mentioned this pull request Aug 5, 2025

Add CUSOLVERRF.jl integration for GPU-accelerated sparse LU factorization #673

Merged

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread test/gpu/Project.toml Outdated

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread Project.toml Outdated

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread Project.toml Outdated

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread Project.toml

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread ext/LinearSolveCUSOLVERRFExt.jl Outdated

ChrisRackauckas commented Aug 5, 2025

View reviewed changes

Comment thread ext/LinearSolveCUSOLVERRFExt.jl Outdated

claude and others added 11 commits August 5, 2025 14:34

Update Project.toml

e40ad85

Update Project.toml

cc7911b

Update factorization.jl

235e333

Update extension_algs.jl

f784d42

Update solvers.md

0ac5d28

Update Project.toml

d7f1f8c

ChrisRackauckas force-pushed the add-cusolverrf-support branch from 8e3ce1a to d7f1f8c Compare August 5, 2025 18:35

Update src/extension_algs.jl

0a075fe

github-actions Bot reviewed Aug 5, 2025

View reviewed changes

ChrisRackauckas added 5 commits August 5, 2025 15:04

Update src/extension_algs.jl

1c1e917

Update Project.toml

b92906c

Update Project.toml

e88bad8

Update Project.toml

82fbc55

Update ext/LinearSolveCUSOLVERRFExt.jl

7a8dac7

ChrisRackauckas added 2 commits August 5, 2025 17:37

Update ext/LinearSolveCUSOLVERRFExt.jl

288d382

Update ext/LinearSolveCUSOLVERRFExt.jl

62bc9ae

github-actions Bot reviewed Aug 5, 2025

View reviewed changes

Update ext/LinearSolveCUSOLVERRFExt.jl

d559e8b

github-actions Bot reviewed Aug 6, 2025

View reviewed changes

ChrisRackauckas added 5 commits August 5, 2025 20:59

Update ext/LinearSolveCUSOLVERRFExt.jl

5175137

Update test/gpu/cusolverrf.jl

f1f3bb8

Update ext/LinearSolveCUSOLVERRFExt.jl

6db7c55

Update ext/LinearSolveCUSOLVERRFExt.jl

b8ca961

Update test/gpu/cusolverrf.jl

6a96db1

github-actions Bot reviewed Aug 6, 2025

View reviewed changes

ChrisRackauckas closed this Aug 9, 2025

ChrisRackauckas deleted the add-cusolverrf-support branch August 9, 2025 14:07


		prob_gpu = LinearProblem(A_gpu, b_gpu)
		sol_gpu = solve(prob_gpu, CUSOLVERRFFactorization())

	A::Union{CuSparseMatrixCSR{Float64, Int32}, SparseMatrixCSC{Float64, <:Integer}},
	A::Union{CuSparseMatrixCSR{Float64, Int32}, SparseMatrixCSC{Float64, <:Integer}},

	RFLU(A_gpu; nrhs=nrhs, symbolic=symbolic)
	RFLU(A_gpu; nrhs = nrhs, symbolic = symbolic)


		function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::LinearSolve.CUSOLVERRFFactorization; kwargs...)
		A = cache.A

	fact = RFLU(A_gpu; nrhs=nrhs, symbolic=alg.symbolic)
	fact = RFLU(A_gpu; nrhs = nrhs, symbolic = alg.symbolic)


		SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)

Uh oh!

Uh oh!

Conversation

ChrisRackauckas commented Jul 28, 2025

Summary

Motivation

Key Features

Implementation Details

Usage Example

Limitations

Testing

Uh oh!

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Jul 28, 2025

Choose a reason for hiding this comment

Uh oh!

oscardssmith commented Jul 28, 2025

Uh oh!

ChrisRackauckas commented Jul 28, 2025

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Aug 5, 2025

Choose a reason for hiding this comment

Uh oh!