Add KKT solver: cliquetrees.

Project.toml

@@ -20,15 +20,18 @@ TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [weakdeps]
+CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
 Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
 
 [extensions]
+CliqueTreesExt = "CliqueTrees"
 HSLExt = "HSL"
 PardisoExt = "Pardiso"
 
 [compat]
 AMD = "0.4, 0.5"
+CliqueTrees = "1.15"
 DataStructures = "0.18"
 GenericLinearAlgebra = "0.3"
 HSL = "0.4, 0.5"

ext/CliqueTreesExt.jl

@@ -0,0 +1,5 @@
+module CliqueTreesExt
+
+include("./directldl_cliquetrees.jl")
+
+end

ext/directldl_cliquetrees.jl

@@ -0,0 +1,82 @@
+using CliqueTrees.Multifrontal, SparseArrays, LinearAlgebra, Clarabel
+using CliqueTrees.Multifrontal: FVector, flatindices, setflatindex!
+
+import Clarabel: DefaultInt, AbstractDirectLDLSolver, LinearSolverInfo
+import Clarabel: ldlsolver_constructor, ldlsolver_matrix_shape, ldlsolver_is_available
+import Clarabel: linear_solver_info, update_values!, scale_values!, refactor!, solve!
+
+struct CliqueTreesDirectLDLSolver{T} <: AbstractDirectLDLSolver{T}
+    F::ChordalLDLt{:L, T, Int, FVector{T}, FVector{Int}}
+    P::Vector{Int}
+    reg::DynamicRegularization{T}
+    signs::Vector{Int}
+
+    function CliqueTreesDirectLDLSolver{T}(KKT::SparseMatrixCSC{T}, signs, settings) where {T}
+        A = Symmetric(KKT, :L)
+        F = ChordalLDLt{:L}(A)
+        P = flatindices(F, A)
+
+        reg = DynamicRegularization(;
+            delta=convert(T, settings.dynamic_regularization_delta),
+            epsilon=convert(T, settings.dynamic_regularization_eps),
+        )
+
+        return new{T}(F, P, reg, signs)
+    end
+end
+
+ldlsolver_constructor(::Val{:cliquetrees}) = CliqueTreesDirectLDLSolver
+ldlsolver_matrix_shape(::Val{:cliquetrees}) = :tril
+ldlsolver_is_available(::Val{:cliquetrees}) = true
+
+function linear_solver_info(solver::CliqueTreesDirectLDLSolver{T}) where {T}
+    name = :cliquetrees
+    threads = 1
+    direct = true
+    nnzA = length(solver.P)
+    nnzL = nnz(solver.F)
+    LinearSolverInfo(name, threads, direct, nnzA, nnzL)
+end
+
+function update_values!(
+    solver::CliqueTreesDirectLDLSolver{T},
+    index::AbstractVector{DefaultInt},
+    values::Vector{T}
+) where {T}
+    return
+end
+
+function scale_values!(
+    solver::CliqueTreesDirectLDLSolver{T},
+    index::AbstractVector{DefaultInt},
+    scale::T
+) where {T}
+    return
+end
+
+function refactor!(solver::CliqueTreesDirectLDLSolver{T}, K::SparseMatrixCSC{T}) where {T}
+    F = solver.F
+    P = solver.P
+    reg = solver.reg
+    signs = solver.signs
+    nzval = K.nzval
+
+    fill!(F, zero(T))
+
+    @inbounds for i in eachindex(P)
+        setflatindex!(F, nzval[i], P[i])
+    end
+
+    ldlt!(F; check=false, reg, signs)
+    return issuccess(F)
+end
+
+function solve!(
+    solver::CliqueTreesDirectLDLSolver{T},
+    K::SparseMatrixCSC{T},
+    x::Vector{T},
+    b::Vector{T}
+) where {T}
+    ldiv!(x, solver.F, b)
+    return
+end