Refactor CudaOffloadFactorization into LU and QR variants

docs/src/release_notes.md

@@ -1,5 +1,12 @@
 # Release Notes
 
+## Upcoming Changes
+
+  - `CudaOffloadFactorization` has been split into two algorithms:
+    - `CudaOffloadLUFactorization` - Uses LU factorization for better performance
+    - `CudaOffloadQRFactorization` - Uses QR factorization for better numerical stability
+  - `CudaOffloadFactorization` is now deprecated and will show a warning suggesting to use one of the new algorithms
+
 ## v2.0
 
   - `LinearCache` changed from immutable to mutable. With this, the out of place interfaces like

docs/src/solvers/solvers.md

@@ -23,12 +23,14 @@ use your base system BLAS which can be fast or slow depending on the hardware co
 
 For very large dense factorizations, offloading to the GPU can be preferred. Metal.jl can be used
 on Mac hardware to offload, and has a cutoff point of being faster at around size 20,000 x 20,000
-matrices (and only supports Float32). `CudaOffloadFactorization` can be more efficient at a
-much smaller cutoff, possibly around size 1,000 x 1,000 matrices, though this is highly dependent
-on the chosen GPU hardware. `CudaOffloadFactorization` requires a CUDA-compatible NVIDIA GPU.
+matrices (and only supports Float32). `CudaOffloadLUFactorization` and `CudaOffloadQRFactorization` 
+can be more efficient at a much smaller cutoff, possibly around size 1,000 x 1,000 matrices, though 
+this is highly dependent on the chosen GPU hardware. These algorithms require a CUDA-compatible NVIDIA GPU.
 CUDA offload supports Float64 but most consumer GPU hardware will be much faster on Float32
 (many are >32x faster for Float32 operations than Float64 operations) and thus for most hardware
-this is only recommended for Float32 matrices.
+this is only recommended for Float32 matrices. Choose `CudaOffloadLUFactorization` for better 
+performance on well-conditioned problems, or `CudaOffloadQRFactorization` for better numerical 
+stability on ill-conditioned problems.
 
 !!! note
 
@@ -232,9 +234,11 @@ The following are non-standard GPU factorization routines.
 
 !!! note
 
-    Using this solver requires adding the package CUDA.jl, i.e. `using CUDA`
+    Using these solvers requires adding the package CUDA.jl, i.e. `using CUDA`
 
 ```@docs
+CudaOffloadLUFactorization
+CudaOffloadQRFactorization
 CudaOffloadFactorization
 ```
 

docs/src/tutorials/gpu.md

@@ -40,10 +40,17 @@ This computation can be moved to the GPU by the following:
 
 ```julia
 using CUDA # Add the GPU library
-sol = LS.solve(prob, LS.CudaOffloadFactorization())
+sol = LS.solve(prob, LS.CudaOffloadLUFactorization())
 sol.u
 ```
 
+LinearSolve.jl provides two GPU offloading algorithms:
+- `CudaOffloadLUFactorization()` - Uses LU factorization (generally faster for well-conditioned problems)
+- `CudaOffloadQRFactorization()` - Uses QR factorization (more stable for ill-conditioned problems)
+
+!!! warning
+    The old `CudaOffloadFactorization()` is deprecated. Use `CudaOffloadLUFactorization()` or `CudaOffloadQRFactorization()` instead.
+
 ## GPUArray Interface
 
 For more manual control over the factorization setup, you can use the

examples/cuda_offload_example.jl

@@ -0,0 +1,97 @@
+"""
+Example demonstrating the new CudaOffloadLUFactorization and CudaOffloadQRFactorization algorithms.
+
+This example shows how to use the new GPU offloading algorithms for solving linear systems
+with different numerical properties.
+"""
+
+using LinearSolve
+using LinearAlgebra
+using Random
+
+# Set random seed for reproducibility
+Random.seed!(123)
+
+println("CUDA Offload Factorization Examples")
+println("=" ^ 40)
+
+# Create a well-conditioned problem
+println("\n1. Well-conditioned problem (condition number ≈ 10)")
+A_good = rand(100, 100)
+A_good = A_good + 10I  # Make it well-conditioned
+b_good = rand(100)
+prob_good = LinearProblem(A_good, b_good)
+
+println("   Matrix size: $(size(A_good))")
+println("   Condition number: $(cond(A_good))")
+
+# Try to use CUDA if available
+try
+    using CUDA
+
+    # Solve with LU (faster for well-conditioned)
+    println("\n   Solving with CudaOffloadLUFactorization...")
+    sol_lu = solve(prob_good, CudaOffloadLUFactorization())
+    println("   Solution norm: $(norm(sol_lu.u))")
+    println("   Residual norm: $(norm(A_good * sol_lu.u - b_good))")
+
+    # Solve with QR (more stable)
+    println("\n   Solving with CudaOffloadQRFactorization...")
+    sol_qr = solve(prob_good, CudaOffloadQRFactorization())
+    println("   Solution norm: $(norm(sol_qr.u))")
+    println("   Residual norm: $(norm(A_good * sol_qr.u - b_good))")
+
+catch e
+    println("\n   Note: CUDA.jl is not loaded. To use GPU offloading:")
+    println("   1. Install CUDA.jl: using Pkg; Pkg.add(\"CUDA\")")
+    println("   2. Add 'using CUDA' before running this example")
+    println("   3. Ensure you have a CUDA-compatible NVIDIA GPU")
+end
+
+# Create an ill-conditioned problem
+println("\n2. Ill-conditioned problem (condition number ≈ 1e6)")
+A_bad = rand(50, 50)
+# Make it ill-conditioned
+U, S, V = svd(A_bad)
+S[end] = S[1] / 1e6  # Create large condition number
+A_bad = U * Diagonal(S) * V'
+b_bad = rand(50)
+prob_bad = LinearProblem(A_bad, b_bad)
+
+println("   Matrix size: $(size(A_bad))")
+println("   Condition number: $(cond(A_bad))")
+
+try
+    using CUDA
+
+    # For ill-conditioned problems, QR is typically more stable
+    println("\n   Solving with CudaOffloadQRFactorization (recommended for ill-conditioned)...")
+    sol_qr_bad = solve(prob_bad, CudaOffloadQRFactorization())
+    println("   Solution norm: $(norm(sol_qr_bad.u))")
+    println("   Residual norm: $(norm(A_bad * sol_qr_bad.u - b_bad))")
+
+    println("\n   Solving with CudaOffloadLUFactorization (may be less stable)...")
+    sol_lu_bad = solve(prob_bad, CudaOffloadLUFactorization())
+    println("   Solution norm: $(norm(sol_lu_bad.u))")
+    println("   Residual norm: $(norm(A_bad * sol_lu_bad.u - b_bad))")
+
+catch e
+    println("\n   Skipping GPU tests (CUDA not available)")
+end
+
+# Demonstrate the deprecation warning
+println("\n3. Testing deprecated CudaOffloadFactorization")
+try
+    using CUDA
+    println("   Creating deprecated CudaOffloadFactorization...")
+    alg = CudaOffloadFactorization()  # This will show a deprecation warning
+    println("   The deprecated algorithm still works but shows a warning above")
+catch e
+    println("   Skipping deprecation test (CUDA not available)")
+end
+
+println("\n" * "=" ^ 40)
+println("Summary:")
+println("- Use CudaOffloadLUFactorization for well-conditioned problems (faster)")
+println("- Use CudaOffloadQRFactorization for ill-conditioned problems (more stable)")
+println("- The old CudaOffloadFactorization is deprecated")

ext/LinearSolveCUDAExt.jl

@@ -5,7 +5,7 @@ using LinearSolve: LinearSolve, is_cusparse, defaultalg, cudss_loaded, DefaultLi
                    DefaultAlgorithmChoice, ALREADY_WARNED_CUDSS, LinearCache,
                    needs_concrete_A,
                    error_no_cudss_lu, init_cacheval, OperatorAssumptions,
-                   CudaOffloadFactorization,
+                   CudaOffloadFactorization, CudaOffloadLUFactorization, CudaOffloadQRFactorization,
                    SparspakFactorization, KLUFactorization, UMFPACKFactorization
 using LinearSolve.LinearAlgebra, LinearSolve.SciMLBase, LinearSolve.ArrayInterface
 using SciMLBase: AbstractSciMLOperator
@@ -35,6 +35,43 @@ function LinearSolve.error_no_cudss_lu(A::CUDA.CUSPARSE.CuSparseMatrixCSR)
     nothing
 end
 
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadLUFactorization;
+        kwargs...)
+    if cache.isfresh
+        fact = lu(CUDA.CuArray(cache.A))
+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+    y = Array(ldiv!(CUDA.CuArray(cache.u), cache.cacheval, CUDA.CuArray(cache.b)))
+    cache.u .= y
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end
+
+function LinearSolve.init_cacheval(alg::CudaOffloadLUFactorization, A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    ArrayInterface.lu_instance(A)
+end
+
+function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadQRFactorization;
+        kwargs...)
+    if cache.isfresh
+        fact = qr(CUDA.CuArray(cache.A))
+        cache.cacheval = fact
+        cache.isfresh = false
+    end
+    y = Array(ldiv!(CUDA.CuArray(cache.u), cache.cacheval, CUDA.CuArray(cache.b)))
+    cache.u .= y
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
+end
+
+function LinearSolve.init_cacheval(alg::CudaOffloadQRFactorization, A, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    qr(CUDA.CuArray(A))
+end
+
+# Keep the deprecated CudaOffloadFactorization working by forwarding to QR
 function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::CudaOffloadFactorization;
         kwargs...)
     if cache.isfresh

lib/LinearSolveAutotune/src/algorithms.jl

@@ -84,10 +84,10 @@ function get_gpu_algorithms(; skip_missing_algs::Bool = false)
     # CUDA algorithms
     if is_cuda_available()
         try
-            push!(gpu_algs, CudaOffloadFactorization())
-            push!(gpu_names, "CudaOffloadFactorization")
+            push!(gpu_algs, CudaOffloadLUFactorization())
+            push!(gpu_names, "CudaOffloadLUFactorization")
         catch e
-            msg = "CUDA hardware detected but CudaOffloadFactorization not available: $e. Load CUDA.jl package."
+            msg = "CUDA hardware detected but CudaOffloadLUFactorization not available: $e. Load CUDA.jl package."
             if skip_missing_algs
                 @warn msg
             else

src/LinearSolve.jl

@@ -255,6 +255,8 @@ export SimpleGMRES
 
 export HYPREAlgorithm
 export CudaOffloadFactorization
+export CudaOffloadLUFactorization
+export CudaOffloadQRFactorization
 export MKLPardisoFactorize, MKLPardisoIterate
 export PanuaPardisoFactorize, PanuaPardisoIterate
 export PardisoJL

src/extension_algs.jl

@@ -61,23 +61,70 @@ struct HYPREAlgorithm <: SciMLLinearSolveAlgorithm
     end
 end
 
+# Debug: About to define CudaOffloadLUFactorization
+"""
+`CudaOffloadLUFactorization()`
+
+An offloading technique used to GPU-accelerate CPU-based computations using LU factorization.
+Requires a sufficiently large `A` to overcome the data transfer costs.
+
+!!! note
+
+    Using this solver requires adding the package CUDA.jl, i.e. `using CUDA`
+"""
+struct CudaOffloadLUFactorization <: AbstractFactorization
+    function CudaOffloadLUFactorization()
+        ext = Base.get_extension(@__MODULE__, :LinearSolveCUDAExt)
+        if ext === nothing
+            error("CudaOffloadLUFactorization requires that CUDA is loaded, i.e. `using CUDA`")
+        else
+            return new()
+        end
+    end
+end
+
+"""
+`CudaOffloadQRFactorization()`
+
+An offloading technique used to GPU-accelerate CPU-based computations using QR factorization.
+Requires a sufficiently large `A` to overcome the data transfer costs.
+
+!!! note
+
+    Using this solver requires adding the package CUDA.jl, i.e. `using CUDA`
+"""
+struct CudaOffloadQRFactorization <: AbstractFactorization
+    function CudaOffloadQRFactorization()
+        ext = Base.get_extension(@__MODULE__, :LinearSolveCUDAExt)
+        if ext === nothing
+            error("CudaOffloadQRFactorization requires that CUDA is loaded, i.e. `using CUDA`")
+        else
+            return new()
+        end
+    end
+end
+
 """
 `CudaOffloadFactorization()`
 
+!!! warning
+    This algorithm is deprecated. Use `CudaOffloadLUFactorization` or `CudaOffloadQRFactorization()` instead.
+
 An offloading technique used to GPU-accelerate CPU-based computations.
 Requires a sufficiently large `A` to overcome the data transfer costs.
 
 !!! note
 
     Using this solver requires adding the package CUDA.jl, i.e. `using CUDA`
 """
-struct CudaOffloadFactorization <: LinearSolve.AbstractFactorization
+struct CudaOffloadFactorization <: AbstractFactorization
     function CudaOffloadFactorization()
+        Base.depwarn("`CudaOffloadFactorization` is deprecated, use `CudaOffloadLUFactorization` or `CudaOffloadQRFactorization` instead.", :CudaOffloadFactorization)
         ext = Base.get_extension(@__MODULE__, :LinearSolveCUDAExt)
         if ext === nothing
             error("CudaOffloadFactorization requires that CUDA is loaded, i.e. `using CUDA`")
         else
-            return new{}()
+            return new()
         end
     end
 end

test/gpu/cuda.jl

@@ -45,7 +45,7 @@ function test_interface(alg, prob1, prob2)
     return
 end
 
-@testset "$alg" for alg in (CudaOffloadFactorization(), NormalCholeskyFactorization())
+@testset "$alg" for alg in (CudaOffloadLUFactorization(), CudaOffloadQRFactorization(), NormalCholeskyFactorization())
     test_interface(alg, prob1, prob2)
 end
 

test/resolve.jl

@@ -11,6 +11,8 @@ for alg in vcat(InteractiveUtils.subtypes(AbstractDenseFactorization),
     if !(alg in [
            DiagonalFactorization,
            CudaOffloadFactorization,
+           CudaOffloadLUFactorization,
+           CudaOffloadQRFactorization,
            CUSOLVERRFFactorization,
            AppleAccelerateLUFactorization,
            MetalLUFactorization