[WIP] Speedup

.github/workflows/formatting.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        os: [ubuntu-latest]
+        os: [macos-latest]
 
     steps:
       - uses: actions/checkout@v5

include/PSLP/PSLP_API.h

@@ -102,7 +102,7 @@ extern "C"
     } Presolver;
 
     /* The user is responsible for freeing the settings using 'free_settings'. */
-    Settings *default_settings();
+    Settings *default_settings(void);
     void free_settings(Settings *stgs);
     void set_settings_true(Settings *stgs);
     void set_settings_false(Settings *stgs);

include/PSLP/PSLP_stats.h

@@ -47,6 +47,8 @@ extern "C"
         double time_init;
         double time_fast_reductions;
         double time_medium_reductions;
+        double time_ston_cols;
+        double time_dton_rows;
         double time_primal_propagation;
         double time_parallel_rows;
         double time_parallel_cols;

include/core/Activity.h

@@ -132,8 +132,7 @@ double compute_max_act_one_tag(const double *vals, const int *cols, int len,
 void update_activities_bound_change(Activity *activities, const struct Matrix *A,
                                     const struct Bound *bounds, const double *vals,
                                     const int *rows, int len, double old_bound,
-                                    double new_bound, double finite_bound,
-                                    bool lower,
+                                    double new_bound, bool finite_bound, bool lower,
                                     iVec *updated_activities HUGE_BOUND_PARAM);
 
 void update_activities_fixed_col(Activity *activities, const struct Matrix *A,

include/core/Constraints.h

@@ -90,4 +90,8 @@ void delete_inactive_rows(Constraints *constraints);
 */
 void delete_inactive_cols_from_A_and_AT(Constraints *constraints);
 
+#ifdef TESTING
+size_t update_column_map(const int *col_sizes, int *map, size_t n_cols);
+#endif
+
 #endif // CONSTRAINTS_H

include/core/Debugger.h

@@ -28,7 +28,6 @@ struct Lock;
 struct iVec;
 struct PresolveStats;
 
-#include <assert.h>
 #include <stdbool.h>
 #include <stddef.h>
 

include/core/Matrix.h

@@ -108,8 +108,8 @@ empty column of A, will be removed when this function is called on AT.
 To have consistency between A and AT we must update the column indices of A.
 This is done in the end of the function (columns[j] = colsmap[columns[j]]).
 */
-void remove_extra_space(Matrix *A, const int *row_sizes, const int *col_sizes,
-                        bool remove_all, int *col_idxs_map);
+void remove_extra_space(Matrix *A, const int *row_sizes, bool remove_all,
+                        const int *col_idxs_map, size_t new_n_cols);
 
 void print_row_starts(const RowRange *row_ranges, size_t len);
 
@@ -120,6 +120,7 @@ Matrix *random_matrix_new(size_t n_rows, size_t n_cols, double density);
 // rows; otherwise it throws an assertion
 void replace_row_A(Matrix *A, int row, double ratio, double *new_vals, int *cols_new,
                    int new_len);
+
 #endif
 
 #endif // SPARSE_LA_MATRIX_H

include/core/Workspace.h

@@ -46,6 +46,10 @@ typedef struct Work
     // int_vec is used within parallel_rows and dual propagation
     iVec *int_vec;
 
+    // radix_aux is used as auxiliary buffer for radix sort in parallel_rows/cols
+    // TODO: do we really need this? or can we reuse space?
+    int *radix_aux;
+
     Mapping *mappings;
 } Work;
 

include/core/debug_macros.h

@@ -22,6 +22,21 @@
 #include "Tags.h"
 #include <stdbool.h>
 #include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/* Assertion that fires even in release mode (NDEBUG). Only use inside
+   the debugger (Debugger.c / Debugger.h) — never in production code. */
+#define PSLP_ASSERT(cond)                                                           \
+    do                                                                              \
+    {                                                                               \
+        if (!(cond))                                                                \
+        {                                                                           \
+            fprintf(stderr, "PSLP assertion failed: %s [%s:%d]\n", #cond, __FILE__, \
+                    __LINE__);                                                      \
+            abort();                                                                \
+        }                                                                           \
+    } while (0)
 
 #ifdef NDEBUG
 #define HUGE_BOUND_IS_OK
@@ -89,7 +104,7 @@ static inline bool ARRAYS_EQUAL_COLTAG(ColTag *arr1, const ColTag *arr2, int siz
     {                                                                               \
         for (int iii = 0; iii < (len); iii++)                                       \
         {                                                                           \
-            assert(!HAS_TAG(row_tags[rows[iii]], R_TAG_INACTIVE));                  \
+            PSLP_ASSERT(!HAS_TAG(row_tags[rows[iii]], R_TAG_INACTIVE));             \
         }                                                                           \
     } while (0)
 

include/core/radix_sort.h

@@ -0,0 +1,34 @@
+/*
+ * Copyright 2025 Daniel Cederberg
+ *
+ * This file is part of the PSLP project (LP Presolver).
+ *
+ * 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.
+ */
+
+#ifndef RADIX_SORT_H
+#define RADIX_SORT_H
+
+#include <stddef.h>
+
+// LSD radix sort on composite key (sparsity_ID, coeff_hash).
+// Sorts row indices in-place. aux must have space for n ints.
+void radix_sort_rows(int *rows, size_t n, const int *sparsity_IDs,
+                     const int *coeff_hashes, int *aux);
+
+// Parallel version: splits into 4 chunks, sorts in parallel, merges.
+// Falls back to sequential radix_sort_rows for n < 100000.
+void parallel_radix_sort_rows(int *rows, size_t n, const int *sparsity_IDs,
+                              const int *coeff_hashes, int *aux);
+
+#endif /* RADIX_SORT_H */

include/data_structures/Binary_search.h

@@ -0,0 +1,92 @@
+/*
+ * Copyright 2025 Daniel Cederberg
+ *
+ * This file is part of the PSLP project (LP Presolver).
+ *
+ * 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.
+ */
+
+#ifndef CORE_BINARY_SEARCH_H
+#define CORE_BINARY_SEARCH_H
+
+#define BSEARCH_LINEAR_THRESHOLD 8
+
+/* Find the index of 'target' in a sorted array 'arr' of length 'len'.
+   Returns the relative index (0-based) if found, or -1 if not found.
+   Falls back to linear scan for short arrays. */
+static inline int sorted_find(const int *arr, int len, int target)
+{
+    if (len <= BSEARCH_LINEAR_THRESHOLD)
+    {
+        for (int i = 0; i < len; ++i)
+        {
+            if (arr[i] == target) return i;
+            if (arr[i] > target) return -1;
+        }
+        return -1;
+    }
+
+    int lo = 0, hi = len - 1;
+    while (lo <= hi)
+    {
+        int mid = lo + (hi - lo) / 2;
+        if (arr[mid] == target)
+        {
+            return mid;
+        }
+        if (arr[mid] < target)
+        {
+            lo = mid + 1;
+        }
+        else
+        {
+            hi = mid - 1;
+        }
+    }
+    return -1;
+}
+
+/* Find the first index where arr[i] >= target in a sorted array of length 'len'.
+   Returns 'len' if all elements are less than target.
+   Falls back to linear scan for short arrays. */
+static inline int sorted_lower_bound(const int *arr, int len, int target)
+{
+    if (len <= BSEARCH_LINEAR_THRESHOLD)
+    {
+        for (int i = 0; i < len; ++i)
+        {
+            if (arr[i] >= target)
+            {
+                return i;
+            }
+        }
+        return len;
+    }
+
+    int lo = 0, hi = len;
+    while (lo < hi)
+    {
+        int mid = lo + (hi - lo) / 2;
+        if (arr[mid] < target)
+        {
+            lo = mid + 1;
+        }
+        else
+        {
+            hi = mid;
+        }
+    }
+    return lo;
+}
+
+#endif // CORE_BINARY_SEARCH_H

include/explorers/Parallel_cols.h

@@ -25,7 +25,7 @@
 struct Problem;
 
 /* This function finds parallel columns and merges them. It always
-   returns UNCHANGED. When this function finishes, the problem is
+   returns UNCHANGED or UNBNDORINFEAS. When this function finishes, the problem is
    up to date, and no extra synchronization is needed.
 */
 PresolveStatus remove_parallel_cols(struct Problem *prob);

include/explorers/Parallel_rows.h

@@ -55,7 +55,8 @@ void compute_supp_and_coeff_hash(const struct Matrix *A, const RowTag *rowtags,
 */
 void find_parallel_rows(const struct Matrix *A, const RowTag *r_Tags,
                         struct iVec *group_starts, int *parallel_rows,
-                        int *sparsity_IDs, int *coeff_hashes, RowTag INACTIVE_TAG);
+                        int *sparsity_IDs, int *coeff_hashes, RowTag INACTIVE_TAG,
+                        int *radix_aux);
 
 /* This function finds parallel rows and deletes the redundant ones, or declares
    the problem as infeasible if the parallel rows are not consistent. It returns

src/core/Activity.c

@@ -484,8 +484,7 @@ Altered_Activity update_act_bound_change(Activity *act, double coeff,
 void update_activities_bound_change(Activity *activities, const Matrix *A,
                                     const Bound *bounds, const double *vals,
                                     const int *rows, int len, double old_bound,
-                                    double new_bound, double finite_bound,
-                                    bool lower,
+                                    double new_bound, bool finite_bound, bool lower,
                                     iVec *updated_activities HUGE_BOUND_PARAM)
 {
     assert(!IS_HUGE(new_bound) || huge_bound_ok);
@@ -705,15 +704,22 @@ Altered_Activity update_activity_coeff_change(Activity *act, double lb, double u
 #endif
 
     // recompute from scratch if necessary
-    if (act->n_inf_min == 0 && n_inf_min_before == 1)
+    bool recompute_min = (act->n_inf_min == 0 && n_inf_min_before == 1);
+    bool recompute_max = (act->n_inf_max == 0 && n_inf_max_before == 1);
+
+    if (recompute_max && recompute_min)
     {
-        return MIN_ALTERED_RECOMPUTE;
+        assert(new_coeff == 0.0);
+        return MIN_ALTERED_RECOMPUTE | MAX_ALTERED_RECOMPUTE;
     }
-
-    if (act->n_inf_max == 0 && n_inf_max_before == 1)
+    else if (recompute_max)
     {
         return MAX_ALTERED_RECOMPUTE;
     }
+    else if (recompute_min)
+    {
+        return MIN_ALTERED_RECOMPUTE;
+    }
 
     return NO_RECOMPUTE;
 }

src/core/Constraints.c

@@ -302,18 +302,45 @@ static void bounds_shrink(Bound *ptr, int *map, size_t len)
     }
 }
 
+size_t update_column_map(const int *col_sizes, int *map, size_t n_cols)
+{
+    size_t col_count = 0;
+    for (size_t col = 0; col < n_cols; ++col)
+    {
+        if (col_sizes[col] == SIZE_INACTIVE_COL)
+        {
+            map[col] = -1;
+        }
+        else
+        {
+            map[col] = (int) col_count++;
+        }
+    }
+
+    return col_count;
+}
+
 void constraints_clean(Constraints *constraints, Mapping *map, bool remove_all)
 {
     const int *row_sizes = constraints->state->row_sizes;
     const int *col_sizes = constraints->state->col_sizes;
 
-    remove_extra_space(constraints->A, row_sizes, col_sizes, remove_all, map->cols);
-    remove_extra_space(constraints->AT, col_sizes, row_sizes, remove_all, map->rows);
+    size_t new_n_cols = update_column_map(col_sizes, map->cols, constraints->n);
+    remove_extra_space(constraints->A, row_sizes, remove_all, map->cols, new_n_cols);
 
-    dPtr_shrink(constraints->rhs, map->rows, constraints->m);
-    dPtr_shrink(constraints->lhs, map->rows, constraints->m);
+    size_t new_n_rows = update_column_map(row_sizes, map->rows, constraints->m);
+    (void) new_n_rows;
+
+    /* we don't need these for anything so no need to clean it */
+#if defined(TESTING) || !defined(NDEBUG)
+    remove_extra_space(constraints->AT, col_sizes, remove_all, map->rows,
+                       new_n_rows);
     rowTagPtr_shrink(constraints->row_tags, map->rows, constraints->m);
     colTagPtr_shrink(constraints->col_tags, map->cols, constraints->n);
+#endif
+
+    dPtr_shrink(constraints->rhs, map->rows, constraints->m);
+    dPtr_shrink(constraints->lhs, map->rows, constraints->m);
     bounds_shrink(constraints->bounds, map->cols, constraints->n);
     constraints->m = constraints->A->m;
     constraints->n = constraints->A->n;