From af4f3401a7854d0cddb1ba93b530e294b24ba506 Mon Sep 17 00:00:00 2001 From: fpacaud Date: Wed, 10 Apr 2024 14:21:47 +0200 Subject: [PATCH 1/2] [LinearSolvers] Add support for LDL factorization in CHOLMOD --- src/LinearSolvers/cholmod.jl | 60 +++++++++++++++++++++++++++++++++--- test/madnlp_test.jl | 58 +++++++++++++++++++--------------- 2 files changed, 90 insertions(+), 28 deletions(-) diff --git a/src/LinearSolvers/cholmod.jl b/src/LinearSolvers/cholmod.jl index 85570fa5..c93b3b41 100644 --- a/src/LinearSolvers/cholmod.jl +++ b/src/LinearSolvers/cholmod.jl @@ -1,5 +1,5 @@ -# TODO: check options in CHOLMOD @kwdef mutable struct CHOLMODOptions <: AbstractOptions + cholmod_algorithm::LinearFactorization = CHOLESKY end mutable struct CHOLMODSolver{T} <: AbstractLinearSolver{T} @@ -32,7 +32,6 @@ function CHOLMODSolver( ) full.nzval .= 1.0 - # TODO: use AMD permutation here A = CHOLMOD.Sparse(full) inner = CHOLMOD.symbolic(A) @@ -42,7 +41,11 @@ end function factorize!(M::CHOLMODSolver) M.full.nzval .= M.tril_to_full_view # We check the factorization succeeded later in the backsolve - CHOLMOD.cholesky!(M.inner, M.full; check=false) + if M.opt.cholmod_algorithm == LDL + CHOLMOD.ldlt!(M.inner, M.full; check=false) + elseif M.opt.cholmod_algorithm == CHOLESKY + CHOLMOD.cholesky!(M.inner, M.full; check=false) + end return M end @@ -57,8 +60,28 @@ function solve!(M::CHOLMODSolver{T}, rhs::Vector{T}) where T return rhs end +# Utils function to copy the diagonal entries directly from CHOLMOD's factor. +function _get_diagonal!(F::CHOLMOD.Factor{T}, d::Vector{T}) where T + s = unsafe_load(CHOLMOD.typedpointer(F)) + # Wrap in memory the factor LD stored in CHOLMOD. + colptr = unsafe_wrap(Array, s.p, s.n+1, own=false) + nz = unsafe_wrap(Array, s.nz, s.n, own=false) + rowval = unsafe_wrap(Array, s.i, s.nzmax, own=false) + nzvals = unsafe_wrap(Array, Ptr{T}(s.x), s.nzmax, own=false) + # Read LD factor and load diagonal entries + for j in 1:s.n + for c in colptr[j]:colptr[j]+nz[j]-1 + i = rowval[c+1] + 1 # convert to 1-based indexing + if i == j + d[i] = nzvals[c+1] + end + end + end + return d +end + is_inertia(::CHOLMODSolver) = true -function inertia(M::CHOLMODSolver) +function _inertia_cholesky(M::CHOLMODSolver) n = size(M.full, 1) if issuccess(M.inner) return (n, 0, 0) @@ -66,6 +89,35 @@ function inertia(M::CHOLMODSolver) return (0, n, 0) end end +function _inertia_ldlt(M::CHOLMODSolver{T}) where T + n = size(M.full, 1) + if !issuccess(M.inner) + return (0, n, 0) + end + D = M.d + # Extract diagonal elements + _get_diagonal!(M.inner, D) + (pos, zero, neg) = (0, 0, 0) + @inbounds for i in 1:n + d = D[i] + if d > 0 + pos += 1 + elseif d == 0 + zero += 1 + else + neg += 1 + end + end + @assert pos + zero + neg == n + return pos, zero, neg +end +function inertia(M::CHOLMODSolver) + if M.opt.cholmod_algorithm == CHOLESKY + return _inertia_cholesky(M) + elseif M.opt.cholmod_algorithm == LDL + return _inertia_ldlt(M) + end +end input_type(::Type{CHOLMODSolver}) = :csc default_options(::Type{CHOLMODSolver}) = CHOLMODOptions() diff --git a/test/madnlp_test.jl b/test/madnlp_test.jl index 2ed6e760..de23eaa6 100644 --- a/test/madnlp_test.jl +++ b/test/madnlp_test.jl @@ -1,13 +1,27 @@ testset = [ [ - "Umfpack", + "SparseKKTSystem + Umfpack", ()->MadNLP.Optimizer( linear_solver=MadNLP.UmfpackSolver, print_level=MadNLP.ERROR), [] ], [ - "LapackCPU-BUNCHKAUFMAN", + "SparseKKTSystem + InertiaFree", + ()->MadNLP.Optimizer( + inertia_correction_method=MadNLP.InertiaFree, + print_level=MadNLP.ERROR), + [] + ], + [ + "SparseKKTSystem + RelaxBound", + ()->MadNLP.Optimizer( + fixed_variable_treatment=MadNLP.RelaxBound, + print_level=MadNLP.ERROR), + [] + ], + [ + "DenseKKTSystem + LapackCPU-BUNCHKAUFMAN", ()->MadNLP.Optimizer( kkt_system=MadNLP.DenseKKTSystem, linear_solver=MadNLP.LapackCPUSolver, @@ -16,7 +30,7 @@ testset = [ [] ], [ - "LapackCPU-LU", + "DenseKKTSystem + LapackCPU-LU", ()->MadNLP.Optimizer( kkt_system=MadNLP.DenseKKTSystem, linear_solver=MadNLP.LapackCPUSolver, @@ -25,7 +39,7 @@ testset = [ [] ], [ - "LapackCPU-QR", + "DenseKKTSystem + LapackCPU-QR", ()->MadNLP.Optimizer( kkt_system=MadNLP.DenseKKTSystem, linear_solver=MadNLP.LapackCPUSolver, @@ -34,7 +48,7 @@ testset = [ [] ], [ - "LapackCPU-CHOLESKY", + "DenseKKTSystem + LapackCPU-CHOLESKY", ()->MadNLP.Optimizer( kkt_system=MadNLP.DenseKKTSystem, linear_solver=MadNLP.LapackCPUSolver, @@ -43,52 +57,48 @@ testset = [ ["infeasible", "lootsma", "eigmina"] ], [ - "Option: RELAX_BOUND", + "SparseUnreducedKKTSystem", ()->MadNLP.Optimizer( - fixed_variable_treatment=MadNLP.RelaxBound, + kkt_system=MadNLP.SparseUnreducedKKTSystem, print_level=MadNLP.ERROR), [] ], [ - "Option: AUGMENTED KKT SYSTEM", + "SparseUnreducedKKTSystem + InertiaFree", ()->MadNLP.Optimizer( + inertia_correction_method=MadNLP.InertiaFree, kkt_system=MadNLP.SparseUnreducedKKTSystem, print_level=MadNLP.ERROR), [] ], [ - "Option: SPARSE CONDENSED KKT SYSTEM", + "SparseCondensedKKTSystem + CHOLMOD", ()->MadNLP.Optimizer( kkt_system=MadNLP.SparseCondensedKKTSystem, equality_treatment = MadNLP.RelaxEquality, fixed_variable_treatment = MadNLP.RelaxBound, - tol = 1e-4, - print_level=MadNLP.ERROR), - [] - ], - [ - "Option: InertiaFree", - ()->MadNLP.Optimizer( - inertia_correction_method=MadNLP.InertiaFree, - print_level=MadNLP.ERROR), + linear_solver=MadNLP.CHOLMODSolver, + print_level=MadNLP.INFO), [] ], [ - "Option: InertiaFree & AUGMENTED KKT SYSTEM", + "SparseCondensedKKTSystem + CHOLMOD-LDL", ()->MadNLP.Optimizer( - inertia_correction_method=MadNLP.InertiaFree, - kkt_system=MadNLP.SparseUnreducedKKTSystem, - print_level=MadNLP.ERROR), + kkt_system=MadNLP.SparseCondensedKKTSystem, + equality_treatment = MadNLP.RelaxEquality, + fixed_variable_treatment = MadNLP.RelaxBound, + linear_solver=MadNLP.CHOLMODSolver, + cholmod_algorithm=MadNLP.LDL, + print_level=MadNLP.INFO), [] ], [ - "Option: InertiaFree & SPARSE CONDENSED KKT SYSTEM", + "SparseCondensedKKTSystem + InertiaFree", ()->MadNLP.Optimizer( inertia_correction_method=MadNLP.InertiaFree, kkt_system=MadNLP.SparseCondensedKKTSystem, equality_treatment = MadNLP.RelaxEquality, fixed_variable_treatment = MadNLP.RelaxBound, - tol = 1e-4, print_level=MadNLP.ERROR), [] ], From d48e9931049fb4a089ceaefab3e3f9419997d9fd Mon Sep 17 00:00:00 2001 From: fpacaud Date: Fri, 26 Apr 2024 16:26:39 +0200 Subject: [PATCH 2/2] fix test for julia <= v1.9 --- src/LinearSolvers/cholmod.jl | 6 +++++- test/madnlp_test.jl | 33 ++++++++++++++++++++------------- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/src/LinearSolvers/cholmod.jl b/src/LinearSolvers/cholmod.jl index c93b3b41..89c53081 100644 --- a/src/LinearSolvers/cholmod.jl +++ b/src/LinearSolvers/cholmod.jl @@ -23,6 +23,10 @@ function CHOLMODSolver( d = Vector{Float64}(undef,csc.n) full, tril_to_full_view = get_tril_to_full(Float64, csc) + if opt.cholmod_algorithm == LDL && VERSION <= v"1.9" + error("[CHOLMOD] Option `cholmod_algorithm=LDL` is not supported for Julia version <= 1.9") + end + full = SparseMatrixCSC{Float64,Int}( full.m, full.n, @@ -30,7 +34,7 @@ function CHOLMODSolver( Vector{Int64}(full.rowval), full.nzval ) - full.nzval .= 1.0 + fill!(full.nzval, one(T)) A = CHOLMOD.Sparse(full) inner = CHOLMOD.symbolic(A) diff --git a/test/madnlp_test.jl b/test/madnlp_test.jl index de23eaa6..195df4d6 100644 --- a/test/madnlp_test.jl +++ b/test/madnlp_test.jl @@ -72,24 +72,13 @@ testset = [ [] ], [ - "SparseCondensedKKTSystem + CHOLMOD", + "SparseCondensedKKTSystem + CHOLMOD-CHOLESKY", ()->MadNLP.Optimizer( kkt_system=MadNLP.SparseCondensedKKTSystem, equality_treatment = MadNLP.RelaxEquality, fixed_variable_treatment = MadNLP.RelaxBound, linear_solver=MadNLP.CHOLMODSolver, - print_level=MadNLP.INFO), - [] - ], - [ - "SparseCondensedKKTSystem + CHOLMOD-LDL", - ()->MadNLP.Optimizer( - kkt_system=MadNLP.SparseCondensedKKTSystem, - equality_treatment = MadNLP.RelaxEquality, - fixed_variable_treatment = MadNLP.RelaxBound, - linear_solver=MadNLP.CHOLMODSolver, - cholmod_algorithm=MadNLP.LDL, - print_level=MadNLP.INFO), + print_level=MadNLP.ERROR), [] ], [ @@ -104,6 +93,24 @@ testset = [ ], ] +# N.B. Current CHOLMOD interface is supported only starting from Julia v1.10. +if VERSION >= v"1.10" + push!( + testset, + [ + "SparseCondensedKKTSystem + CHOLMOD-LDL", + ()->MadNLP.Optimizer( + kkt_system=MadNLP.SparseCondensedKKTSystem, + equality_treatment = MadNLP.RelaxEquality, + fixed_variable_treatment = MadNLP.RelaxBound, + linear_solver=MadNLP.CHOLMODSolver, + cholmod_algorithm=MadNLP.LDL, + print_level=MadNLP.ERROR), + [] + ] + ) +end + for (name,optimizer_constructor,exclude) in testset test_madnlp(name,optimizer_constructor,exclude)