[LinearSolvers] Add support for LDL factorization in CHOLMOD (#321)

* [LinearSolvers] Add support for LDL factorization in CHOLMOD

* fix test for julia <= v1.9

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}
@@ -23,16 +23,19 @@ 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,
         Vector{Int64}(full.colptr),
         Vector{Int64}(full.rowval),
         full.nzval
     )
-    full.nzval .= 1.0
+    fill!(full.nzval, one(T))
 
-    # TODO: use AMD permutation here
     A = CHOLMOD.Sparse(full)
     inner = CHOLMOD.symbolic(A)
 
@@ -42,7 +45,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,15 +64,64 @@ 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)
     else
         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()
 

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,57 +57,60 @@ testset = [
         ["infeasible", "lootsma", "eigmina", "lp_examodels_issue75"]
     ],
     [
-        "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-CHOLESKY",
         ()->MadNLP.Optimizer(
             kkt_system=MadNLP.SparseCondensedKKTSystem,
             equality_treatment = MadNLP.RelaxEquality,
             fixed_variable_treatment = MadNLP.RelaxBound,
-            tol = 1e-4,
+            linear_solver=MadNLP.CHOLMODSolver,
             print_level=MadNLP.ERROR),
         []
     ],
     [
-        "Option: InertiaFree",
-        ()->MadNLP.Optimizer(
-            inertia_correction_method=MadNLP.InertiaFree,
-            print_level=MadNLP.ERROR),
-        []
-    ],
-    [
-        "Option: InertiaFree & AUGMENTED KKT SYSTEM",
-        ()->MadNLP.Optimizer(
-            inertia_correction_method=MadNLP.InertiaFree,
-            kkt_system=MadNLP.SparseUnreducedKKTSystem,
-            print_level=MadNLP.ERROR),
-        []
-    ],
-    [
-        "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),
         []
     ],
 ]
 
+# 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)