Integrate CliqueTrees.jl

.github/workflows/test.yml

@@ -103,6 +103,37 @@ jobs:
   #     - uses: codecov/codecov-action@v5
   #       with:
   #         file: lcov.info
+  test-cliquetrees:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: ['1']
+        julia-arch: [x64]
+        os: [ubuntu-latest, macos-latest]
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/install-juliaup@v2
+        with:
+          channel: ${{ matrix.julia-version }}
+      - uses: julia-actions/cache@v1
+      - name: test MadNLPCliqueTrees
+        shell: julia --project=./lib/MadNLPCliqueTrees --color=yes {0}
+        run: |
+          using Pkg
+          Pkg.Registry.update()
+          Pkg.develop(path=".")
+          Pkg.instantiate()
+          using MadNLP
+          PATH_MADNLP = dirname(pathof(MadNLP))
+          Pkg.develop(path=joinpath(PATH_MADNLP, "..", "lib", "MadNLPTests"))
+          Pkg.test("MadNLPCliqueTrees", coverage=true)
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: lib/MadNLPCliqueTrees/src
+      - uses: codecov/codecov-action@v5
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
   # need to run one on cuda and one on amdgpu
   test-gpu:
     strategy:

lib/MadNLPCliqueTrees/Project.toml

@@ -0,0 +1,21 @@
+name = "MadNLPCliqueTrees"
+uuid = "b01c5a2c-fd8a-4b14-b48d-9b33b3498a9c"
+version = "0.1.0"
+
+[deps]
+CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[compat]
+CliqueTrees = "1.16.2"
+MadNLP = "0.9"
+julia = "1.10"
+
+[extras]
+MadNLPTests = "b52a2a03-04ab-4a5f-9698-6a2deff93217"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test", "MadNLPTests"]

lib/MadNLPCliqueTrees/src/MadNLPCliqueTrees.jl

@@ -0,0 +1,153 @@
+module MadNLPCliqueTrees
+
+using CliqueTrees
+using CliqueTrees: EliminationAlgorithm, DEFAULT_ELIMINATION_ALGORITHM
+using CliqueTrees.Multifrontal
+using CliqueTrees.Multifrontal: FChordalLDLt, FChordalCholesky, AbstractRegularization, NoRegularization
+using LinearAlgebra
+using LinearAlgebra: PivotingStrategy
+using SparseArrays
+
+import MadNLP:
+    MadNLP,
+    @kwdef,
+    MadNLPLogger,
+    AbstractOptions,
+    AbstractLinearSolver,
+    LinearFactorization,
+    LDL,
+    CHOLESKY,
+    set_options!,
+    introduce,
+    factorize!,
+    solve_linear_system!,
+    improve!,
+    is_inertia,
+    inertia,
+    input_type,
+    is_supported,
+    default_options
+
+@kwdef mutable struct CliqueTreesOptions <: AbstractOptions
+    cliquetrees_algorithm::LinearFactorization = LDL
+    cliquetrees_ordering::EliminationAlgorithm = DEFAULT_ELIMINATION_ALGORITHM # AMF
+    cliquetrees_strategy::PivotingStrategy = choose_strategy(cliquetrees_algorithm)
+    cliquetrees_regularization::AbstractRegularization = NoRegularization()
+end
+
+function choose_strategy(alg::LinearFactorization)
+    if alg == CHOLESKY
+        strategy = NoPivot()
+    elseif alg == LDL
+        strategy = RowMaximum()
+    else
+        error()
+    end
+
+    return strategy
+end
+
+mutable struct CliqueTreesSolver{T, F <: Factorization{T}} <: AbstractLinearSolver{T}
+    tril::SparseMatrixCSC{T, Int32}
+    F::F
+    signs::Vector{Int}
+    opt::CliqueTreesOptions
+    logger::MadNLPLogger
+end
+
+function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{LDL}) where T
+    S = Symmetric(tril, :L)
+    F = ChordalLDLt{:L}(S; alg = opt.cliquetrees_ordering)::FChordalLDLt{:L, T, Int32}
+    return F
+end
+
+function _build_factorization(tril::SparseMatrixCSC{T, Int32}, opt::CliqueTreesOptions, ::Val{CHOLESKY}) where T
+    S = Symmetric(tril, :L)
+    F = ChordalCholesky{:L}(S; alg = opt.cliquetrees_ordering)::FChordalCholesky{:L, T, Int32}
+    return F
+end
+
+function CliqueTreesSolver(
+    tril::SparseMatrixCSC{T, Int32};
+    opt = CliqueTreesOptions(),
+    logger = MadNLPLogger(),
+    pos = size(tril, 1),
+) where T
+    signs = fill(-1, size(tril, 1)); signs[1:pos] .= 1
+    F = _build_factorization(tril, opt, Val(opt.cliquetrees_algorithm))
+    return CliqueTreesSolver{T, typeof(F)}(tril, F, signs, opt, logger)
+end
+
+function factorize!(M::CliqueTreesSolver{T, <:ChordalLDLt}) where T
+    ldlt!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; signs=M.signs, reg=M.opt.cliquetrees_regularization, check=false)
+    return M
+end
+
+function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalLDLt}, rhs::Vector{T}) where T
+    if issuccess(M.F)
+        ldiv!(M.F, rhs)
+    end
+
+    return rhs
+end
+
+function inertia(M::CliqueTreesSolver{T, <:ChordalLDLt}) where T
+    d = M.F.d
+    pos = 0; zer = 0; neg = 0
+    @inbounds for di in d
+        if di > 0
+            pos += 1
+        elseif di == 0
+            zer += 1
+        else
+            neg += 1
+        end
+    end
+    return pos, zer, neg
+end
+
+introduce(::CliqueTreesSolver{T, <:ChordalLDLt}) where T =
+    "CliqueTrees/LDLᵀ v$(pkgversion(CliqueTrees))"
+
+function factorize!(M::CliqueTreesSolver{T, <:ChordalCholesky}) where T
+    cholesky!(copy!(M.F, M.tril), M.opt.cliquetrees_strategy; check=false)
+    return M
+end
+
+function solve_linear_system!(M::CliqueTreesSolver{T, <:ChordalCholesky}, rhs::Vector{T}) where T
+    if issuccess(M.F)
+        ldiv!(M.F, rhs)
+    end
+
+    return rhs
+end
+
+function inertia(M::CliqueTreesSolver{T, <:ChordalCholesky}) where T
+    n = size(M.tril, 1)
+
+    if issuccess(M.F)
+        return (n, 0, 0)
+    else
+        return (0, n, 0)
+    end
+end
+
+introduce(::CliqueTreesSolver{T, <:ChordalCholesky}) where T =
+    "CliqueTrees/Cholesky v$(pkgversion(CliqueTrees))"
+
+is_inertia(::CliqueTreesSolver) = true
+improve!(::CliqueTreesSolver) = false
+input_type(::Type{<:CliqueTreesSolver}) = :csc
+default_options(::Type{<:CliqueTreesSolver}) = CliqueTreesOptions()
+is_supported(::Type{<:CliqueTreesSolver}, ::Type{T}) where T <: AbstractFloat = true
+
+export CliqueTreesSolver, CliqueTreesOptions
+
+for name in names(MadNLP, all = true)
+    if Base.isexported(MadNLP, name)
+        @eval using MadNLP: $(name)
+        @eval export $(name)
+    end
+end
+
+end # module MadNLPCliqueTrees

lib/MadNLPCliqueTrees/test/runtests.jl

@@ -0,0 +1,88 @@
+using Test, MadNLP, MadNLPCliqueTrees, MadNLPTests, CliqueTrees
+using CliqueTrees.Multifrontal: DynamicRegularization, GMW81, SE99
+using SparseArrays
+
+@testset "MadNLPCliqueTrees test" begin
+    # Default LDL algorithm
+    MadNLPTests.test_linear_solver(CliqueTreesSolver, Float32)
+    MadNLPTests.test_linear_solver(CliqueTreesSolver, Float64)
+
+    # Cholesky algorithm on positive definite matrix
+    @testset "Cholesky algorithm" begin
+        for T in (Float32, Float64)
+            row = Int32[1,2,2]; col = Int32[1,1,2]; val = T[1., .1, 2.]
+            b = T[1.0, 3.0]
+            sol = T[0.8542713567839195, 1.4572864321608041]
+            csc = sparse(row, col, val, 2, 2)
+            opt = CliqueTreesOptions(cliquetrees_algorithm=MadNLP.CHOLESKY)
+            M = CliqueTreesSolver(csc; opt=opt)
+            MadNLP.factorize!(M)
+            @test MadNLP.is_inertia(M)
+            @test MadNLP.inertia(M) == (2, 0, 0)
+            x = MadNLP.solve_linear_system!(M, copy(b))
+            @test MadNLPTests.solcmp(x, sol)
+        end
+    end
+
+    test_madnlp(
+        "CliqueTrees",
+        () -> MadNLP.Optimizer(linear_solver = CliqueTreesSolver, print_level = MadNLP.ERROR),
+        ["eigmina"],
+    )
+
+    @testset "Elimination orderings" begin
+        for alg in (AMF(), MMD(), MF())
+            test_madnlp(
+                "CliqueTrees/$(nameof(typeof(alg)))",
+                () -> MadNLP.Optimizer(
+                    linear_solver = CliqueTreesSolver,
+                    print_level = MadNLP.ERROR,
+                    cliquetrees_ordering = alg,
+                ),
+                ["eigmina"],
+            )
+        end
+    end
+
+    @testset "Dynamic Regularization" begin
+        reg = DynamicRegularization()
+
+        test_madnlp(
+            "CliqueTrees/$(nameof(typeof(reg)))",
+            () -> MadNLP.Optimizer(
+                linear_solver = CliqueTreesSolver,
+                print_level = MadNLP.ERROR,
+                cliquetrees_regularization = reg,
+            ),
+            [],
+        )
+    end
+
+    @testset "GMW81" begin
+        reg = GMW81()
+
+        test_madnlp(
+            "CliqueTrees/$(nameof(typeof(reg)))",
+            () -> MadNLP.Optimizer(
+                linear_solver = CliqueTreesSolver,
+                print_level = MadNLP.ERROR,
+                cliquetrees_regularization = reg,
+            ),
+            ["unbounded"],
+        )
+    end
+
+    @testset "SE99" begin
+        reg = SE99()
+
+        test_madnlp(
+            "CliqueTrees/$(nameof(typeof(reg)))",
+            () -> MadNLP.Optimizer(
+                linear_solver = CliqueTreesSolver,
+                print_level = MadNLP.ERROR,
+                cliquetrees_regularization = reg,
+            ),
+            ["eigmina"],
+        )
+    end
+end

src/KKT/Dense/augmented.jl

@@ -81,7 +81,7 @@ function create_kkt_system(
     fill!(diag_hess, zero(T))
 
     quasi_newton = create_quasi_newton(hessian_approximation, cb, n; options=qn_options)
-    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver)
+    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver, pos = n + ns)
 
     return DenseKKTSystem(
         hess, jac, quasi_newton,

src/KKT/Dense/condensed.jl

@@ -94,7 +94,7 @@ function create_kkt_system(
     ind_ineq_shifted = cb.ind_ineq .+ n .+ ns
 
     quasi_newton = create_quasi_newton(hessian_approximation, cb, n; options=qn_options)
-    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver)
+    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver, pos = n + m - ns)
 
     return DenseCondensedKKTSystem(
         hess, jac, quasi_newton, jac_ineq,

src/KKT/Sparse/augmented.jl

@@ -126,7 +126,7 @@ function create_kkt_system(
     hess_com, hess_csc_map = coo_to_csc(hess_raw)
 
     _linear_solver = linear_solver(
-        aug_com; opt = opt_linear_solver
+        aug_com; opt = opt_linear_solver, pos = n_tot
     )
 
     return SparseKKTSystem(

src/KKT/Sparse/condensed.jl

@@ -116,7 +116,7 @@ function create_kkt_system(
         jt_csc
     )
 
-    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver)
+    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver, pos = n_tot)
     ext = get_sparse_condensed_ext(VT, hess_com, jptr, jt_csc_map, hess_csc_map)
     return SparseCondensedKKTSystem(
         hess, hess_raw, hess_com, hess_csc_map,

src/KKT/Sparse/scaled_augmented.jl

@@ -160,7 +160,7 @@ function create_kkt_system(
     hess_com, hess_csc_map = coo_to_csc(hess_raw)
 
     _linear_solver = linear_solver(
-        aug_com; opt = opt_linear_solver
+        aug_com; opt = opt_linear_solver, pos = n_tot
     )
 
     return ScaledSparseKKTSystem(

src/KKT/Sparse/unreduced.jl

@@ -145,7 +145,7 @@ function create_kkt_system(
     jac_com, jac_csc_map = coo_to_csc(jac_raw)
     hess_com, hess_csc_map = coo_to_csc(hess_raw)
 
-    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver)
+    _linear_solver = linear_solver(aug_com; opt = opt_linear_solver, pos = n_tot)
     return SparseUnreducedKKTSystem(
         hess, jac_callback, jac, quasi_newton, reg, pr_diag, du_diag,
         l_diag, u_diag, l_lower, u_lower, l_lower_aug, u_lower_aug,

src/LinearSolvers/cholmod.jl

@@ -17,7 +17,7 @@ end
 
 function CHOLMODSolver(
     csc::SparseMatrixCSC{T};
-    opt=CHOLMODOptions(), logger=MadNLPLogger(),
+    opt=CHOLMODOptions(), logger=MadNLPLogger(), pos=nothing,
 ) where T
     p = Vector{Float64}(undef,csc.n)
     d = Vector{Float64}(undef,csc.n)

src/LinearSolvers/lapack.jl

@@ -22,6 +22,7 @@ mutable struct LapackCPUSolver{T, MT} <: AbstractLinearSolver{T}
         A::MT;
         opt=LapackOptions(),
         logger=MadNLPLogger(),
+        pos=nothing,
     ) where {MT <: AbstractMatrix}
         T = eltype(A)
         m,n = size(A)

src/LinearSolvers/ldl.jl

@@ -13,7 +13,7 @@ end
 
 function LDLSolver(
     tril::SparseMatrixCSC{T};
-    opt=LDLFactorizationsOptions(), logger=MadNLPLogger(),
+    opt=LDLFactorizationsOptions(), logger=MadNLPLogger(), pos=nothing,
 ) where T
     # TODO: convert tril to triu, not full
     full, tril_to_full_view = get_tril_to_full(T,tril)

src/LinearSolvers/mumps.jl

@@ -164,7 +164,7 @@ end
 # ---------------------------------------------------------------------------------------
 
 function MumpsSolver(csc::SparseMatrixCSC{T,Int32};
-    opt=MumpsOptions(), logger=MadNLPLogger(),
+    opt=MumpsOptions(), logger=MadNLPLogger(), pos=nothing,
 ) where T
 
     I,J = findIJ(csc)