GitHub

KKT solver

KKTSolver is the front end for the linear systems of an interior-point method: one real symmetric matrix with a fixed sparsity pattern, factorized again with new values at every iteration, solved for one right-hand side at a time, and the inertia read after each factorization. It chooses the pipeline (the fused pipeline of FusedSolver where the backend runs it, else the stock DirectSolver), sets the parameters of the chosen factorization and drives the phases. The constructor runs the analysis and the first factorization (with deferred_analysis = true the first factorize! does, on the values present then); afterwards the loop calls only factorize!, solve! and inertia. MadNLP's SparseDirectSolver binding is built on it (MadNLP).

using SparseDirectSolver, SparseArrays, CUDA, CUDA.CUSPARSE

# K: a symmetric quasi-definite KKT matrix on the host, its lower triangle on the GPU
Kd = CuSparseMatrixCSC(tril(K))
S = KKTSolver(Kd, 'L'; factorization = :ldl)   # analysis and first factorization
SparseDirectSolver.pipeline(S)                 # :fused on CUDA, AMDGPU and oneAPI

Kd.nzVal .= new_values                         # same pattern, the next iteration
factorize!(S)
getparam(S, "info") == 0 || error("factorization failed")
npos, nzero, nneg = SparseDirectSolver.inertia(S)
xb = CuVector(b)
solve!(S, xb)                                  # xb holds b on entry, x on return

On the CPU backend and on Metal the same calls run the stock pipeline; pivoting = :bk selects the Bunch–Kaufman pivoting of the stock L D Lᵀ on every backend (Symmetric indefinite systems), and fused = false forces the stock pipeline. Both pipelines run the static L D Lᵀ with matching (the "matching_alg" of the analysis: "algo5" and "algo6" scale the matrix and give the 2×2 pivot pairs of the ordering) and pivot replacement (pivot_epsilon, pivot_sign: a pivot below pivot_epsilon becomes ±pivot_epsilon, with the requested sign). On the fused pipeline a NaN or infinite pivot always fails the factorization (info ≠ 0), where the stock kernel replaces a NaN pivot and accepts an infinite one. Matching and the pivot pairs are decided from the values of the analysis: when the solver is built before the matrix holds its first values, as an interior-point method builds it, pass deferred_analysis = true. ir_n_steps = k follows every solve! with k steps of plain iterative refinement on either pipeline, the residual taken with the caller's matrix values (unscaled, without the replaced pivots), which recovers the accuracy that pivot replacement costs; the steps run without a host synchronization.

sgn = vcat(ones(Int8, n_primal), fill(Int8(-1), n_dual))   # the sign a replaced pivot takes
S = KKTSolver(Kd, 'L'; factorization = :ldl, matching = "algo5", pivot_epsilon = 1e-10,
              pivot_sign = sgn, ir_n_steps = 2,             # two refinement steps per solve
              deferred_analysis = true)                     # options checked, nothing analyzed yet
Kd.nzVal .= first_values
factorize!(S)                                  # analysis (matching, pairs) and first factorization
Kd.nzVal .= new_values
factorize!(S)                                  # refactorization

The performance settings are an SDSTuning passed as tuning (Fused pipeline).

inertia and pipeline are public but not exported (MadNLP has an inertia of its own): write SparseDirectSolver.inertia or using SparseDirectSolver: inertia.

SparseDirectSolver.KKTSolver — Type
KKTSolver(A, view::Char = 'L'; kwargs...)
KKTSolver(rowptr, colval, nzval, view::Char = 'L'; kwargs...)

Solver for a real symmetric matrix that is factorized many times with the same sparsity pattern and solved for one right-hand side at a time, as the KKT systems of an interior-point method are. A is anything DirectSolver accepts (a CSR, a SparseMatrixCSC, a vendor CSR or CSC matrix such as CuSparseMatrixCSC, ROCSparseMatrixCSC or oneSparseMatrixCSC) and view the stored triangle ('L', 'U' or 'F'); a CSC matrix, SparseMatrixCSC included, is read without a copy, as the CSR of its transpose (a SparseMatrixCSC with spare capacity in rowval/nzval raises InvalidValueError, see csr_of_transpose). The backend is the one of the arrays.

The constructor chooses the pipeline, sets the parameters, runs the analysis and factorizes the values A holds; with deferred_analysis = true it stops after the parameters, and the first factorize! runs the analysis on the values present then. The values alias the caller's nzval: write the new values into it and call factorize!, then solve! and inertia.

K = KKTSolver(A, 'L'; factorization = :ldl)   # analysis and first factorization
nonzeros(A) .= new_values                      # same pattern
factorize!(K)
getparam(K, "info") == 0 || error("factorization failed")
npos, nzero, nneg = SparseDirectSolver.inertia(K)
solve!(K, xb)                                  # xb holds b, then x

Keywords:

  • factorization::Symbol = :cholesky: :cholesky factors A = L Lᵀ (structure "SPD"), :ldl factors A = L D Lᵀ (structure "S") and gives the inertia from D.
  • pivoting::Symbol = :none: :ldl only. :none factors in the pivot order of the analysis (static pivoting: a pivot below pivot_epsilon is replaced); :bk uses Bunch–Kaufman pivoting inside the fronts, which the stock pipeline provides.
  • matching::AbstractString = "none": :ldl only, the "matching_alg" of the analysis ("none" or "default" off, "algo1"–"algo6"; "algo5" and "algo6" also scale the matrix and choose the 2×2 pivot pairs of the ordering from the cycles of the matching, with either pivoting; with "algo1"–"algo4" and pivoting = :none the stock analysis looks for no pairs, as without matching), on both pipelines. It reads the values present at the analysis, so pass deferred_analysis = true when the matrix holds no real values yet.
  • fused::Union{Nothing, Bool} = nothing: true runs the fused pipeline (a FusedSolver), false the stock one (a DirectSolver). nothing selects the fused pipeline where it runs: pivoting = :none on a backend with forward progress between workgroups (CUDA, AMDGPU, oneAPI). The CPU backend, Metal and pivoting = :bk get the stock pipeline.
  • pivot_threshold = nothing: :bk only, the "pivot_threshold" of the acceptance test.
  • pivot_epsilon = nothing: the "pivot_epsilon" below which a pivot is replaced (nothing: 1e-13 for Float64, 1e-5 for Float32, as in cuDSS), on both pipelines. With :ldl and pivoting = :none a pivot d with |d| < pivot_epsilon becomes ±pivot_epsilon; 0 replaces none, and a zero pivot then fails the factorization. On the fused pipeline a NaN or infinite pivot always fails the factorization (info ≠ 0).
  • pivot_sign = nothing: :ldl only, a vector of n entries in (-1, 0, 1), the "pivot_sign" a replaced pivot takes (an interior-point method passes 1 for the primal and -1 for the dual rows; 0 or nothing: the sign of the pivot), on both pipelines.
  • ir_n_steps::Union{Nothing, Integer} = nothing: the "ir_n_steps" of the solve, on both pipelines: every solve! is followed by that many steps of plain iterative refinement against the matrix values the caller passed (not the scaled values of a matching, not the replaced pivots), which recovers accuracy lost to static pivot replacement. nothing or 0: no refinement. The steps run without a host synchronization; getparam(K, "ir_n_steps") reports the steps of the last solve.
  • reordering::String = "algo4": the "reordering_alg" of the analysis ("algo4" nested dissection, "algo3" AMD, "default" the choice of the cost model).
  • user_perm = nothing: a 1-based permutation of length n that replaces the reordering.
  • tuning::Union{Nothing, SDSTuning} = nothing: the performance settings, see SDSTuning. The fused pipeline takes default_tuning of the backend for nothing; the stock pipeline keeps the defaults of DirectSolver for nothing.
  • deferred_analysis::Bool = false: true leaves the analysis and the first factorization to the first factorize!, for a matrix whose values are not assembled yet when the solver is built (the matching and the 2×2 pivot pairs are decided from the values of the analysis). Until then getparam(K, "info") is 0, and solve! and inertia raise FactorizationError.

nothing keeps the default of the parameter. The inner solver is the field inner; getparam is forwarded to it, and pipeline tells which one runs.

Raises InvalidValueError for an unknown factorization, pivoting or matching, for :bk without :ldl or with fused = true, for matching without :ldl, for pivot_threshold without :bk and for pivot_sign without :ldl; NotSupportedError for a complex element type, a uniform batch and fused = true on a backend that cannot run the fused pipeline. Every option is checked by the constructor, also with deferred_analysis = true.

source
SparseDirectSolver.factorize! — Method
factorize!(K::KKTSolver) -> K

Factorize the values the matrix of K holds now, with the analysis of the constructor. With deferred_analysis = true the first call runs the analysis (ordering, matching, 2×2 pivot pairs) on these values and then the first factorization; the later calls refactorize with that analysis. The call returns when the factorization is complete. A failed factorization does not throw: getparam(K, "info") is then the column of the first failed pivot, and 0 otherwise.

source
SparseDirectSolver.solve! — Method
solve!(K::KKTSolver, xb::AbstractVector) -> xb
solve!(K::KKTSolver, x::AbstractVector, b::AbstractVector) -> x

Solve A x = b with the last factorization; with one vector, xb holds b on entry and x on return. The vectors are on the backend of K. The call returns when the solution is complete. Before the first factorization (deferred_analysis = true) it raises FactorizationError.

source
SparseDirectSolver.inertia — Function
inertia(numeric) -> (npos, nneg)

Inertia of D of an LDLᵀ/LDLᴴ factorization (PLAN §1.4 inertia). With perturbed pivots this is the inertia of A + E; read it together with npivots.

source
inertia(K::KKTSolver) -> (npos, nzero, nneg)

The numbers of positive, zero and negative eigenvalues of the matrix K factorized last, which sum to n: for factorization = :ldl from the signs of D (2×2 blocks included), for :cholesky the number of positive pivots and 0 negative ones. It is the inertia of the factorized matrix, so a replaced pivot counts with the sign it was given. Valid when getparam(K, "info") is 0; before the first factorization (deferred_analysis = true) it raises FactorizationError.

source