diff --git a/Project.toml b/Project.toml index f38ccf43..f50bb653 100644 --- a/Project.toml +++ b/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" diff --git a/ext/CliqueTreesExt.jl b/ext/CliqueTreesExt.jl new file mode 100644 index 00000000..9da117be --- /dev/null +++ b/ext/CliqueTreesExt.jl @@ -0,0 +1,5 @@ +module CliqueTreesExt + +include("./directldl_cliquetrees.jl") + +end diff --git a/ext/directldl_cliquetrees.jl b/ext/directldl_cliquetrees.jl new file mode 100644 index 00000000..530c1575 --- /dev/null +++ b/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