Solver

SparseDirectSolver

Portable sparse direct solver (LLᵀ/LLᴴ, LDLᵀ/LDLᴴ, LDU) for GPUs on KernelAbstractions.jl, with the parameter names and phases of CUDSS.jl. See PLAN.md for the design.

AbstractDirectSolver{T, INT} <: LinearAlgebra.Factorization{T}

Supertype of the solvers of SparseDirectSolver.jl (≅ AbstractCudssSolver).

DirectSolver{T, INT, M, B, SY, NU, WS, RF, MS} <: AbstractDirectSolver{T, INT}

Sparse direct solver handle (≅ CudssSolver), PLAN §3.1–§3.2.

DirectSolver(A::CSR, structure::String, view::Char; index = A.index)
DirectSolver(rowptr, colval, nzval, structure::String, view::Char; index = 'O')
DirectSolver(A::SparseMatrixCSC, structure::String, view::Char; index = 'O')
DirectSolver(A::CuSparseMatrixCSR, structure, view; index = 'O')   # CUDA extension
DirectSolver(A::CuSparseMatrixCSC, structure, view; index = 'O')   # CUDA extension

structure is "G", "S", "H", "SPD" or "HPD" and view is 'L', 'U' or 'F' (the triangle of the matrix that is read; 'F' reads the lower triangle of a symmetric or Hermitian structure, and for "G" the view is ignored and every stored entry is read, as in cuDSS), as in CUDSS.jl; index ('O'/'Z') is the base of the row pointers and column indices. The CSR arrays are wrapped without copies and live on a KernelAbstractions backend (the solver's backend), except for a SparseMatrixCSC, which is converted to host CSR arrays once. A CSC matrix (CuSparseMatrixCSC, or a CSR with transposed = true, e.g. from csr_of_transpose) is read as the CSR of its transpose with the view flipped ('L' ↔ 'U'); for complex Hermitian matrices that transpose is the conjugate, and the solve phases conjugate the right-hand side and the solution (as for solve_mode, see solve_conjugated).

Implemented at this point: structures "SPD" (real T), "HPD", "S" (LDLᵀ; complex symmetric for complex T), "H" (LDLᴴ) and "G" (LU on the symmetric pattern of A + Aᵀ with row pivoting inside the fully-summed block of each front, P_r P A Pᵀ = L D U; any view, ignored), single matrices and uniform batches (see below), the phases "reordering", "symbolic_factorization", "analysis", "factorization", "refactorization", "solve", "solve_fwd_perm", "solve_fwd", "solve_diag", "solve_bwd", "solve_bwd_perm" and "solve_refinement" (execute!), and the Schur complement mode with "solve_fwd_schur" and "solve_bwd_schur" (below).

Uniform batch (PLAN §1.6, §3.5, ≅ CUDSS.jl's uniform batch): nbatch matrices with the pattern of rowptr/colval and values nzval, either a vector of nbatch · nnz entries (member after member) or an nnz × nbatch matrix; nbatch is deduced from the values (length(nzval) ÷ length(colval)), and "ubatch_size" must be 0 (deduced) or equal to it. One analysis serves every member; each member has its own factor, "info", "inertia", "npivots" and "pivot_stats" (vectors over the members). Right-hand sides and solutions hold nrhs columns per member, member after member: a strided vector of n · nrhs · nbatch entries, an n × (nrhs · nbatch) matrix, an n × nrhs × nbatch array or a MatrixDescriptor with nbatch. "ubatch_index" (0-based member, -1 = all) and "ubatch_mask" (one 0/1 flag per member) restrict the factorization and solve phases to a subset of the members: the other members' factors and solutions are left untouched. "pivot_sign" and the other options apply to every member; the 2×2 pivot pairs of "S"/"H" ("pivot_pairs") come from the first member's values.

Fields: A (the current CSR, re-pointed by update!), structure, view (as given), options (Options, set through setparam!), backend, nbatch, fresh_factorization (true until the first "factorization" after an analysis; the LinearAlgebra layer then switches cholesky! to "refactorization"), info (the "info" data parameter), the analysis (ordering, host host_symbolic, device symbolic), the numeric storage numeric, the solve workspace, the refinement storage refinement (allocated by the first solve that refines) and ir_steps (the "ir_n_steps" data parameter: refinement steps of the last solve, -1 before one) and schur (the SchurState of an analysis with schur_mode = 1, else nothing), host_matching (the host Matching of matching_alg ≠ "default", else nothing), analysis_colval (the column indices the analysis runs on: those of (Dr A Dc)[:, q] for "G" with matching, else host_colval) and matching (the device MatchingState, or nothing).

Matching and scaling (PLAN §1.3, M9, "matching_alg" = "algo1"–"algo6" before the analysis; compute_matching, from the first batch member's values): "G" factors M = (Dr A Dc)[:, q], whose diagonal holds the matched entries (job 5: |mᵢᵢ| = 1 ≥ |mᵢⱼ|), and the symmetric structures factor D A D (same inertia as A; for "S"/"H" with "algo5"/"algo6" the 2×2 pivot pairs come from the cycles of the matching, matching_pairs; jobs 1–4 have no scaling and keep the pairs of "pivot_pairs" without matching). The solve phases apply the permutation and the scalings ("solve_fwd_perm": Dr, "solve_bwd_perm": Dc Q, as in cuDSS) and refinement measures the residual of A itself; "diag" and the pivot statistics are those of the scaled matrix. Not with the Schur complement mode.

Schur complement mode (PLAN §3.6, ≅ CUDSS.jl's Schur complement): with "schur_mode" = 1 and "user_schur_indices" (n flags 0/1) set before the analysis, the flagged rows and columns are ordered last and the factorization stops before them: getparam(solver, "schur_shape") and getparam(solver, "schur_matrix") export S = A₂₂ − A₂₁ A₁₁⁻¹ A₁₂ (2 = the flagged block), and the solve phases are "solve_fwd_schur", "solve_diag" and "solve_bwd_schur" (see execute!). Single matrices only.

execute!(phase::String, solver::DirectSolver, X, B; asynchronous = true) -> nothing

Execute phase on solver (≅ cudss(phase, solver, x, b)), PLAN §1.2:

• "reordering": fill-reducing ordering on the host (reordering_alg, user_perm, …);
• "symbolic_factorization": supernodes, schedule, layout and device maps; allocates the numeric storage and the solve workspace (needs "reordering");
• "analysis": both;
• "factorization": numeric factorization with the current values of the matrix (needs the analysis); sets "info" (0 or, for "SPD"/"HPD", the original column of the first non-positive pivot; LDLᵀ/LDLᴴ and LU perturb tiny pivots instead and report 0, see "npivots", except for an exactly zero pivot with pivot_epsilon = 0, whose original column they report) and clears fresh_factorization;
• "refactorization": the same with the analysis and storage reused, "info" reset first (needs a previous "factorization");
• "solve": X = op(A)⁻¹ B (solve_mode 0: A, 1: Aᵀ, 2: Aᴴ), followed by ir_n_steps steps of iterative refinement; the sub-phases "solve_fwd_perm" (permute B into the workspace), "solve_fwd" (forward sweep), "solve_diag" (D⁻¹ of LDLᵀ/LDLᴴ, the identity for Cholesky), "solve_bwd" (backward sweep), "solve_bwd_perm" (inverse permutation into X) and "solve_refinement" (refine X in place, refine!) do the same in six calls with the same X, B (bitwise equal to "solve").

X and B are only read by the solve phases (pass anything, e.g. nothing, to the others): vectors, n × nrhs matrices, strided vectors or MatrixDescriptors (row-major when transposed) on the solver's backend; X === B is allowed in every phase but "solve_refinement" (it needs the original B; "solve" keeps a copy of B when they alias). A failed factorization does not throw (as in cuDSS): check getparam(solver, "info"); solving with it gives garbage.

Refinement (PLAN §2.5): each step computes the residual R = B - op(A) X with the solver's current matrix values (KA SpMV over the full pattern), and, when ir_tol > 0, stops once ‖Rₖ‖₂ ≤ ir_tol ‖Bₖ‖₂ for every right-hand side k (this costs one host synchronization per step; ir_tol = 0, the default, never stops early and never synchronizes), then adds the correction op(A)⁻¹ R. The steps performed are the data parameter "ir_n_steps" (getparam). With ir_mode = "fgmres" (needs using Krylov, else NotSupportedError) the refinement is FGMRES on the correction system, right-preconditioned by the factorization (fgmres_refine!): at most ir_n_steps iterations (one SpMV and one solve each, as a plain step), stopping once every relative residual is below ir_tol (ir_tol = 0: all ir_n_steps iterations); "ir_n_steps" then reports the FGMRES iterations.

The flag "user_host_interrupt" (a Threads.Atomic{Bool}, setparam!) is polled at the start of the analysis phases, between the launch groups of the factorization and between refinement steps: when set, the phase raises InterruptedError. An interrupted factorization leaves the solver analyzed (the next phase must be "factorization"); an interrupted refinement leaves the last completed iterate in X and reports the steps completed in "ir_n_steps".

Schur complement mode (schur_mode = 1 at the analysis, PLAN §3.6, as in the CUDSS.jl docs): "factorization"/"refactorization" factor A₁₁ (the rows and columns outside "user_schur_indices") and assemble the Schur complement S (getparam "schur_matrix"). The solve phases then read B and write X (X === B allowed), with intermediate vectors in factor order, where the last ns entries are the Schur block in increasing original order:

• "solve_fwd_schur" (permutation and forward sweep): the last ns entries of X are the condensed right-hand side bₛ = b₂ − A₂₁ A₁₁⁻¹ b₁;
• "solve_diag": D⁻¹ of LDLᵀ/LDLᴴ on the factored part (the identity on the Schur block, and for Cholesky and LU, whose diagonal is applied by "solve_bwd_schur" as in cuDSS);
• "solve_bwd_schur" (backward sweep and inverse permutation): with the solution x₂ of S x₂ = bₛ in the last ns entries of B, X is the solution of A x = b.

solve_mode 1/2 solve with Aᵀ/Aᴴ (condensed system Sᵀ/Sᴴ); the other solve phases raise NotSupportedError in Schur mode, and the Schur phases without it FactorizationError. Unknown phase strings raise ArgumentError. Executing a phase before the phases it depends on raises FactorizationError. With asynchronous = false the backend is synchronized before returning (KernelAbstractions.synchronize). Phase summaries are logged (see SparseDirectSolver.set_log_level!).

analyze!(solver, X = nothing, B = nothing; asynchronous = true) -> solver

execute!("analysis", solver, X, B; asynchronous).

factorize!(numeric, symbolic, nzval; impl = :auto, opts = Options()) -> info::Int

Structures "S"/"H": factorize_ldlt! with the pivoting options opts (impl selects its regime-C GEMMs); structure "G": factorize_lu! (impl is not used). Both return info (nonzero only for a zero pivot with pivot_epsilon = 0). Structures "SPD"/"HPD": multifrontal Cholesky P A Pᵀ = L Lᴴ on the device: nzval are the stored values of A (same pattern, view and index base as the analysis) on the backend of numeric, and symbolic has its maps on that backend (adapt(backend, symbolic, INT)). For every launch group of the schedule, in order: a regime-A group (budget class) is one factorize_subtrees! launch (one workgroup per subtree: its fronts assembled, factored and extend-added in local memory, panels and the root's contribution block written to global memory); a regime-B group (width class W ≤ 64) is one factorize_fronts_b! launch (fused assembly, F11 Cholesky in local memory, trsm, syrk/herk into the front's contribution block); a regime-C group is one launch each of zero_fronts!, scatter_A! and extend_add! over its fronts, then per front potrf_info!, trsm! and herk!/syrk! (into the workspace, then pack_add! into its packed contribution block) through the dense interface (impl is passed on: :auto, :generic, :vendor, :ka, resolved once per call; :auto means :ka when the analysis used factorization_alg = "algo1", no vendor calls). Then cholesky_stats! and one read of the reduced status.

Returns info = 0 on success, else the original column of the smallest non-positive pivot of the factor (as ref_factorize!); fronts above a failed one hold garbage. Uniform batch (numeric.nbatch > 1, PLAN §3.5): the active members of numeric (set_members!) are factorized by the same launches with the batch index as one more grid dimension (nzval holds the values of every member, member after member), a regime-C front by strided-batched potrf/trsm/gemm over each run of consecutive members (vendor potrf/trsm on the member pointers built at allocation; per member with impl = :generic); the return value is the first nonzero member info (see member_info! for all of them). opts.user_host_interrupt is polled before every launch group (a host read): when set, InterruptedError is raised and the factor is incomplete. Structure "SPD" (real) or "HPD". Nothing is allocated on the device; the panels are bitwise reproducible for a fixed impl and backend (regimes A and B are deterministic by construction).

factorize!(solver, X = nothing, B = nothing; asynchronous = true) -> solver

execute!("factorization", solver, X, B; asynchronous).

refactorize!(solver, X = nothing, B = nothing; asynchronous = true) -> solver

execute!("refactorization", solver, X, B; asynchronous).

solve!(solver, X, B; asynchronous = true) -> X

execute!("solve", solver, X, B; asynchronous).

update!(desc::MatrixDescriptor, x::AbstractArray) -> desc

Point desc at x without copying (≅ cudss_update(matrix, x)). x must hold nrows * ncols * nbatch entries of type T: either as a strided vector, or with the descriptor's shape ((nrows, ncols), (ncols, nrows) when transposed, plus a trailing nbatch dimension for 3-D arrays; a uniform batch also accepts the nrows × (ncols nbatch) matrix with the same memory, e.g. n × nbatch for MatrixDescriptor(T, n; nbatch), as CUDSS.jl). Anything else throws InvalidValueError.

update!(solver::DirectSolver, A)
update!(solver::DirectSolver, rowptr, colval, nzval)

Point solver at new matrix data without copying (≅ cudss_update): a CSR of the same type, a SparseMatrixCSC (converted on the host), a CuSparseMatrixCSR/CuSparseMatrixCSC (CUDA extension), or raw CSR arrays (with the index base and orientation of the solver's matrix). The sparsity pattern must be the one of the analysis (only the sizes are checked, as in cuDSS); new values are used by the next "factorization" or "refactorization". Raises InvalidValueError for a different size, number of stored entries, index base, orientation or array type.

update!(solver::DirectSolver, A::CuSparseMatrixCSR)
update!(solver::DirectSolver, A::CuSparseMatrixCSC)

Point solver at the arrays of A (≅ cudss_update(solver, A)), see update!.

update!(solver::DirectSolver, A::ROCSparseMatrixCSR)
update!(solver::DirectSolver, A::ROCSparseMatrixCSC)

Point solver at the arrays of A, see update!.

setparam!(solver::DirectSolver, name::String, value)

Set a configuration or user-input data parameter of solver (≅ cudss_set): everything setparam!(::Options, …) accepts, stored in solver.options and used by the next phase that reads it, plus the data parameter "info" (an integer, or for a uniform batch an integer for every member or a vector of nbatch integers; CUDSS.jl resets it before a refactorization). Uniform batch: "ubatch_size" (0 or the batch size), "ubatch_index" (0-based member, -1 = all) and "ubatch_mask" (nbatch flags 0/1, or nothing) select the members of the next phases; the index and the mask are checked against the batch size when a phase runs. "pivot_sign" (PLAN §1.7: a vector of n entries in (-1, 0, 1), host or device, or nothing) is checked against the size of the matrix; the next "factorization"/"refactorization" copies it to the device; it steers the perturbation of LDLᵀ/LDLᴴ only and is ignored for "G" (LU perturbs with the sign of the pivot) and the Cholesky structures. Setting "ir_n_steps" (the number of refinement steps requested) makes getparam report that value again until the next solve. "schur_matrix" (Schur complement mode, after "reordering"/"analysis") registers the destination that getparam(solver, "schur_matrix") fills (≅ cudss_set(solver, "schur_matrix", S) then cudss_get): a dense ns × ns matrix on the solver's backend (or a MatrixDescriptor holding one), a CSR (or a backend CSR matrix, e.g. CuSparseMatrixCSR) with ns + 1 row pointers and room for exactly the entries of the exported pattern, or a tuple (csr, view) with view 'L', 'U' or 'F' (the default) to export one triangle; nothing unregisters it. Computed data parameters ("lu_nnz", "diag", "perm_row", …) cannot be set (ArgumentError): read them with getparam or getparam!, which replaces cuDSS's set-buffer-then-get protocol.

getparam(solver::DirectSolver, name::String)

Value of a parameter of solver (≅ cudss_get). Configuration and user-input parameters come from solver.options (getparam(::Options, …)). The data parameters computed by the solver:

Uniform batch (nbatch > 1): "info", "npivots", "inertia" and "pivot_stats" are vectors with one entry per batch member, "diag" is the n · nbatch vector of the members' diagonals one after the other; the analysis outputs are shared by the members. A factorization restricted by "ubatch_index"/"ubatch_mask" updates the active members only: the others keep the "info", "npivots", "inertia" and "pivot_stats" of their last factorization. "perm_row" of an LU batch is batch member 1's row order, by design: cuDSS returns the vector outputs of a uniform batch ("diag") for member 1 only and does not implement "perm_row" at all (cuDSS 0.8). "diag" is the deliberate exception and returns every member.

The pivot statistics are reduced on the device (reduce_stats!) and copied to the host when read (one synchronization). The reordering permutation after "reordering" alone is the ordering algorithm's; "symbolic_factorization" composes it with the supernodal renumbering, and from then on both permutations are the one the factor uses. Data parameters of later tasks ("nd_partition_tree", "hybrid_device_memory_min") raise NotSupportedError; reading one before the phase that computes it raises FactorizationError, and the matching outputs without matching InvalidValueError.

getparam!(buffer, solver::DirectSolver, name::String) -> buffer

Write the vector-valued data parameter name ("perm_reorder_row", "perm_reorder_col", "perm_row", "perm_col", "diag", "memory_estimates", "user_perm") of getparam into buffer, a host or device vector of the right length (converted to its element type), without the C-style set-buffer-then-get protocol of CUDSS.jl. "schur_matrix": buffer is any destination setparam! accepts for it (dense matrix, CSR, (csr, view)), filled without registering it.

Options(; kwargs...)

Every configuration parameter of PLAN §1.3 and §1.7 plus the user-provided data parameters (user_perm, user_schur_indices, user_nd_partition_tree, user_host_interrupt, ubatch_mask, pivot_sign), in typed fields with the defaults below. Keyword arguments are applied through setparam!, so Options(ir_n_steps = 2) validates like setparam!(opts, "ir_n_steps", 2).

Vectors are stored as host copies; user_host_interrupt is stored by reference because it is polled while a phase runs.

The last five rows are analysis tuning knobs of the schedule (PLAN §2.3 step 5, TASKS T07). They are not cuDSS parameter strings, so setparam! does not know them; set them with Options(; regime_c_width = 128) (validated) or by assigning the field.

setparam!(opts::Options, name::String, value)

Set the configuration or user-input data parameter name (a CUDSS.jl name from CONFIG_PARAMETERS/DATA_PARAMETERS or a name from EXTRA_PARAMETERS) after validating value.

Accepted values:

• algorithms ("reordering_alg", "factorization_alg", "solve_alg", "matching_alg", "pivot_epsilon_alg"): "default", "algoN" or the integer N;
• "pivot_type": 'A', 'N', 'D', 'L', 'B' ('C'/'R' raise NotSupportedError); LDLᵀ/LDLᴴ: 'A', 'B', 'L' Bunch–Kaufman, 'D' 1×1 only, 'N' none; LU ("G", lu_pivoting): 'A', 'B', 'L' in-block threshold row pivoting, 'N', 'D' diagonal pivots;
• integer and flag parameters: an Integer in the documented range;
• "ir_tol", "pivot_threshold", "pivot_epsilon": a finite real ≥ 0 (nothing resets "pivot_epsilon" to default_pivot_epsilon);
• "ir_mode": "ir" or "fgmres"; "schedule": "auto", "subtree+level", "syncfree";
• "pivot_pairs" (beyond cuDSS): "default" (for "S"/"H", a row with a zero or negligible diagonal whose pivot would be structurally zero is ordered together with a 2×2 pivot partner, see zero_pivot_pairs!), "all" (every such row, see pivot_pairs; about 2× nnz(L) on KKT systems) or "none"; ignored for the other structures, with user_perm and with the natural ordering; with matching_alg = "algo5"/"algo6" the pairs come from the cycles of the scaled symmetric matching instead (matching_pairs). Pairs are decided from the values present at "analysis" (the first batch member): an all-zero nzval gives no pairs and an undefined one arbitrary pairs, so run "analysis" after the first assembly of the matrix (MadNLP: after the first KKT assembly); "all" reads the values as well (its candidates and partners come from them);
• "pivot_threshold" also shapes the ordering with pivot_pairs = "default": partners are accepted at analysis down to pivot_threshold · maxₖ |aᵢₖ| (zero_pivot_pairs!); a value set after "analysis" changes the in-front pivot test but not the pairs already chosen;
• "pivot_pair_tolerance" (beyond cuDSS): a finite real ≥ 0, the relative diagonal size τ below which a row is a 2×2 candidate (|aᵢᵢ| ≤ τ maxⱼ≠ᵢ |aᵢⱼ|, on the scaled matrix with matching; default PIVOT_PAIR_TOLERANCE);
• "factor_precision": Float32, Float64 or nothing;
• "amalgamation": a NamedTuple with any of max_width, zero_fraction, min_width;
• "user_perm", "user_nd_partition_tree": an integer vector (host or device), or nothing;
• "user_schur_indices", "ubatch_mask": a vector of 0/1 flags, or nothing;
• "pivot_sign": a vector with entries in (-1, 0, 1), or nothing;
• "user_host_interrupt": a Threads.Atomic{Bool}, or nothing.

Unknown names raise ArgumentError; values of the wrong type or range raise InvalidValueError; "device_count" ≠ 1, "device_indices", "comm_device", "comm_host" and global pivoting raise NotSupportedError. Parameters computed by a solver ("lu_nnz", "inertia", …) cannot be set (ArgumentError). The hybrid-mode parameters are stored but have no effect yet; a non-default value warns once. reordering_alg = "algo1"/"algo2" (COLAMD-based orderings with global pivoting) is accepted and warns once that the symmetric-pattern path is used.

getparam(opts::Options, name::String)

Value of the configuration or user-input data parameter name, in the form setparam! accepts: algorithm and mode strings ("algo3", "fgmres"), Char for "pivot_type", Int/Float64 for numbers, a copy of stored vectors, nothing for unset optional values ("pivot_epsilon", "factor_precision", user vectors). "device_count" is always 1.

Unknown names raise ArgumentError, as do data parameters computed by a solver (query those on the solver). "device_indices", "comm_device" and "comm_host" raise NotSupportedError.

CONFIG_PARAMETERS

The configuration parameter names of CUDSS.jl (CUDSS_CONFIG_PARAMETERS), verbatim.

DATA_PARAMETERS

The data parameter names of CUDSS.jl (CUDSS_DATA_PARAMETERS), verbatim.

CUDSS08_DATA_PARAMETERS

cuDSS 0.8 data parameters that CUDSS.jl does not list in DATA_PARAMETERS (PLAN §1.4): "ir_n_steps" (refinement steps performed; the configuration parameter of the same name is the number requested), "ubatch_mask" and "flops".

EXTRA_PARAMETERS

Parameters beyond cuDSS (PLAN §1.7): the data parameters "pivot_sign" (input) and "pivot_stats" (output), and the configuration parameters "ir_mode", "factor_precision", "amalgamation", "schedule", "pivot_pairs" and "pivot_pair_tolerance".

default_pivot_epsilon(T) -> Float64

Default pivot_epsilon for element type T, as in cuDSS: 1e-5 for Float32/ComplexF32, 1e-13 for Float64/ComplexF64.

SparseDirectSolver.set_log_level!(level) -> Int

Set the verbosity of the solver's log messages: 0/"none" (default; the messages are @debug records), 1/"info" (phase summaries: analysis, factorization, refinement steps performed, as @info), 2/"debug" (also the relative residual of every refinement step). The initial value comes from the environment variable SDS_LOG_LEVEL. Returns the previous level.

SparseDirectSolverError <: Exception

Abstract supertype of the exceptions raised by SparseDirectSolver.jl.

NotSupportedError(msg)

A parameter, value or feature that SparseDirectSolver.jl does not provide (for example "device_count" > 1, or global pivoting pivot_type = 'C').

InvalidValueError(msg)

A value of the wrong type or outside the admissible range for a parameter, structure, view, index base or phase.

FactorizationError(info, msg = "")

A factorization that could not be completed, or a phase called out of order. info follows the "info" data parameter: k > 0 is the first failed pivot (1-based, original numbering); for uniform batches it is a vector with one entry per member.

InterruptedError(msg = "interrupted through "user_host_interrupt"")

Raised when the "user_host_interrupt" flag is set while a phase runs. The solver stays usable: the phase can be executed again once the flag is cleared.