CuClarabel.jl interface (cvxpy#2812)

PTNobel · web-flow · commit a6bdefb9fb70 · 2025-06-26T13:29:22.000-07:00
* Work on interface

* Finishes interface

* Deletes dead code

* Fixes CI

* Work on CuClarabel

* Starts docs

* Adds docs deletes dead code
diff --git a/cvxpy/__init__.py b/cvxpy/__init__.py
@@ -63,6 +63,7 @@
 from cvxpy.settings import (
     CBC as CBC,
     CLARABEL as CLARABEL,
+    CUCLARABEL as CUCLARABEL,
     COPT as COPT,
     CPLEX as CPLEX,
     CPP_CANON_BACKEND as CPP_CANON_BACKEND,
diff --git a/cvxpy/reductions/solvers/conic_solvers/cuclarabel_conif.py b/cvxpy/reductions/solvers/conic_solvers/cuclarabel_conif.py
@@ -0,0 +1,215 @@
+"""
+Copyright 2022, the CVXPY Authors
+
+Licensed under the Apache License, Version 2.0 (the "License");
+
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+"""
+import numpy as np
+import scipy.sparse as sp
+
+import cvxpy.settings as s
+from cvxpy.constraints import SOC, ExpCone, PowCone3D
+from cvxpy.reductions.solution import Solution, failure_solution
+from cvxpy.reductions.solvers import utilities
+from cvxpy.reductions.solvers.conic_solvers.conic_solver import ConicSolver
+from cvxpy.utilities.citations import CITATION_DICT
+
+
+def dims_to_solver_cones(jl, cone_dims):
+
+    jl.seval("""cones = Clarabel.SupportedCone[]""")
+
+    # assume that constraints are presented
+    # in the preferred ordering of SCS.
+
+    if cone_dims.zero > 0:
+        jl.push_b(jl.cones, jl.Clarabel.ZeroConeT(cone_dims.zero))
+
+    if cone_dims.nonneg > 0:
+        jl.push_b(jl.cones, jl.Clarabel.NonnegativeConeT(cone_dims.nonneg))
+
+    for dim in cone_dims.soc:
+        jl.push_b(jl.cones, jl.Clarabel.SecondOrderConeT(dim))
+
+    for dim in cone_dims.psd:
+        jl.push_b(jl.cones, jl.Clarabel.PSDTriangleConeT(dim))
+
+    for _ in range(cone_dims.exp):
+        jl.push_b(jl.cones, jl.Clarabel.ExponentialConeT())
+
+    for pow in cone_dims.p3d:
+        jl.push_b(jl.cones, jl.Clarabel.PowerConeT(pow))
+
+
+class CUCLARABEL(ConicSolver):
+    """An interface for the Clarabel solver.
+    """
+
+    # Solver capabilities.
+    MIP_CAPABLE = False
+    SUPPORTED_CONSTRAINTS = ConicSolver.SUPPORTED_CONSTRAINTS \
+        + [SOC, ExpCone, PowCone3D]
+
+    STATUS_MAP = {
+                    "SOLVED": s.OPTIMAL,
+                    "PRIMAL_INFEASIBLE": s.INFEASIBLE,
+                    "DUAL_INFEASIBLE": s.UNBOUNDED,
+                    "ALMOST_SOLVED": s.OPTIMAL_INACCURATE,
+                    "ALMOST_PRIMAL_INFEASIBLE": s.INFEASIBLE_INACCURATE,
+                    "Almost_DUAL_INFEASIBLE": s.UNBOUNDED_INACCURATE,
+                    "MAX_ITERATIONS": s.USER_LIMIT,
+                    "MAX_TIME": s.USER_LIMIT,
+                    "NUMERICAL_ERROR": s.SOLVER_ERROR,
+                    "INSUFFICIENT_PROGRESS": s.SOLVER_ERROR
+                }
+
+    # Order of exponential cone arguments for solver.
+    EXP_CONE_ORDER = [0, 1, 2]
+
+    def name(self):
+        """The name of the solver.
+        """
+        return 'CUCLARABEL'
+
+    def import_solver(self) -> None:
+        """Imports the solver.
+        """
+        # Does NOT import juliacall because that starts a Julia interpreter
+        # and this method is called on CVXPY importing and that's too heavy.
+        import cupy  # noqa F401
+
+    def supports_quad_obj(self) -> bool:
+        """Clarabel supports quadratic objective with any combination
+        of conic constraints.
+        """
+        return True
+
+    @staticmethod
+    def extract_dual_value(result_vec, offset, constraint):
+        """Extracts the dual value for constraint starting at offset.
+        """
+        return utilities.extract_dual_value(result_vec, offset, constraint)
+
+    def invert(self, solution, inverse_data):
+        """Returns the solution to the original problem given the inverse_data.
+        """
+
+        attr = {}
+        status = self.STATUS_MAP[str(solution.status)]
+        attr[s.SOLVE_TIME] = solution.solve_time
+        attr[s.NUM_ITERS] = solution.iterations
+        # more detailed statistics here when available
+        # attr[s.EXTRA_STATS] = solution.extra.FOO
+
+        if status in s.SOLUTION_PRESENT:
+            primal_val = solution.obj_val
+            opt_val = primal_val + inverse_data[s.OFFSET]
+            primal_vars = {
+                inverse_data[CUCLARABEL.VAR_ID]: np.array(solution.x)
+            }
+            eq_dual_vars = utilities.get_dual_values(
+                np.array(solution.z[:inverse_data[ConicSolver.DIMS].zero]),
+                self.extract_dual_value,
+                inverse_data[CUCLARABEL.EQ_CONSTR]
+            )
+            ineq_dual_vars = utilities.get_dual_values(
+                np.array(solution.z[inverse_data[ConicSolver.DIMS].zero:]),
+                self.extract_dual_value,
+                inverse_data[CUCLARABEL.NEQ_CONSTR]
+            )
+            dual_vars = {}
+            dual_vars.update(eq_dual_vars)
+            dual_vars.update(ineq_dual_vars)
+            return Solution(status, opt_val, primal_vars, dual_vars, attr)
+        else:
+            return failure_solution(status, attr)
+
+    def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):
+        """Returns the result of the call to the solver.
+
+        Parameters
+        ----------
+        data : dict
+            Data generated via an apply call.
+        warm_start : Bool
+            Whether to warm_start Clarabel.
+        verbose : Bool
+            Control the verbosity.
+        solver_opts : dict
+            Clarabel-specific solver options.
+
+        Returns
+        -------
+        The result returned by a call to clarabel.solve().
+        """
+        import cupy
+        from cupyx.scipy.sparse import csr_matrix as cucsr_matrix
+        from juliacall import Main as jl
+        jl.seval('using Clarabel, LinearAlgebra, SparseArrays')
+        jl.seval('using CUDA, CUDA.CUSPARSE')
+
+        A = data[s.A]
+        b = data[s.B]
+        q = data[s.C]
+
+        if s.P in data:
+            P = data[s.P]
+        else:
+            nvars = q.size
+            P = sp.csr_array((nvars, nvars))
+
+        P = sp.triu(P).tocsr()
+
+        cones = data[ConicSolver.DIMS]
+
+        Pgpu = cucsr_matrix(P)
+        qgpu = cupy.array(q)
+
+        Agpu = cucsr_matrix(A)
+        bgpu = cupy.array(b)
+
+        if Pgpu.nnz != 0:
+            jl.P = jl.Clarabel.cupy_to_cucsrmat(
+                jl.Float64, int(Pgpu.data.data.ptr), int(Pgpu.indices.data.ptr),
+                int(Pgpu.indptr.data.ptr), *Pgpu.shape, Pgpu.nnz)
+        else:
+            jl.seval(f"""
+            P = CuSparseMatrixCSR(sparse(Float64[], Float64[], Float64[], {nvars}, {nvars}))
+            """)
+        jl.q = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(qgpu.data.ptr), qgpu.size)
+
+        jl.A = jl.Clarabel.cupy_to_cucsrmat(
+                jl.Float64, int(Agpu.data.data.ptr), int(Agpu.indices.data.ptr),
+                int(Agpu.indptr.data.ptr), *Agpu.shape, Agpu.nnz)
+        jl.b = jl.Clarabel.cupy_to_cuvector(jl.Float64, int(bgpu.data.ptr), bgpu.size)
+
+
+        dims_to_solver_cones(jl, cones)
+
+        results = jl.seval("""
+        settings = Clarabel.Settings(direct_solve_method = :cudss)
+        solver   = Clarabel.Solver(P,q,A,b,cones, settings)
+        Clarabel.solve!(solver)
+        """)
+        return results
+    
+    def cite(self, data):
+        """Returns bibtex citation for the solver.
+
+        Parameters
+        ----------
+        data : dict
+            Data generated via an apply call.
+        """
+        return CITATION_DICT["CUCLARABEL"]
diff --git a/cvxpy/reductions/solvers/defines.py b/cvxpy/reductions/solvers/defines.py
@@ -22,6 +22,7 @@
 from cvxpy.reductions.solvers.conic_solvers.clarabel_conif import CLARABEL as CLARABEL_con
 from cvxpy.reductions.solvers.conic_solvers.copt_conif import COPT as COPT_con
 from cvxpy.reductions.solvers.conic_solvers.cplex_conif import CPLEX as CPLEX_con
+from cvxpy.reductions.solvers.conic_solvers.cuclarabel_conif import CUCLARABEL as CUCLARABEL_con
 from cvxpy.reductions.solvers.conic_solvers.cvxopt_conif import CVXOPT as CVXOPT_con
 
 # Conic interfaces
@@ -61,7 +62,7 @@
                      GLPK_MI_con(), CBC_con(), CLARABEL_con(), SCS_con(), SDPA_con(),
                      GUROBI_con(), MOSEK_con(), CPLEX_con(), NAG_con(), XPRESS_con(),
                      SCIP_con(), SCIPY_con(), HIGHS_con(), GLOP_con(), PDLP_con(),
-                     QOCO_con(), ECOS_BB_con()]
+                     QOCO_con(), CUCLARABEL_con(), ECOS_BB_con()]
 solver_qp_intf = [OSQP_qp(),
                   GUROBI_qp(),
                   CPLEX_qp(),
@@ -84,7 +85,7 @@
                  s.CPLEX, s.GUROBI, s.COPT, s.GLPK, s.NAG,
                  s.GLPK_MI, s.CBC, s.CVXOPT, s.XPRESS, s.DIFFCP,
                  s.SCIP, s.SCIPY, s.HIGHS, s.GLOP, s.PDLP, s.QOCO, 
-                 s.ECOS_BB]
+                 s.CUCLARABEL, s.ECOS_BB]
 QP_SOLVERS = [s.OSQP,
               s.GUROBI,
               s.CPLEX,
diff --git a/cvxpy/settings.py b/cvxpy/settings.py
@@ -92,13 +92,15 @@
 SCIP = "SCIP"
 SCIPY = "SCIPY"
 CLARABEL = "CLARABEL"
+CUCLARABEL = "CUCLARABEL"
 DAQP = "DAQP"
 HIGHS = "HIGHS"
 MPAX = "MPAX"
 SOLVERS = [CLARABEL, ECOS, CVXOPT, GLOP, GLPK, GLPK_MI,
            SCS, SDPA, GUROBI, OSQP, CPLEX,
            MOSEK, CBC, COPT, XPRESS, PIQP, PROXQP, QOCO,
-           NAG, PDLP, SCIP, SCIPY, DAQP, HIGHS, MPAX]
+           NAG, PDLP, SCIP, SCIPY, DAQP, HIGHS, MPAX,
+           CUCLARABEL]
 
 # Xpress-specific items
 XPRESS_IIS = "XPRESS_IIS"
diff --git a/cvxpy/tests/test_conic_solvers.py b/cvxpy/tests/test_conic_solvers.py
@@ -543,6 +543,114 @@ def test_clarabel_sdp_2(self) -> None:
         sth.check_complementarity(places)
         sth.check_dual_domains(places)
 
+@unittest.skipUnless('CUCLARABEL' in INSTALLED_SOLVERS, 'CLARABEL is not installed.')
+class TestCuClarabel(BaseTest):
+
+    """ Unit tests for Clarabel. """
+    def setUp(self) -> None:
+
+        self.x = cp.Variable(2, name='x')
+        self.y = cp.Variable(3, name='y')
+
+        self.A = cp.Variable((2, 2), name='A')
+        self.B = cp.Variable((2, 2), name='B')
+        self.C = cp.Variable((3, 2), name='C')
+
+    def test_clarabel_parameter_update(self) -> None:
+        """Test warm start.
+        """
+        x = cp.Variable(2)
+        P = cp.Parameter(nonneg=True),
+        A = cp.Parameter(4)
+        b = cp.Parameter(2, nonneg=True)
+        q = cp.Parameter(2)
+
+        def update_parameters(P, A, b, q):
+            P[0].value = np.random.rand()
+            A.value = np.random.randn(4)
+            b.value = np.random.rand(2)
+            q.value = np.random.randn(2)
+
+        prob = cp.Problem(
+                cp.Minimize(P[0]*cp.square(x[0]) + cp.quad_form(x, np.ones([2, 2])) + q.T @ x),
+                [A[0] * x[0] + A[1] * x[1] == b[0],
+                 A[2] * x[0] + A[3] * x[1] <= b[1]]
+            )
+
+        update_parameters(P, A, b, q)
+        result1 = prob.solve(solver=cp.CLARABEL, warm_start=False)
+        result2 = prob.solve(solver=cp.CLARABEL, warm_start=True)
+        self.assertAlmostEqual(result1, result2)
+
+        update_parameters(P, A, b, q)
+        result1 = prob.solve(solver=cp.CLARABEL, warm_start=True)
+        result2 = prob.solve(solver=cp.CLARABEL, warm_start=False)
+        self.assertAlmostEqual(result1, result2)
+
+        # consecutive solves, no data update
+        result1 = prob.solve(solver=cp.CLARABEL, warm_start=False)
+        self.assertAlmostEqual(result1, result2)
+
+
+    def test_clarabel_lp_0(self) -> None:
+        StandardTestLPs.test_lp_0(solver=cp.CUCLARABEL)
+
+    def test_clarabel_nonstandard_name(self) -> None:
+        # Test that solver name with non-standard capitalization works.
+        StandardTestLPs.test_lp_0(solver="CuClarabel")
+
+    def test_clarabel_lp_1(self) -> None:
+        StandardTestLPs.test_lp_1(solver='CUCLARABEL')
+
+    def test_clarabel_lp_2(self) -> None:
+        StandardTestLPs.test_lp_2(solver='CUCLARABEL')
+
+    def test_clarabel_lp_3(self) -> None:
+        StandardTestLPs.test_lp_3(solver='CUCLARABEL')
+
+    def test_clarabel_lp_4(self) -> None:
+        StandardTestLPs.test_lp_4(solver='CUCLARABEL')
+
+    def test_clarabel_lp_5(self) -> None:
+        StandardTestLPs.test_lp_5(solver='CUCLARABEL')
+
+    def test_clarabel_qp_0(self) -> None:
+        StandardTestQPs.test_qp_0(solver='CUCLARABEL')
+
+    def test_clarabel_qp_0_linear_obj(self) -> None:
+        StandardTestQPs.test_qp_0(solver='CUCLARABEL', use_quad_obj=False)
+
+    def test_clarabel_socp_0(self) -> None:
+        StandardTestSOCPs.test_socp_0(solver='CUCLARABEL')
+
+    def test_clarabel_socp_1(self) -> None:
+        StandardTestSOCPs.test_socp_1(solver='CUCLARABEL')
+
+    def test_clarabel_socp_2(self) -> None:
+        StandardTestSOCPs.test_socp_2(solver='CUCLARABEL')
+
+    def test_clarabel_socp_3(self) -> None:
+        # axis 0
+        breakpoint()
+        StandardTestSOCPs.test_socp_3ax0(solver='CUCLARABEL')
+        # axis 1
+        StandardTestSOCPs.test_socp_3ax1(solver='CUCLARABEL')
+
+    def test_clarabel_expcone_1(self) -> None:
+        StandardTestECPs.test_expcone_1(solver='CUCLARABEL')
+
+    def test_clarabel_exp_soc_1(self) -> None:
+        StandardTestMixedCPs.test_exp_soc_1(solver='CUCLARABEL')
+
+    def test_clarabel_pcp_0(self) -> None:
+        StandardTestSOCPs.test_socp_0(solver='CUCLARABEL')
+
+    def test_clarabel_pcp_1(self) -> None:
+        StandardTestSOCPs.test_socp_1(solver='CUCLARABEL')
+
+    def test_clarabel_pcp_2(self) -> None:
+        StandardTestSOCPs.test_socp_2(solver='CUCLARABEL')
+
 
 @unittest.skipUnless('MOSEK' in INSTALLED_SOLVERS, 'MOSEK is not installed.')
 class TestMosek(unittest.TestCase):
diff --git a/doc/source/tutorial/solvers/index.rst b/doc/source/tutorial/solvers/index.rst
@@ -97,6 +97,8 @@ The table below shows the types of problems the supported solvers can handle.
 +----------------+----+----+------+-----+-----+-----+-----+
 | `MPAX`_        | X  | X  |      |     |     |     |     |
 +----------------+----+----+------+-----+-----+-----+-----+
+| `CUCLARABEL`_  | X  | X  | X    |     | X   | X   |     |
++----------------+----+----+------+-----+-----+-----+-----+
 | `CVXOPT`_      | X  | X  | X    | X   |     |     |     |
 +----------------+----+----+------+-----+-----+-----+-----+
 | `SDPA`_ \*\*\* | X  | X  | X    | X   |     |     |     |
