Lhongpei · Lhongpei · Aug 12, 2026 · Aug 12, 2026 · Aug 12, 2026
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -3,6 +3,11 @@
 # -----------------------------------------------------------------------------
 cmake_minimum_required(VERSION 3.20)
 
+set(_PDHCG_CUDA_ARCHITECTURES_EXPLICIT FALSE)
+if(DEFINED CMAKE_CUDA_ARCHITECTURES OR DEFINED ENV{CUDAARCHS})
+    set(_PDHCG_CUDA_ARCHITECTURES_EXPLICIT TRUE)
+endif()
+
 project(pdhcg LANGUAGES C CXX CUDA)
 
 set(PDHCG_VERSION_MAJOR 0)
@@ -26,9 +31,17 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
     set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
 endif()
 
-if(NOT DEFINED CMAKE_CUDA_ARCHITECTURES)
-  set(CMAKE_CUDA_ARCHITECTURES 60 70 75 80 86 89 90)
+if(NOT _PDHCG_CUDA_ARCHITECTURES_EXPLICIT)
+    if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 13.0)
+        set(_PDHCG_DEFAULT_CUDA_ARCHITECTURES all)
+    else()
+        set(_PDHCG_DEFAULT_CUDA_ARCHITECTURES 60 70 75 80 86 89 90)
+    endif()
+    set(CMAKE_CUDA_ARCHITECTURES "${_PDHCG_DEFAULT_CUDA_ARCHITECTURES}"
+        CACHE STRING "CUDA architectures to build" FORCE)
+    unset(_PDHCG_DEFAULT_CUDA_ARCHITECTURES)
 endif()
+unset(_PDHCG_CUDA_ARCHITECTURES_EXPLICIT)
 
 # Global Compile flags (corresponding to CFLAGS/NVCCFLAGS)
 add_compile_options(-fPIC -O3 -Wall -Wextra -g)
@@ -162,6 +175,7 @@ endif()
 set(CORE_LINK_LIBS PUBLIC
   CUDA::cudart
   CUDA::cublas
+  CUDA::cusolver
   CUDA::cusparse
   ZLIB::ZLIB
 )
@@ -323,7 +337,7 @@ if(PDHCG_BUILD_TESTS)
 
         # Input-driven inspection/export tools are built but require explicit
         # command-line data, so they are not zero-argument CTest cases.
-        if(TEST_NAME STREQUAL "test_distributed_conic" AND PDHCG_COMPILE_DISTRIBUTED)
+        if(TEST_NAME MATCHES "^test_distributed_(conic|psd)$" AND PDHCG_COMPILE_DISTRIBUTED)
             add_test(NAME ${TEST_NAME}
                      COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 ${MPIEXEC_PREFLAGS}
                              $<TARGET_FILE:${TEST_NAME}> ${MPIEXEC_POSTFLAGS})

diff --git a/README.md b/README.md
@@ -27,7 +27,7 @@ PDHCG solves convex quadratic conic programs in the following form, with a spars
 - $Q$ is the sparse symmetric quadratic component (optional).
 - $R \in \mathbb{R}^{k\times n}$ is a tall low-rank factor (optional, $k$ = rank).
 - $D \in \mathbb{R}^{k\times k}$ is an optional middle matrix that scales / weights / signs the low-rank term. When omitted it defaults to the identity, recovering the standard $Q + R^\top R$ formulation. $D$ may be **diagonal, sparse, dense, or indefinite** — the backend auto-detects the cheapest runtime representation.
-- Standard SOC, Rotated SOC, Exponential, and Power cones are supported both on variable blocks and through native affine constraints $Fx + g \in \mathcal{K}_a$.
+- Standard SOC, Rotated SOC, Exponential, Power, and positive-semidefinite (PSD) cones are supported both on variable blocks and through native affine constraints $Fx + g \in \mathcal{K}_a$.
 
 
 ## Installation (C++ Executable)
@@ -99,7 +99,7 @@ Solver Parameters:
 | --pock_chambolle_alpha | double | Value for Pock-Chambolle step size parameter $\alpha$. | 1.0 |
 | --no_pock_chambolle | flag | Disable Pock-Chambolle rescaling (enabled by default). | false |
 | --no_bound_obj_rescaling | flag | Disable bound objective rescaling (enabled by default). | false |
-| --no_cone_preserving_scaling | flag | Keep coordinate-wise scaling within cone blocks. | false |
+| --no_cone_preserving_scaling | flag | Disable cone-preserving scaling. | false |
 | --sv_max_iter | int | Max iterations for singular value estimation (Power Method). | 5000 |
 | --sv_tol | double | Tolerance for singular value estimation. | 1e-4 |
 | --eval_freq | int | Frequency of termination criteria evaluation (in iterations). | 200 |
@@ -132,7 +132,9 @@ The block scale is
 | Pock-Chambolle | `d_rms` | `sqrt(d_max * d_rms)` |
 
 PDHCG applies the rule to both variable and affine cone blocks. Setting
-`--no_cone_preserving_scaling` bypasses block aggregation.
+`--no_cone_preserving_scaling` bypasses block aggregation. PSD blocks
+remain cone-compatible through a structured element-wise scale
+$s_{ij}=d_i d_j$.
 
 **Distributed Options** (only available when built with `-DPDHCG_COMPILE_DISTRIBUTED=ON`):
 | Option | Type | Description | Default |
@@ -238,9 +240,8 @@ problem = cp.Problem(cp.Minimize(x), [x >= 1])
 problem.solve(solver="PDHCG", eps=1e-6)
 ```
 
-The backend supports quadratic objectives and CVXPY Zero, NonNeg, SOC,
-ExpCone, and PowCone3D constraints. PSD and mixed-integer models are not
-supported.
+The backend supports quadratic objectives and CVXPY Zero, NonNeg, SOC, PSD,
+ExpCone, and PowCone3D constraints.
 
 ## Citation
 If you use this software or method in your research, please cite our paper:

diff --git a/distributed/distributed_solver.cu b/distributed/distributed_solver.cu
@@ -19,7 +19,7 @@ limitations under the License.
 #include "distributed_utils.h"
 #include "internal_types.h"
 #include "pdhcg.h"
-#include "pdhcg_kernels.cuh"
+#include "pdhcg_kernels.h"
 #include "pdhg_core_op.h"
 #include "permute.h"
 #include "preconditioner.h"

diff --git a/docs/C_API.md b/docs/C_API.md
@@ -46,7 +46,7 @@ The objective minimized is `0.5 * x^T (Q + R^T D R) x + c^T x + c0`. `Q`, `R`, a
 - `var_lb`: Variable lower bounds. If `NULL`, defaults to all `-INFINITY`.
 - `var_ub`: Variable upper bounds. If `NULL`, defaults to all `+INFINITY`.
 - `objective_constant`: Scalar constant term added to the objective value. If `NULL`, defaults to `0.0`.
-- `num_var_cones`, `var_cones`: Optional conic variable blocks. Supported types are SOC, rotated SOC, exponential, and power cones. Pass `0` and `NULL` when no variable cones are present. See [Types](c/types.md) for slot layouts and `set_cone_fixed` in [Functions](c/functions.md) for pinning individual slots.
+- `num_var_cones`, `var_cones`: Optional conic variable blocks. Supported types are SOC, rotated SOC, exponential, power, and positive-semidefinite cones. Pass `0` and `NULL` when no variable cones are present. See [Types](c/types.md) for slot layouts and `set_cone_fixed` in [Functions](c/functions.md) for pinning individual non-PSD slots.
 - `affine_cone_matrix_desc`: Matrix descriptor for `F` in the native constraint `F x + affine_cone_offset in K`. Pass `NULL` when no affine cone rows are present.
 - `affine_cone_offset`: Offset vector with one entry per row of `F`. Pass `NULL` for zero offsets.
 - `num_affine_cones`, `affine_cones`: Cone blocks covering every row of `F`. Each `start_idx` is relative to `F`; blocks must be disjoint. Affine cone descriptors must set `is_fixed = NULL`.

diff --git a/docs/c/functions.md b/docs/c/functions.md
@@ -92,7 +92,7 @@ int set_cone_fixed(
 
 Pins one slot of cone `cone_idx` to `value`. Allocates the `is_fixed` flag array on first use and also writes `primal_start[start_idx + slot] = value` so the projection sees the constant. During preprocessing, that slot is also converted to equal lower and upper bounds. Typical use: fix the `y` slot of an exponential cone (e.g. Fisher-market entropy term with `y = 1`).
 
-Variable SOC, rotated-SOC, exponential, and power cones support every fixed-slot pattern whose intersection with the cone is nonempty. The solver validates the section before preprocessing and rejects empty or non-finite sections. Projection and stationarity residuals use the same weighted fixed-section operator, including diagonal quadratic objectives and large SOC/rotated-SOC blocks.
+Variable SOC, rotated-SOC, exponential, and power cones support every fixed-slot pattern whose intersection with the cone is nonempty. The solver validates the section before preprocessing and rejects empty or non-finite sections. Projection and stationarity residuals use the same weighted fixed-section operator, including diagonal quadratic objectives and large SOC/rotated-SOC blocks. `set_cone_fixed` rejects PSD blocks.
 
 **Parameters:**
 

diff --git a/docs/c/types.md b/docs/c/types.md
@@ -103,7 +103,9 @@ typedef enum {
   CONE_ROTATED_SOC = 0,
   CONE_STANDARD_SOC = 1,
   CONE_EXPONENTIAL = 2,
-  CONE_POWER = 3
+  CONE_POWER = 3,
+  CONE_PSD = 4,
+  NUM_CONE_TYPES = 5
 } cone_type_t;
 ```
 
@@ -113,6 +115,12 @@ typedef enum {
 | `CONE_ROTATED_SOC` | `\|\|v\|\|^2 <= 2 s t`, `s, t >= 0` | `v` (`v_dim`), `s`, `t` |
 | `CONE_EXPONENTIAL` | `y * exp(x / y) <= z`, `y > 0` | `x`, `y`, `z` (`v_dim` must be 1) |
 | `CONE_POWER` | `x^alpha * y^(1-alpha) >= \|z\|`, `x,y >= 0` | `x`, `y`, `z` (`v_dim` must be 1) |
+| `CONE_PSD` | Symmetric matrix `X` is positive semidefinite | `svec(X)` (`v_dim` is the matrix order; length `v_dim * (v_dim + 1) / 2`) |
+
+PSD `svec` stores the lower triangle in column-major order. Diagonal entries
+are unchanged and off-diagonal entries are multiplied by `sqrt(2)`.
+Distributed partitioning keeps each PSD block on one GPU and never splits its
+`svec` coordinates across devices.
 
 ## Cone Spec
 
@@ -128,12 +136,12 @@ typedef struct {
 
 Input descriptor for a single cone block. In `var_cones`, `start_idx` indexes
 the variable vector; in `affine_cones`, it indexes rows of the separately
-supplied affine matrix `F`. The
-slot count is `v_dim + 2` for SOC/RSOC and `3` for
-exponential/power cones. Power cones require `power_alpha` in `(0,1)`.
-Variable cones may provide an `is_fixed` array of `slot_count` bytes. Every
-mathematically nonempty fixed-slot pattern is supported for all four cone
-types. Affine cones must set `is_fixed` to NULL.
+supplied affine matrix `F`. The slot count is `v_dim + 2` for SOC/RSOC, `3`
+for exponential/power cones, and `v_dim * (v_dim + 1) / 2` for PSD cones.
+Power cones require `power_alpha` in `(0,1)`. Variable non-PSD cones may
+provide an `is_fixed` array of `slot_count` bytes. Every mathematically
+nonempty fixed-slot pattern is supported for those cone types. PSD blocks do
+not support fixed slots, and affine cones must set `is_fixed` to NULL.
 
 ## Cone Blocks
 

diff --git a/docs/index.md b/docs/index.md
@@ -26,13 +26,13 @@ Where:
 - $c$ is the linear objective vector
 - $\ell_c, u_c$ are constraint bounds
 - $\ell_v, u_v$ are variable bounds
-- $\mathcal{K}_a$ and $\mathcal{K}_v$ are products of Standard SOC, Rotated SOC, Exponential, or Power cones
+- $\mathcal{K}_a$ and $\mathcal{K}_v$ are products of Standard SOC, Rotated SOC, Exponential, Power, or positive-semidefinite cones
 
 ## Key Features
 
 - **GPU Acceleration**: Fully leverages NVIDIA CUDA for extreme-scale QP problems
 - **Flexible Problem Structure**: Supports sparse, low-rank, and middle-weighted low-rank ($R^\top D R$) quadratic terms — alone or combined
-- **Conic constraints**: fully GPU-accelerated SOC, Rotated SOC, Exponential, and Power cone projection on variable blocks or native affine maps $Fx + g$
+- **Conic constraints**: fully GPU-accelerated SOC, Rotated SOC, Exponential, Power, and positive-semidefinite cone projection on variable blocks or native affine maps $Fx + g$
 - **High Performance**: Competitive with commercial solvers on large-scale problems
 - **SpMVOp Auto-Detection**: Automatically uses cuSPARSE SpMVOp on CUDA 13+ while falling back to standard SpMV on CUDA 12.x
 - **Multi-GPU Distributed Solving**: Supports parallel solving across multiple GPUs via MPI and NCCL (optional, enabled at compile time)

diff --git a/docs/python/model.md b/docs/python/model.md
@@ -36,9 +36,9 @@ require one Python dict per cone.
 
 | Field | Type | Notes |
 |---|---|---|
-| `types` | scalar or `int32[K]` | `ConeType.SOC`, `RSOC`, `EXP`, or `POWER`. |
+| `types` | scalar or `int32[K]` | `ConeType.SOC`, `RSOC`, `EXP`, `POWER`, or `PSD`. |
 | `starts` | `int32[K]` | First variable index or affine row of each block. |
-| `v_dims` | scalar or `int32[K]` | Length of `v`; defaults to 1. |
+| `v_dims` | scalar or `int32[K]` | Length of `v`, or matrix order for PSD; defaults to 1. |
 | `power_alphas` | scalar or `float64[K]` | Required in `(0, 1)` for power cones. |
 | `fixed_mask` | optional `uint8[N]` | Ambient-coordinate mask for fixed variable slots. Values come from the primal warm start. |
 
@@ -48,6 +48,14 @@ Slot layout per cone:
 - `rsoc`: `v[0..v_dim-1], s, t` with `||v||^2 <= 2 s t`, `s, t >= 0`.
 - `exp`: `x, y, z` with `y * exp(x / y) <= z`, `y > 0`.
 - `power`: `x, y, z` with `x^alpha * y^(1-alpha) >= |z|`, `x, y >= 0`.
+- `psd`: `svec(X)` for an order-`v_dim` symmetric matrix `X >= 0`. `svec`
+  stores the lower triangle in column-major order, leaves diagonal entries
+  unchanged, and multiplies off-diagonal entries by `sqrt(2)`.
+
+PSD blocks do not support entries in `fixed_mask`; express fixed matrix entries
+as ordinary linear equalities instead.
+In distributed solves, each PSD block remains on one GPU; permutation and
+partitioning never split its `svec` coordinates across devices.
 
 ```python
 import numpy as np

diff --git a/docs/python/quickstart.md b/docs/python/quickstart.md
@@ -100,8 +100,8 @@ print(problem.status, value, x.value)
 ```
 
 The backend preserves CVXPY's primal and dual conventions. It supports
-quadratic objectives and Zero, NonNeg, SOC, ExpCone, and PowCone3D
-constraints. PSD and mixed-integer models are not supported.
+quadratic objectives and Zero, NonNeg, SOC, PSD, ExpCone, and PowCone3D
+constraints. Mixed-integer models are not supported.
 
 ## Model Creation
 

diff --git a/include/pdhcg_types.h b/include/pdhcg_types.h
@@ -66,14 +66,15 @@ extern "C"
         CONE_STANDARD_SOC = 1,
         CONE_EXPONENTIAL = 2,
         CONE_POWER = 3, /* 3-dim: x^alpha * y^(1-alpha) >= |z|, x,y >= 0 */
-        NUM_CONE_TYPES = 4
+        CONE_PSD = 4,   /* svec(X), X symmetric positive semidefinite */
+        NUM_CONE_TYPES = 5
     } cone_type_t;
 
     typedef struct
     {
         int num_cones;
         int *start_idx;      /* [num_cones] */
-        int *v_dim;          /* [num_cones] */
+        int *v_dim;          /* [num_cones]; PSD stores the matrix order */
         cone_type_t *type;   /* [num_cones] */
         double *power_alpha; /* [num_cones]; alpha in (0,1) for CONE_POWER, else unused */
         int fixed_mask_size; /* number of entries in is_fixed; zero when no mask is stored */
@@ -83,10 +84,10 @@ extern "C"
     typedef struct
     {
         cone_type_t type;
-        int start_idx; /* variable index, or row of F for affine cones */
-        int v_dim;
+        int start_idx;        /* variable index, or row of F for affine cones */
+        int v_dim;            /* vector dimension for SOC/RSOC; matrix order for PSD */
         double power_alpha;   /* required for CONE_POWER (in (0,1)); ignored otherwise */
-        const char *is_fixed; /* variable cones only; must be NULL for affine cones */
+        const char *is_fixed; /* variable non-PSD cones only; must be NULL for affine cones */
     } cone_spec_t;
 
     typedef struct

diff --git a/internal/cone_kernel_ops.h b/internal/cone_kernel_ops.h
@@ -0,0 +1,108 @@
+/*
+Copyright 2026 Hongpei Li
+
+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.
+*/
+
+#pragma once
+
+#include "internal_types.h"
+
+typedef void (*cone_proj_launcher_t)(double *primal,
+                                     const double *variable_rescaling,
+                                     double *warm_start,
+                                     const int *start_idx,
+                                     const int *v_dim,
+                                     const double *power_alpha,
+                                     const char *is_fixed,
+                                     int count);
+
+typedef void (*cone_dual_res_launcher_t)(double *dual_residual,
+                                         double *complementarity_residual,
+                                         const double *objective_vector,
+                                         const double *dual_product,
+                                         const double *variable_rescaling,
+                                         const double *primal_solution,
+                                         double *warm_start,
+                                         const int *start_idx,
+                                         const int *v_dim,
+                                         const double *power_alpha,
+                                         const char *is_fixed,
+                                         int count);
+
+typedef void (*cone_proj_diag_q_launcher_t)(double *pdhg_primal,
+                                            double *reflected_primal,
+                                            const double *current_primal,
+                                            const double *variable_rescaling,
+                                            const double *q_diag,
+                                            double tau,
+                                            double *warm_start,
+                                            const int *start_idx,
+                                            const int *v_dim,
+                                            const double *power_alpha,
+                                            const char *is_fixed,
+                                            int count);
+
+typedef struct
+{
+    cone_proj_launcher_t project[NUM_PROJ_METHODS];
+    cone_proj_diag_q_launcher_t project_diag_q[NUM_PROJ_METHODS];
+    cone_dual_res_launcher_t dual_residual[NUM_PROJ_METHODS];
+} cone_kernel_ops_t;
+
+extern const cone_kernel_ops_t pdhcg_rsoc_cone_kernel_ops;
+extern const cone_kernel_ops_t pdhcg_soc_cone_kernel_ops;
+extern const cone_kernel_ops_t pdhcg_exp_cone_kernel_ops;
+extern const cone_kernel_ops_t pdhcg_power_cone_kernel_ops;
+
+void launch_block_projected_mapping_only_dual(double *dual_residual,
+                                              double *complementarity_residual,
+                                              const double *objective_vector,
+                                              const double *dual_product,
+                                              const double *variable_rescaling,
+                                              const double *primal_solution,
+                                              double *warm_start,
+                                              const int *start_idx,
+                                              const int *v_dim,
+                                              const double *power_alpha,
+                                              const char *is_fixed,
+                                              int count);
+
+void launch_grid_projected_mapping_only_dual(double *dual_residual,
+                                             double *complementarity_residual,
+                                             const double *objective_vector,
+                                             const double *dual_product,
+                                             const double *variable_rescaling,
+                                             const double *primal_solution,
+                                             double *warm_start,
+                                             const int *start_idx,
+                                             const int *v_dim,
+                                             const double *power_alpha,
+                                             const char *is_fixed,
+                                             int count);
+
+void launch_cone_reflection(cone_proj_method_t method,
+                            double *reflected_primal,
+                            const double *pdhg_primal,
+                            const double *current_primal,
+                            const int *start_idx,
+                            const int *v_dim,
+                            int count);
+
+void launch_cone_dual_slack(cone_proj_method_t method,
+                            double *dual_slack,
+                            const double *objective_vector,
+                            const double *dual_product,
+                            const int *start_idx,
+                            const int *v_dim,
+                            int count);