Add geometric-mean scaling

README.md

@@ -108,7 +108,8 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the
 | `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` |
 | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` |
 | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` |
-| `--l_inf_ruiz_iter` | `int` | Iterations for L-inf Ruiz rescaling| `10` |
+| `--geo_mean_iter` | `int` | Iterations of geometric-mean scaling | `12` |
+| `--l_inf_ruiz_iter` | `int` | Iterations of L-inf Ruiz rescaling| `10` |
 | `--no_pock_chambolle` | `flag` | Disable Pock-Chambolle rescaling | `enabled` |
 | `--pock_chambolle_alpha` | `float` | Value for Pock-Chambolle alpha | `1.0` |
 | `--no_bound_obj_rescaling` | `flag` | Disable bound objective rescaling | `enabled` |

include/cupdlpx_types.h

@@ -109,6 +109,7 @@ extern "C"
         bool presolve;
         double matrix_zero_tol;
         double infinite_bound;
+        int geometric_mean_iterations;
     } pdhg_parameters_t;
 
     typedef struct

internal/cuda_to_hip.h

@@ -70,6 +70,11 @@ limitations under the License.
 // Device synchronization
 #define cudaDeviceSynchronize hipDeviceSynchronize
 
+// Warp-level primitives. AMD wavefronts are 64 lanes, so CUDA's 32-bit full
+// mask would be rejected; the width-limited non-sync shuffle has the same
+// sub-group semantics on HIP.
+#define __shfl_down_sync(mask, var, delta, width) __shfl_down(var, delta, width)
+
 // ----------------------------------------------------------------------------
 // CUDA Graph API -> HIP Graph API
 // ----------------------------------------------------------------------------

python/README.md

@@ -163,7 +163,8 @@ Below is a list of commonly used parameters, their internal keys, and descriptio
 | `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. |
 | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. |
 | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. |
-| `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations for L∞ Ruiz scaling. Improves numerical conditioning. |
+| `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. |
+| `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations of L∞ Ruiz scaling. Improves numerical conditioning. |
 | `UsePCAlpha` | `has_pock_chambolle_alpha` | bool | `True` | Whether to use the Pock–Chambolle α step size adjustment. |
 | `PCAlpha` | `pock_chambolle_alpha` | float | `1.0` | Value of the Pock–Chambolle α parameter. |
 | `BoundObjRescaling` | `bound_objective_rescaling` | bool | `True` | Whether to rescale the objective vector during preprocessing. |

python/cupdlpx/PDLP.py

@@ -42,6 +42,7 @@
     "OptimalityTol": "eps_optimal_relative",
     "FeasibilityTol": "eps_feasible_relative",
     # scaling / step size
+    "GeoMeanIters": "geometric_mean_iterations",
     "RuizIters": "l_inf_ruiz_iterations",
     "UsePCAlpha": "has_pock_chambolle_alpha",
     "PCAlpha": "pock_chambolle_alpha",

python/cupdlpx/model.py

@@ -43,6 +43,7 @@
     {
         "termination_evaluation_frequency",
         "iteration_limit",
+        "geometric_mean_iterations",
         "l_inf_ruiz_iterations",
         "sv_max_iter",
     }

python_bindings/_core_bindings.cpp

@@ -279,6 +279,7 @@ static py::dict get_default_params_py()
     d["iteration_limit"] = p.termination_criteria.iteration_limit;
 
     // rescaling
+    d["geometric_mean_iterations"] = p.geometric_mean_iterations;
     d["l_inf_ruiz_iterations"] = p.l_inf_ruiz_iterations;
     d["has_pock_chambolle_alpha"] = p.has_pock_chambolle_alpha;
     d["pock_chambolle_alpha"] = p.pock_chambolle_alpha;
@@ -389,6 +390,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
     geti("iteration_limit", p->termination_criteria.iteration_limit);
 
     // rescaling
+    geti("geometric_mean_iterations", p->geometric_mean_iterations);
     geti("l_inf_ruiz_iterations", p->l_inf_ruiz_iterations);
     getb("has_pock_chambolle_alpha", p->has_pock_chambolle_alpha);
     getf("pock_chambolle_alpha", p->pock_chambolle_alpha);
@@ -423,6 +425,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
         throw std::invalid_argument("termination_evaluation_frequency must be positive.");
     if (p->termination_criteria.iteration_limit < 0)
         throw std::invalid_argument("iteration_limit must be nonnegative.");
+    if (p->geometric_mean_iterations < 0)
+        throw std::invalid_argument("geometric_mean_iterations must be nonnegative.");
     if (p->l_inf_ruiz_iterations < 0)
         throw std::invalid_argument("l_inf_ruiz_iterations must be nonnegative.");
     if (p->sv_max_iter <= 0)

src/cli.c

@@ -176,9 +176,12 @@ void print_usage(const char *prog_name)
     fprintf(stderr,
             "      --eps_feas <tolerance>          "
             "Relative feasibility tolerance (default: 1e-4).\n");
+    fprintf(stderr,
+            "      --geo_mean_iter <int>           "
+            "Iterations of geometric-mean scaling (default: 12).\n");
     fprintf(stderr,
             "      --l_inf_ruiz_iter <int>         "
-            "Iterations for L-inf Ruiz rescaling (default: 10).\n");
+            "Iterations of L-inf Ruiz rescaling (default: 10).\n");
     fprintf(stderr,
             "      --no_pock_chambolle             "
             "Disable Pock-Chambolle rescaling (default: enabled).\n");
@@ -231,6 +234,7 @@ int main(int argc, char *argv[])
                                            {"eps_feas", required_argument, 0, 1004},
                                            {"eps_feas_polish", required_argument, 0, 1006},
                                            {"feasibility_polishing", no_argument, 0, 'f'},
+                                           {"geo_mean_iter", required_argument, 0, 1018},
                                            {"l_inf_ruiz_iter", required_argument, 0, 1007},
                                            {"pock_chambolle_alpha", required_argument, 0, 1008},
                                            {"no_pock_chambolle", no_argument, 0, 1009},
@@ -276,6 +280,9 @@ int main(int argc, char *argv[])
             case 'f': // --feasibility_polishing
                 params.feasibility_polishing = true;
                 break;
+            case 1018: // --geo_mean_iter
+                params.geometric_mean_iterations = atoi(optarg);
+                break;
             case 1007: // --l_inf_ruiz_iter
                 params.l_inf_ruiz_iterations = atoi(optarg);
                 break;

src/preconditioner.cu

@@ -23,6 +23,10 @@ limitations under the License.
 
 #define SCALING_EPSILON 1e-12
 
+#define GEOMETRIC_MEAN_MULTIPLIER_MAX 1e20
+
+#define MAX_LOG2_VECTOR_WIDTH 5
+
 __global__ void scale_variables_kernel(double *__restrict__ objective_vector,
                                        double *__restrict__ variable_lower_bound,
                                        double *__restrict__ variable_upper_bound,
@@ -75,6 +79,18 @@ __global__ void scale_objective_kernel(double *__restrict__ objective_vector,
                                        double constraint_scale,
                                        double objective_scale);
 __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables);
+__global__ void max_csr_row_length_kernel(const int *__restrict__ row_ptr, int num_rows, int *__restrict__ result);
+__global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr,
+                                                const int *__restrict__ col_ind,
+                                                const double *__restrict__ matrix_vals,
+                                                int num_rows,
+                                                int log2_width,
+                                                const double *__restrict__ other_multiplier,
+                                                double *__restrict__ multiplier);
+__global__ void geometric_mean_finalize_kernel(const double *__restrict__ multiplier,
+                                               double *__restrict__ scaling_factors,
+                                               double *__restrict__ cumulative_rescaling,
+                                               int n);
 static void scale_problem(pdhg_solver_state_t *state,
                           double *constraint_rescaling,
                           double *variable_rescaling,
@@ -87,6 +103,13 @@ static void ruiz_rescaling(pdhg_solver_state_t *state,
                            double *variable_rescaling,
                            double *inverse_constraint_rescaling,
                            double *inverse_variable_rescaling);
+static void geometric_mean_rescaling(pdhg_solver_state_t *state,
+                                     int num_iterations,
+                                     rescale_info_t *rescale_info,
+                                     double *constraint_rescaling,
+                                     double *variable_rescaling,
+                                     double *inverse_constraint_rescaling,
+                                     double *inverse_variable_rescaling);
 static void pock_chambolle_rescaling(pdhg_solver_state_t *state,
                                      const double alpha,
                                      rescale_info_t *rescale_info,
@@ -158,6 +181,88 @@ static void ruiz_rescaling(pdhg_solver_state_t *state,
     }
 }
 
+static int log2_vector_width(long long num_nonzeros, int num_rows, int longest_row)
+{
+    if (longest_row >= 512)
+        return MAX_LOG2_VECTOR_WIDTH;
+    if (num_rows <= 0)
+        return 0;
+    const long long mean_nnz = (num_nonzeros + num_rows - 1) / num_rows;
+    int log2_width = 0;
+    while (log2_width < MAX_LOG2_VECTOR_WIDTH && (1LL << log2_width) < mean_nnz)
+        ++log2_width;
+    return log2_width;
+}
+
+static void longest_csr_rows(const pdhg_solver_state_t *state, int *longest)
+{
+    int *device_longest = NULL;
+    CUDA_CHECK(cudaMalloc(&device_longest, 2 * sizeof(int)));
+    CUDA_CHECK(cudaMemset(device_longest, 0, 2 * sizeof(int)));
+
+    max_csr_row_length_kernel<<<state->num_blocks_dual, THREADS_PER_BLOCK>>>(
+        state->constraint_matrix->row_ptr, state->num_constraints, device_longest);
+    max_csr_row_length_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
+        state->constraint_matrix_t->row_ptr, state->num_variables, device_longest + 1);
+
+    CUDA_CHECK(cudaMemcpy(longest, device_longest, 2 * sizeof(int), cudaMemcpyDeviceToHost));
+    CUDA_CHECK(cudaFree(device_longest));
+}
+
+// Geometric-mean scaling: each row/column is divided by sqrt(min * max) of the
+// absolute values of its nonzeros, alternating between rows and columns.
+// Reference: J.A. Tomlin, "On Scaling Linear Programming Problems,"
+// Mathematical Programming Study 4 (1975) 146-166.
+static void geometric_mean_rescaling(pdhg_solver_state_t *state,
+                                     int num_iterations,
+                                     rescale_info_t *rescale_info,
+                                     double *constraint_rescaling,
+                                     double *variable_rescaling,
+                                     double *inverse_constraint_rescaling,
+                                     double *inverse_variable_rescaling)
+{
+    const int num_constraints = state->num_constraints;
+    const int num_variables = state->num_variables;
+
+    fill_ones_kernel<<<state->num_blocks_dual, THREADS_PER_BLOCK>>>(inverse_constraint_rescaling, num_constraints);
+    fill_ones_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(inverse_variable_rescaling, num_variables);
+
+    int longest[2] = {0, 0};
+    longest_csr_rows(state, longest);
+    const long long num_nonzeros = state->constraint_matrix->num_nonzeros;
+    const int row_log2 = log2_vector_width(num_nonzeros, num_constraints, longest[0]);
+    const int col_log2 = log2_vector_width(num_nonzeros, num_variables, longest[1]);
+    const int row_blocks =
+        (int)((((long long)num_constraints << row_log2) + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK);
+    const int col_blocks = (int)((((long long)num_variables << col_log2) + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK);
+
+    for (int iter = 0; iter < num_iterations; ++iter)
+    {
+        geometric_mean_iteration_kernel<<<row_blocks, THREADS_PER_BLOCK>>>(state->constraint_matrix->row_ptr,
+                                                                           state->constraint_matrix->col_ind,
+                                                                           state->constraint_matrix->val,
+                                                                           num_constraints,
+                                                                           row_log2,
+                                                                           inverse_variable_rescaling,
+                                                                           inverse_constraint_rescaling);
+        geometric_mean_iteration_kernel<<<col_blocks, THREADS_PER_BLOCK>>>(state->constraint_matrix_t->row_ptr,
+                                                                           state->constraint_matrix_t->col_ind,
+                                                                           state->constraint_matrix_t->val,
+                                                                           num_variables,
+                                                                           col_log2,
+                                                                           inverse_constraint_rescaling,
+                                                                           inverse_variable_rescaling);
+    }
+
+    geometric_mean_finalize_kernel<<<state->num_blocks_dual, THREADS_PER_BLOCK>>>(
+        inverse_constraint_rescaling, constraint_rescaling, rescale_info->con_rescale, num_constraints);
+    geometric_mean_finalize_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
+        inverse_variable_rescaling, variable_rescaling, rescale_info->var_rescale, num_variables);
+
+    scale_problem(
+        state, constraint_rescaling, variable_rescaling, inverse_constraint_rescaling, inverse_variable_rescaling);
+}
+
 static void pock_chambolle_rescaling(pdhg_solver_state_t *state,
                                      const double alpha,
                                      rescale_info_t *rescale_info,
@@ -260,6 +365,20 @@ rescale_info_t *rescale_problem(const pdhg_parameters_t *params, pdhg_solver_sta
     CUDA_CHECK(cudaMalloc(&inverse_constraint_rescaling, num_constraints * sizeof(double)));
     CUDA_CHECK(cudaMalloc(&inverse_variable_rescaling, num_variables * sizeof(double)));
 
+    if (params->geometric_mean_iterations > 0)
+    {
+        if (params->verbose)
+        {
+            printf("  Geometric-mean scaling (%d iterations)\n", params->geometric_mean_iterations);
+        }
+        geometric_mean_rescaling(state,
+                                 params->geometric_mean_iterations,
+                                 rescale_info,
+                                 constraint_rescaling,
+                                 variable_rescaling,
+                                 inverse_constraint_rescaling,
+                                 inverse_variable_rescaling);
+    }
     if (params->l_inf_ruiz_iterations > 0)
     {
         if (params->verbose)
@@ -496,3 +615,63 @@ __global__ void fill_ones_kernel(double *__restrict__ x, int num_variables)
     if (i < num_variables)
         x[i] = 1.0;
 }
+
+__global__ void max_csr_row_length_kernel(const int *__restrict__ row_ptr, int num_rows, int *__restrict__ result)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < num_rows)
+        atomicMax(result, row_ptr[i + 1] - row_ptr[i]);
+}
+
+__global__ void geometric_mean_iteration_kernel(const int *__restrict__ row_ptr,
+                                                const int *__restrict__ col_ind,
+                                                const double *__restrict__ matrix_vals,
+                                                int num_rows,
+                                                int log2_width,
+                                                const double *__restrict__ other_multiplier,
+                                                double *__restrict__ multiplier)
+{
+    const int width = 1 << log2_width;
+    const int groups_per_block = blockDim.x >> log2_width;
+    const int row = blockIdx.x * groups_per_block + (threadIdx.x >> log2_width);
+    const int lane = threadIdx.x & (width - 1);
+
+    double lo = INFINITY, hi = 0.0;
+    if (row < num_rows)
+    {
+        const int start = row_ptr[row], end = row_ptr[row + 1];
+        for (int k = start + lane; k < end; k += width)
+        {
+            const double scaled = fabs(matrix_vals[k]) * other_multiplier[col_ind[k]];
+            if (!(scaled > 0.0 && isfinite(scaled)))
+                continue;
+            lo = fmin(lo, scaled);
+            hi = fmax(hi, scaled);
+        }
+    }
+
+    for (int offset = width >> 1; offset > 0; offset >>= 1)
+    {
+        lo = fmin(lo, __shfl_down_sync(0xffffffffu, lo, offset, width));
+        hi = fmax(hi, __shfl_down_sync(0xffffffffu, hi, offset, width));
+    }
+
+    if (lane != 0 || row >= num_rows || hi <= 0.0)
+        return;
+
+    const double updated = 1.0 / (sqrt(lo) * sqrt(hi));
+    multiplier[row] = fmin(fmax(updated, 1.0 / GEOMETRIC_MEAN_MULTIPLIER_MAX), GEOMETRIC_MEAN_MULTIPLIER_MAX);
+}
+
+__global__ void geometric_mean_finalize_kernel(const double *__restrict__ multiplier,
+                                               double *__restrict__ scaling_factors,
+                                               double *__restrict__ cumulative_rescaling,
+                                               int n)
+{
+    for (int t = blockIdx.x * blockDim.x + threadIdx.x; t < n; t += blockDim.x * gridDim.x)
+    {
+        const double divisor = 1.0 / multiplier[t];
+        cumulative_rescaling[t] *= divisor;
+        scaling_factors[t] = divisor;
+    }
+}

src/utils.cu

@@ -301,6 +301,7 @@ bool should_do_adaptive_restart(pdhg_solver_state_t *solver_state,
 
 void set_default_parameters(pdhg_parameters_t *params)
 {
+    params->geometric_mean_iterations = 12;
     params->l_inf_ruiz_iterations = 10;
     params->has_pock_chambolle_alpha = true;
     params->pock_chambolle_alpha = 1.0;
@@ -610,6 +611,7 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro
         printf("  optimality_norm    : %s\n", params->optimality_norm == NORM_TYPE_L_INF ? "L_inf" : "L2");
     }
 
+    PRINT_DIFF_INT("geo_mean_iter", params->geometric_mean_iterations, default_params.geometric_mean_iterations);
     PRINT_DIFF_INT("l_inf_ruiz_iter", params->l_inf_ruiz_iterations, default_params.l_inf_ruiz_iterations);
     PRINT_DIFF_DBL("pock_chambolle_alpha", params->pock_chambolle_alpha, default_params.pock_chambolle_alpha);
     PRINT_DIFF_BOOL(

test/test_api_surface.py

@@ -12,8 +12,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import os
-
 import numpy as np
 import scipy.sparse as sp
 import pytest
@@ -412,18 +410,3 @@ def test_optimize_rejects_bad_sense(base_lp_data):
     model._model_sense = 999  # bypass the property to hit optimize()'s guard
     with pytest.raises(ValueError):
         model.optimize()
-
-
-def test_read_mps_roundtrip():
-    """read() loads an MPS file into a usable Model."""
-    path = os.path.join(os.path.dirname(__file__), "cplex2.mps")
-    model = read(path)
-    assert model.num_vars > 0
-    assert model.num_constrs > 0
-    assert model.A.shape == (model.num_constrs, model.num_vars)
-    assert model.ModelSense in (PDLP.MINIMIZE, PDLP.MAXIMIZE)
-
-
-def test_read_mps_missing_file():
-    with pytest.raises(FileNotFoundError):
-        read("this_file_does_not_exist.mps")