Linear Solver API#

cunls/linear_solver hosts linear-system abstractions (block-Jacobi PCG, cuDSS integration, dense pivoted LDLT, dense Cholesky, and dense QR solvers) behind a common interface. Every backend accepts scalar CSR; all but cuDSS also accept block BSR directly, which is what SupportsBlockStorage() reports.

This page documents the C++ API only: pycunls does not expose the linear solver classes. From Python, the backend is selected through the sparse_linear_solver_type field of pycunls.MinimizerOptions (see pycunls.SparseLinearSolverType).

SparseLinearSolverType#

Enum in cunls/linear_solver/sparse_linear_solver.h:

  • BlockSparsePCG — iterative block-Jacobi preconditioned conjugate gradient solver. Default backend for Gauss-Newton and Levenberg-Marquardt: outperforms cuDSS on most SBA / PGO workloads (see profile/PCG_RESULTS.md).

  • cuDSS — NVIDIA cuDSS sparse direct solver. Pick when each Solve sees a tiny system and PCG’s per-iter kernel-launch overhead dominates.

  • DenseLDLT — densifies the coefficient matrix and solves with a custom CUDA pivoted LDLT kernel.

  • DenseCholesky — densifies the coefficient matrix and solves with cuSOLVER Cholesky (requires SPD).

  • DenseQR — densifies the coefficient matrix and solves with cuSOLVER QR.

SparseLinearSolverConfig#

Struct in sparse_linear_solver.h. Only the member corresponding to the chosen SparseLinearSolverType is used:

  • block_sparse_pcg_options - [in] BlockSparsePCG-specific knobs; ignored when a different backend is selected.

  • cudss_solver_options - [in] cuDSS-specific options; ignored when a different backend is selected.

  • Dense backends take no extra configuration.

SparseLinearSolver#

Abstract base (cunls/linear_solver/sparse_linear_solver_base.h).

bool Initialize(
cudaStream_t stream,
const Problem &problem,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Performs setup work (at minimum symbolic analysis; some modes also perform an initial factorization) for the given sparsity pattern. Must be called before Solve() and re-called whenever the matrix structure changes.

Solvers may inspect problem to specialize their setup (e.g. BlockSparsePCGSolver reads each state batch’s TangentSize to build the block-Jacobi layout). Solvers that do not need it (cuDSS, dense backends) simply ignore it; for callers that operate on a raw matrix without a cuNLS problem context, pass a default-constructed Problem.

Both rhs and result must be pre-allocated to the same number of elements as the number of rows in spd_matrix; the solver does not resize them.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • problem – [in] Originating problem; used by problem-aware backends to derive per-solver structure.

  • spd_matrix – [in] Symmetric matrix in CSR format (backend-specific definiteness requirements may apply).

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

bool Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Performs factorization and solve phases to compute the solution x of A x = b. Both rhs and result must be pre-allocated to the same number of elements as the number of rows in spd_matrix; the solver does not resize them. Call Initialize() first for the current sparsity pattern.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • spd_matrix – [in] Symmetric matrix in CSR format (backend-specific definiteness requirements may apply).

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

void DisableSafetyChecks()#

Advisory request to skip post-factorization safety checks; backends with no such phase (cuDSS, which reports through its own status, and PCG, which has no factorization) ignore it. By default (safety checks enabled), dense solvers validate every factorization and solve step: Cholesky checks cuSOLVER devInfo after potrf and potrs; QR inspects the diagonal of R for rank deficiency; LDLT performs in-kernel pivot and diagonal checks. Failures cause Solve() to return false with a diagnostic via LogError(). Calling this method skips all of the above (no device-to-host memcpy, no stream synchronization, no in-kernel validation), which can noticeably reduce per-iteration latency for small systems but may produce silently incorrect results for singular or ill-conditioned matrices. Normally called by the minimizer when MinimizerOptions::disable_safety_checks is true.

BlockSparsePCGOptions#

Struct in cunls/linear_solver/block_sparse_pcg_solver.h. Convergence and preconditioner-layout knobs for BlockSparsePCGSolver.

  • block_size - [in] Uniform diagonal tile size; used only when block_layout is empty. Must divide the matrix dimension. block_size = 1 reduces the preconditioner to scalar Jacobi. Default: 6.

  • block_layout - [in] Optional std::vector<std::pair<int, int>> describing a heterogeneous block-diagonal layout: each (count_i, size_i) element contributes count_i * size_i rows with size_i-square diagonal tiles. When empty (default), Initialize derives the layout automatically from the Problem’s state batches (segment per non-empty batch with size = TangentSize(), count = NumActiveStates() - NumConstStates()).

  • max_iterations - [in] PCG iteration cap. Default: 200.

  • relative_tolerance - [in] Stop when ||r_k|| <= relative_tolerance * ||b||. Default: 1e-3.

  • absolute_tolerance - [in] Stop also when ||r_k|| <= absolute_tolerance. Default: 1e-30.

  • pivot_floor - [in] Floor on |D_{kk}| during the per-tile LDLT pivoting; keeps the preconditioner invertible on near-singular diagonal tiles. Default: 1e-12.

  • check_period - [in] Number of PCG iterations between host-side convergence polls. Higher values reduce host syncs but may do up to check_period - 1 extra iterations after convergence. Default: 4.

BlockSparsePCGSolver#

Concrete implementation in block_sparse_pcg_solver.h. Solves H x = b for symmetric positive-(semi-)definite H with the standard preconditioned conjugate gradient recurrence (Saad, Iterative Methods for Sparse Linear Systems, §9.2) and a block-Jacobi preconditioner formed from the dense diagonal tiles of H. Each tile is factored independently with an in-shared-memory LDLT. SpMV is delegated to cuSPARSE.

All scalar quantities (alpha, beta, <p, q>, <r, r>, <r, z>) live on the device for the whole inner loop. Only the residual norm is copied to the host, and only once every check_period iterations.

BlockSparsePCGSolver(
BlockSparsePCGOptions options = BlockSparsePCGOptions()
)#
Parameters:

options – [in] Convergence / preconditioner-layout knobs. See BlockSparsePCGOptions.

Returns:

Constructor has no return value.

bool BlockSparsePCGSolver::Initialize(
cudaStream_t stream,
const Problem &problem,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Builds the cuSPARSE SpMV plan, allocates the PCG scratch vectors, and (when options.block_layout is empty) derives the block-Jacobi layout from problem.GetStateBatches().

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • problem – [in] Source of the block-Jacobi layout when options.block_layout is empty. Pass a default-constructed Problem when calling on a raw matrix.

  • spd_matrix – [in] Coefficient matrix H in CSR format.

  • rhs – [in] Right-hand-side vector b.

  • result – [out] Solution vector x.

Returns:

true on success; false on dimension mismatch or invalid layout.

bool BlockSparsePCGSolver::Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Rebuilds the block-Jacobi factor from the current values of H and runs up to options.max_iterations PCG iterations with the zero initial guess x_0 = 0.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • spd_matrix – [in] Coefficient matrix H (same structure as in Initialize).

  • rhs – [in] Right-hand-side vector b.

  • result – [out] Solution vector x.

Returns:

true on success.

int BlockSparsePCGSolver::LastIterations() const#
Returns:

Number of PCG iterations consumed by the most recent Solve call. Useful for convergence diagnostics.

cuDSSLinearSolverMode#

Enum in cunls/linear_solver/cudss_sparse_linear_solver.h:

  • SlowInitFastSolve — Initialize performs symbolic analysis only; Solve performs a full factorization followed by the solve phase. Best when the matrix values change frequently relative to its structure.

  • FastInitSlowSolve — Initialize performs symbolic analysis and an initial factorization; Solve performs refactorization followed by the solve phase. Best when the sparsity pattern is stable and fast initialization is more important than per-solve speed.

cuDSSLinearSolverOptions#

Struct fields:

  • mode - [in] Trade-off between setup and repeated solve speed (see cuDSSLinearSolverMode).

  • nthreads - [in] Host thread count for cuDSS host-side stages.

  • threading_lib_path - [in] Optional threading runtime path for multi-threaded host processing. Empty string disables multi-threading.

cuDSSLinearSolver#

Concrete implementation in cudss_sparse_linear_solver.h.

cuDSSLinearSolver(
cuDSSLinearSolverOptions options = cuDSSLinearSolverOptions()
)#
Parameters:

options – [in] cuDSS mode/threading configuration.

Returns:

Constructor has no return value.

bool cuDSSLinearSolver::Initialize(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

For SlowInitFastSolve: runs the cuDSS symbolic analysis phase only. For FastInitSlowSolve: runs symbolic analysis followed by an initial factorization so that subsequent cuDSSLinearSolver::Solve() calls can use cheaper refactorization.

Both rhs and result must be pre-allocated to the number of rows in spd_matrix; the solver does not resize them. The internal cuDSS handle is lazily initialized on the first call using the provided stream.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations. Also used to lazily initialize the internal cuDSS handle on first call.

  • spd_matrix – [in] SPD matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

DenseLDLTSolver#

Concrete implementation in dense_linear_solver.h.

bool DenseLDLTSolver::Initialize(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Validates dimensions and allocates reusable dense buffers used by the dense pivoted LDLT pipeline.

Parameters:

  • stream – [in] CUDA stream (unused; buffers are allocated synchronously).

  • spd_matrix – [in] Symmetric matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

bool DenseLDLTSolver::Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Converts CSR to dense row-major storage, computes a CUDA pivoted LDLT factorization P^T A P = L D L^T, then solves in the permuted system. Handles symmetric matrices including indefinite ones (not limited to SPD). Kernel success/failure (singular pivot, zero diagonal) is propagated back to the host via a single async device-to-pinned-host copy after both kernels complete.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • spd_matrix – [in] Symmetric matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false on dimension mismatch, singular pivot during factorization, or zero diagonal during the solve phase.

DenseCholeskySolver#

Concrete implementation in dense_cholesky_solver.h. Uses the cuSOLVER cusolverDnSpotrf / cusolverDnSpotrs routines (Cholesky factorization A = L L^T). Requires the input matrix to be symmetric positive-definite.

bool DenseCholeskySolver::Initialize(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Validates dimensions and pre-allocates internal dense buffers and the cuSOLVER workspace for the given matrix size.

Parameters:

  • stream – [in] CUDA stream used to query the cuSOLVER workspace size.

  • spd_matrix – [in] SPD matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

bool DenseCholeskySolver::Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Converts CSR to dense, performs Cholesky factorization via cusolverDnSpotrf, then solves via cusolverDnSpotrs. When safety checks are enabled, inspects cuSOLVER devInfo after both potrf and potrs: returns false if the matrix is not positive-definite (devInfo > 0 from potrf) or if potrs reports an invalid parameter (devInfo < 0).

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • spd_matrix – [in] SPD matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false on dimension mismatch, non-positive-definite matrix, or invalid parameter from potrs.

DenseQRSolver#

Concrete implementation in dense_qr_solver.h. Uses the cuSOLVER cusolverDnSgeqrf / cusolverDnSormqr routines for QR factorization (A = Q R) followed by a cuBLAS cublasStrsm upper-triangular solve. Works for any non-singular square matrix (not limited to SPD).

bool DenseQRSolver::Initialize(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Validates dimensions and pre-allocates internal dense buffers, Householder tau vector, and the cuSOLVER workspace for the given matrix size.

Parameters:

  • stream – [in] CUDA stream used to query the cuSOLVER workspace size.

  • spd_matrix – [in] Matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

bool DenseQRSolver::Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

Converts CSR to dense, performs QR factorization via cusolverDnSgeqrf, applies Q^T to the RHS via cusolverDnSormqr, then solves the upper-triangular system R x = Q^T b via cublasStrsm.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations.

  • spd_matrix – [in] Matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false on dimension mismatch or singular matrix (cuSOLVER devInfo != 0).

bool cuDSSLinearSolver::Solve(
cudaStream_t stream,
const CSRSparseMatrix &spd_matrix,
const dvector<float> &rhs,
dvector<float> &result
)#

For SlowInitFastSolve: performs full factorization then solve. For FastInitSlowSolve: performs refactorization (reusing the symbolic analysis from Initialize) then solve.

Both rhs and result must be pre-allocated to the number of rows in spd_matrix; the solver does not resize them and returns false on any dimension mismatch.

Parameters:

  • stream – [in] CUDA stream for asynchronous GPU operations. Also used to lazily initialize the internal cuDSS handle if needed.

  • spd_matrix – [in] SPD matrix in CSR format.

  • rhs – [in] Right-hand-side vector b (size must equal matrix rows).

  • result – [out] Solution vector x (size must equal matrix rows).

Returns:

true on success; false if a dimension mismatch is detected.

Factory#

Function in sparse_linear_solver.h:

SparseLinearSolverPtr CreateSparseLinearSolver(
SparseLinearSolverType type,
const SparseLinearSolverConfig &config
)#
Parameters:

  • type – [in] Backend type to instantiate.

  • config – [in] Backend-specific configuration blob.

Returns:

[out] Heap-allocated solver instance.