• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

dance858 / PSLP / 21523923575

30 Jan 2026 05:03PM UTC coverage: 89.013% (+0.03%) from 88.988%
21523923575

Pull #34

github

web-flow
Merge 30140dc5c into 8137ce2f7
Pull Request #34: fix-windows-build

1249 of 1400 branches covered (89.21%)

Branch coverage included in aggregate %.

222 of 232 new or added lines in 20 files covered. (95.69%)

1 existing line in 1 file now uncovered.

3855 of 4334 relevant lines covered (88.95%)

3375.4 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.26
/src/core/Presolver.c
1
/*
2
 * Copyright 2025 Daniel Cederberg
3
 *
4
 * This file is part of the PSLP project (LP Presolver).
5
 *
6
 * Licensed under the Apache License, Version 2.0 (the "License");
7
 * you may not use this file except in compliance with the License.
8
 * You may obtain a copy of the License at
9
 *
10
 *     http://www.apache.org/licenses/LICENSE-2.0
11
 *
12
 * Unless required by applicable law or agreed to in writing, software
13
 * distributed under the License is distributed on an "AS IS" BASIS,
14
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15
 * See the License for the specific language governing permissions and
16
 * limitations under the License.
17
 */
18

19
#include "Activity.h"
20
#include "Constraints.h"
21
#include "CoreTransformations.h"
22
#include "DTonsEq.h"
23
#include "Debugger.h"
24
#include "Locks.h"
25
#include "Matrix.h"
26
#include "Memory_wrapper.h"
27
#include "Numerics.h"
28
#include "PSLP_API.h"
29
#include "PSLP_sol.h"
30
#include "PSLP_stats.h"
31
#include "PSLP_status.h"
32
#include "Parallel_cols.h"
33
#include "Parallel_rows.h"
34
#include "Postsolver.h"
35
#include "Primal_propagation.h"
36
#include "Problem.h"
37
#include "RowColViews.h"
38
#include "SimpleReductions.h"
39
#include "Simple_dual_fix.h"
40
#include "State.h"
41
#include "StonCols.h"
42
#include "Tags.h"
43
#include "Timer.h"
44
#include "Workspace.h"
45
#include "dVec.h"
46
#include "glbopts.h"
47
#include "iVec.h"
48
#include "pslp_thread.h"
49
#include <stdio.h>
50
#include <stdlib.h>
51

52
#define NNZ_REMOVED_FAST 0.95
53
#define NNZ_REMOVED_CYCLE 0.95
54

55
typedef uint8_t Complexity;
56

57
enum Complexity
58
{
59
    FAST = 0,
60
    MEDIUM = 1 << 0
61
};
62

63
PresolveStats *init_stats(size_t n_rows, size_t n_cols, size_t nnz)
99✔
64
{
65
    PresolveStats *stats = (PresolveStats *) ps_calloc(1, sizeof(PresolveStats));
99✔
66
    stats->n_rows_original = n_rows;
99✔
67
    stats->n_cols_original = n_cols;
99✔
68
    stats->nnz_original = nnz;
99✔
69
    return stats;
99✔
70
}
71

72
void free_settings(Settings *stgs)
8✔
73
{
74
    PS_FREE(stgs);
8✔
75
}
8✔
76

77
Settings *default_settings()
103✔
78
{
79
    Settings *stgs = (Settings *) ps_malloc(1, sizeof(Settings));
103✔
80
    RETURN_PTR_IF_NULL(stgs, NULL);
103✔
81
    stgs->ston_cols = true;
103✔
82
    stgs->dton_eq = true;
103✔
83
    stgs->parallel_rows = true;
103✔
84
    stgs->parallel_cols = true;
103✔
85
    stgs->primal_propagation = true;
103✔
86
    stgs->dual_fix = true;
103✔
87
    // stgs->clean_small_coeff = false;
88
    stgs->finite_bound_tightening = true;
103✔
89
    stgs->relax_bounds = true;
103✔
90
    stgs->max_shift = 10;
103✔
91
    stgs->max_time = 60.0;
103✔
92
    stgs->verbose = true;
103✔
93
    return stgs;
103✔
94
}
95

96
// Helper struct for parallel initialization work
97
typedef struct
98
{
99
    Matrix *A;
100
    Work *work;
101
    size_t n_cols;
102
    size_t n_rows;
103
    const double *lbs;
104
    const double *ubs;
105
    double *lhs_copy;
106
    double *rhs_copy;
107
    Bound *bounds;
108
    ColTag *col_tags;
109
    RowTag *row_tags;
110
    Lock *locks;
111
    Activity *activities;
112
    int *row_sizes;
113
} ParallelInitData;
114

115
// thread function for initializing data
116
static void *init_thread_func(void *arg)
99✔
117
{
118
    ParallelInitData *data = (ParallelInitData *) arg;
99✔
119

120
    data->row_tags = new_rowtags(data->lhs_copy, data->rhs_copy, data->n_rows);
99✔
121

122
    for (int i = 0; i < data->n_cols; i++)
573✔
123
    {
124
        data->bounds[i].lb = data->lbs[i];
474✔
125
        data->bounds[i].ub = data->ubs[i];
474✔
126

127
        if (IS_NEG_INF(data->lbs[i]))
474✔
128
        {
129
            UPDATE_TAG(data->col_tags[i], C_TAG_LB_INF);
49✔
130
        }
131

132
        if (IS_POS_INF(data->ubs[i]))
474✔
133
        {
134

135
            UPDATE_TAG(data->col_tags[i], C_TAG_UB_INF);
287✔
136
        }
137
    }
138

139
    data->locks = new_locks(data->A, data->row_tags);
99✔
140
    data->activities = new_activities(data->A, data->col_tags, data->bounds);
99✔
141
    count_rows(data->A, data->row_sizes);
99✔
142

143
    return NULL;
99✔
144
}
145

146
void set_settings_true(Settings *stgs)
25✔
147
{
148
    stgs->ston_cols = true;
25✔
149
    stgs->dton_eq = true;
25✔
150
    stgs->parallel_rows = true;
25✔
151
    stgs->parallel_cols = true;
25✔
152
    stgs->primal_propagation = true;
25✔
153
    stgs->dual_fix = true;
25✔
154
    // stgs->clean_small_coeff = true;
155
    stgs->finite_bound_tightening = true;
25✔
156
    stgs->relax_bounds = true;
25✔
157
    stgs->verbose = true;
25✔
158
}
25✔
159

160
void set_settings_false(Settings *stgs)
17✔
161
{
162
    stgs->ston_cols = false;
17✔
163
    stgs->dton_eq = false;
17✔
164
    stgs->parallel_rows = false;
17✔
165
    stgs->parallel_cols = false;
17✔
166
    stgs->primal_propagation = false;
17✔
167
    stgs->dual_fix = false;
17✔
168
    // stgs->clean_small_coeff = false;
169
    stgs->finite_bound_tightening = false;
17✔
170
    stgs->relax_bounds = false;
17✔
171
    stgs->verbose = false;
17✔
172
}
17✔
173

174
typedef struct clean_up_scope
175
{
176
    Matrix *A, *AT;
177
    double *lhs_copy, *rhs_copy, *c_copy;
178
    int *row_sizes, *col_sizes;
179
    RowTag *row_tags;
180
    Lock *locks;
181
    Activity *activities;
182
    State *data;
183
    Constraints *constraints;
184
    Presolver *presolver;
185
    Objective *obj;
186
    ColTag *col_tags;
187
    Bound *bounds;
188
    Work *work;
189
} clean_up_scope;
190

191
void presolver_clean_up(clean_up_scope scope)
×
192
{
193
    free_matrix(scope.A);
×
194
    free_matrix(scope.AT);
×
195
    PS_FREE(scope.lhs_copy);
×
196
    PS_FREE(scope.rhs_copy);
×
197
    PS_FREE(scope.c_copy);
×
198
    PS_FREE(scope.row_sizes);
×
199
    PS_FREE(scope.col_sizes);
×
200
    PS_FREE(scope.row_tags);
×
201
    PS_FREE(scope.col_tags);
×
202
    PS_FREE(scope.bounds);
×
203
    free_activities(scope.activities);
×
204
    free_locks(scope.locks);
×
205
    free_work(scope.work);
×
206
    free_state(scope.data);
×
207
    free_constraints(scope.constraints);
×
208
    free_objective(scope.obj);
×
209
    free_presolver(scope.presolver);
×
210
}
×
211

212
Presolver *new_presolver(const double *Ax, const int *Ai, const int *Ap, size_t m,
99✔
213
                         size_t n, size_t nnz, const double *lhs, const double *rhs,
214
                         const double *lbs, const double *ubs, const double *c,
215
                         const Settings *stgs)
216
{
217
    Timer timer;
218
    clock_gettime(CLOCK_MONOTONIC, &timer.start);
99✔
219
    Matrix *A = NULL, *AT = NULL;
99✔
220
    double *lhs_copy = NULL, *rhs_copy = NULL, *c_copy = NULL;
99✔
221
    int *row_sizes = NULL, *col_sizes = NULL;
99✔
222
    RowTag *row_tags = NULL;
99✔
223
    Lock *locks = NULL;
99✔
224
    Activity *activities = NULL;
99✔
225
    State *data = NULL;
99✔
226
    Constraints *constraints = NULL;
99✔
227
    Presolver *presolver = NULL;
99✔
228
    Objective *obj = NULL;
99✔
229
    ColTag *col_tags = NULL;
99✔
230
    Bound *bounds = NULL;
99✔
231
    Work *work = NULL;
99✔
232

233
    //  ---------------------------------------------------------------------------
234
    //   Copy data and allocate memory. For an exact specification of the
235
    //   workspace, see the workspace class. The problem object owns all the
236
    //   memory that is allocated in this block of code (and therefore frees it).
237
    //  ---------------------------------------------------------------------------
238
    size_t n_rows = m;
99✔
239
    size_t n_cols = n;
99✔
240
    lhs_copy = (double *) ps_malloc(n_rows, sizeof(double));
99✔
241
    rhs_copy = (double *) ps_malloc(n_rows, sizeof(double));
99✔
242
    c_copy = (double *) ps_malloc(n_cols, sizeof(double));
99✔
243
    col_tags = (ColTag *) ps_calloc(n_cols, sizeof(ColTag));
99✔
244
    bounds = (Bound *) ps_malloc(n_cols, sizeof(Bound));
99✔
245
    work = new_work(n_rows, n_cols);
99✔
246
    row_sizes = (int *) ps_malloc(n_rows, sizeof(int));
99✔
247
    col_sizes = (int *) ps_malloc(n_cols, sizeof(int));
99✔
248

249
    if (!lhs_copy || !rhs_copy || !c_copy || !col_tags || !bounds || !work ||
99✔
250
        !row_sizes || !col_sizes)
99✔
251
    {
252
        goto cleanup;
×
253
    }
254

255
    memcpy(lhs_copy, lhs, n_rows * sizeof(double));
99✔
256
    memcpy(rhs_copy, rhs, n_rows * sizeof(double));
99✔
257
    memcpy(c_copy, c, n_cols * sizeof(double));
99✔
258

259
    // ---------------------------------------------------------------------------
260
    //  Build bounds, row tags, A and AT. We first build A, then in parallel
261
    //  we build AT (main thread) and some other things (second thread).
262
    // ---------------------------------------------------------------------------
263
    A = matrix_new_no_extra_space(Ax, Ai, Ap, n_rows, n_cols, nnz);
99✔
264
    if (!A) goto cleanup;
99✔
265

266
    // commented out for now because does not seem helpful
267
    // if (stgs->clean_small_coeff)
268
    // {
269
    //     clean_small_coeff_A(A, bounds, row_tags, col_tags, rhs_copy, lhs_copy);
270
    // }
271

272
    ps_thread_t thread_id;
273
    ParallelInitData parallel_data = {A,    work,     n_cols,   n_rows,   lbs,
99✔
274
                                      ubs,  lhs_copy, rhs_copy, bounds,   col_tags,
275
                                      NULL, NULL,     NULL,     row_sizes};
276

277
    ps_thread_create(&thread_id, NULL, init_thread_func, &parallel_data);
99✔
278

279
    // Main thread: Transpose A and count rows
280
    AT = transpose(A, work->iwork_n_cols);
99✔
281
    if (!AT)
99✔
282
    {
NEW
283
        ps_thread_join(&thread_id, NULL);
×
284
        goto cleanup;
×
285
    }
286
    count_rows(AT, col_sizes);
99✔
287

288
    // sync threads
289
    ps_thread_join(&thread_id, NULL);
99✔
290

291
    row_tags = parallel_data.row_tags;
99✔
292
    if (!row_tags) goto cleanup;
99✔
293
    locks = parallel_data.locks;
99✔
294
    activities = parallel_data.activities;
99✔
295
    if (!locks || !activities) goto cleanup;
99✔
296

297
    // ---------------------------------------------------------------------------
298
    //  Initialize internal data and constraints
299
    // ---------------------------------------------------------------------------
300
    data = new_state(row_sizes, col_sizes, locks, n_rows, n_cols, activities, work,
99✔
301
                     row_tags);
302

303
    if (!data) goto cleanup;
99✔
304
    constraints =
305
        constraints_new(A, AT, lhs_copy, rhs_copy, bounds, data, row_tags, col_tags);
99✔
306
    if (!constraints) goto cleanup;
99✔
307

308
    // ---------------------------------------------------------------------------
309
    //             Allocate the actual presolver
310
    // ---------------------------------------------------------------------------
311
    presolver = (Presolver *) ps_malloc(1, sizeof(Presolver));
99✔
312
    obj = objective_new(c_copy);
99✔
313
    if (!presolver || !obj) goto cleanup;
99✔
314
    presolver->stgs = stgs;
99✔
315
    presolver->prob = new_problem(constraints, obj);
99✔
316
    presolver->stats = init_stats(A->m, A->n, nnz);
99✔
317
    // presolver->stats->nnz_removed_trivial = nnz - A->nnz; due to clean_small_coeff
318
    presolver->reduced_prob =
99✔
319
        (PresolvedProblem *) ps_calloc(1, sizeof(PresolvedProblem));
99✔
320
    DEBUG(run_debugger(constraints, false));
99✔
321

322
    // ---------------------------------------------------------------------------
323
    //           Allocate space for returning the solution
324
    // ---------------------------------------------------------------------------
325
    presolver->sol = (Solution *) ps_malloc(1, sizeof(Solution));
99✔
326
    if (!presolver->sol) goto cleanup;
99✔
327
    presolver->sol->x = (double *) ps_malloc(n_cols, sizeof(double));
99✔
328
    presolver->sol->y = (double *) ps_malloc(n_rows, sizeof(double));
99✔
329
    presolver->sol->z = (double *) ps_malloc(n_cols, sizeof(double));
99✔
330
    presolver->sol->dim_x = n_cols;
99✔
331
    presolver->sol->dim_y = n_rows;
99✔
332
    if (!presolver->sol->x || !presolver->sol->y || !presolver->sol->z)
99✔
333
    {
334
        goto cleanup;
×
335
    }
336

337
    clock_gettime(CLOCK_MONOTONIC, &timer.end);
99✔
338
    presolver->stats->time_init = GET_ELAPSED_SECONDS(timer);
99✔
339

340
    return presolver;
99✔
341

342
cleanup:
×
343
{
344
    struct clean_up_scope scope = {
×
345
        A,         AT,       lhs_copy, rhs_copy,   c_copy, row_sizes,
346
        col_sizes, row_tags, locks,    activities, data,   constraints,
347
        presolver, obj,      col_tags, bounds,     work};
348
    presolver_clean_up(scope);
×
349
}
350
    return NULL;
×
351
}
352

353
/* This function updates the termination criterion for the presolver.
354
   We terminate if one of the following conditions is satisfied:
355
   1. The problem is infeasible or unbounded.
356
   2. A cycle reduced the number of non-zeros by less than
357
     (1 - NNZ_REMOVED_CYCLE)%.
358
   3. A maximum time limit is reached.
359
*/
360
static inline bool update_termination(size_t nnz_after_cycle,
87✔
361
                                      size_t nnz_before_cycle, Complexity complexity,
362
                                      PresolveStatus status, double max_time,
363
                                      Timer outer_timer)
364
{
365
    if (HAS_STATUS(status, UNBNDORINFEAS))
87✔
366
    {
367
        return true;
2✔
368
    }
369

370
    if (complexity == MEDIUM &&
85✔
371
        (double) nnz_after_cycle >= NNZ_REMOVED_CYCLE * (double) nnz_before_cycle)
38✔
372
    {
373
        return true;
20✔
374
    }
375

376
    clock_gettime(CLOCK_MONOTONIC, &outer_timer.end);
65✔
377

378
    if (GET_ELAPSED_SECONDS(outer_timer) >= max_time)
65✔
379
    {
380
        printf("Maximum time limit of %.2f seconds reached.\n", max_time);
×
381
        return true;
×
382
    }
383

384
    return false;
65✔
385
}
386

387
/* This function updates the complexity class according to
388
   the following rules:
389
    * If the current complexity is fast, the next complexity is
390
      fast if sufficiently many non-zeros were removed, otherwise
391
      it is medium.
392
    * If the current complexity is medium, the next complexity is
393
      fast.
394
*/
395
static inline Complexity update_complexity(Complexity curr_complexity,
87✔
396
                                           size_t nnz_before_phase,
397
                                           size_t nnz_after_phase)
398
{
399
    if (curr_complexity == FAST)
87✔
400
    {
401
        if ((double) nnz_after_phase < NNZ_REMOVED_FAST * (double) nnz_before_phase)
48✔
402
        {
403
            return FAST;
9✔
404
        }
405

406
        return MEDIUM;
39✔
407
    }
408
    else if (curr_complexity == MEDIUM)
39✔
409
    {
410
        return FAST;
39✔
411
    }
412

413
    // should never reach here
414
    assert(false);
×
415
    return FAST;
416
}
417

418
/* This function exhaustively removes singleton rows from the problem. It
419
   also eliminates empty rows and empty columns. If the problem is feasible
420
   and bounded, it returns UNCHANGED (even if the problem was reduced).
421
*/
422
static inline PresolveStatus run_trivial_explorers(Problem *prob,
211✔
423
                                                   const Settings *stgs)
424
{
425
    PresolveStatus status = UNCHANGED;
211✔
426

427
    // TODO: should we just run this once in the beginning of presolve?
428
    // Very important question.
429
    status = remove_variables_with_close_bounds(prob);
211✔
430
    RETURN_IF_INFEASIBLE(status);
211✔
431

432
    if (stgs->dual_fix)
209✔
433
    {
434
        //   after simple dual fix there can be new empty rows and singleton rows,
435
        //   and empty columns (if rows are marked as inactive due to a variable
436
        //   being set to inf)
437
        status |= remove_empty_cols(prob);
157✔
438
        status |= simple_dual_fix(prob);
157✔
439
        RETURN_IF_UNBNDORINFEAS(status);
157✔
440
    }
441

442
    do
443
    {
444
        status = remove_ston_rows(prob);
217✔
445
        RETURN_IF_INFEASIBLE(status);
217✔
446
    } while (status == REDUCED);
216✔
447

448
    status = remove_empty_rows(prob->constraints);
206✔
449
    RETURN_IF_INFEASIBLE(status);
206✔
450

451
    status = remove_empty_cols(prob);
206✔
452
    RETURN_IF_UNBNDORINFEAS(status);
206✔
453

454
    assert(prob->constraints->state->ston_rows->len == 0);
206✔
455
    assert(prob->constraints->state->empty_rows->len == 0);
206✔
456
    assert(prob->constraints->state->empty_cols->len == 0);
206✔
457
    return UNCHANGED;
206✔
458
}
459

460
/* This function applies the fastest reductions, which are doubleton equality
461
   rows, singleton columns, and simple dual fix. It returns UNCHANGED even if
462
   the problem was reduced (provided that the problem does not seem infeasible or
463
   bounded). */
464
static inline PresolveStatus run_fast_explorers(Problem *prob, const Settings *stgs)
48✔
465
{
466
    assert(prob->constraints->state->ston_rows->len == 0);
48✔
467
    assert(prob->constraints->state->empty_rows->len == 0);
48✔
468
    assert(prob->constraints->state->empty_cols->len == 0);
48✔
469
    PresolveStatus status = UNCHANGED;
48✔
470

471
    if (stgs->ston_cols)
48✔
472
    {
473
        status |= remove_ston_cols(prob);
45✔
474

475
        // after removing singleton columns, there can be new empty columns
476
        // and singleton rows, but no empty rows
477
        assert(prob->constraints->state->empty_rows->len == 0);
45✔
478
        status |= run_trivial_explorers(prob, stgs);
45✔
479
    }
480

481
    if (stgs->dton_eq)
48✔
482
    {
483
        status |= remove_dton_eq_rows(prob, stgs->max_shift);
48✔
484
        // after removing doubleton equality rows, there can be new empty rows,
485
        // new singleton rows, and new empty columns
486
        status |= run_trivial_explorers(prob, stgs);
48✔
487
    }
488
    return status;
48✔
489
}
490

491
/* This function applies the medium complexity presolvers, which are domain
492
   propagation and parallel rows. It returns UNCHANGED even if the problem
493
   was reduced (provided that the problem is infeasible and bounded). */
494
static inline PresolveStatus
495
run_medium_explorers(Problem *prob, const Settings *stgs, PresolveStats *stats)
39✔
496
{
497
    assert(prob->constraints->state->ston_rows->len == 0);
39✔
498
    assert(prob->constraints->state->empty_rows->len == 0);
39✔
499
    assert(prob->constraints->state->empty_cols->len == 0);
39✔
500
    DEBUG(verify_no_duplicates_sort(prob->constraints->state->updated_activities));
39✔
501
    size_t nnz_before;
502
    size_t *nnz = &prob->constraints->A->nnz;
39✔
503
    PresolveStatus status = UNCHANGED;
39✔
504
    Timer timer;
505

506
    if (stgs->primal_propagation)
39✔
507
    {
508
        // stats
509
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
15✔
510
        nnz_before = *nnz;
15✔
511

512
        // the actual reduction
513
        status |= check_activities(prob);
15✔
514
        status |= propagate_primal(prob, stgs->finite_bound_tightening);
15✔
515

516
        // after dom prop propagation there can be new empty and ston rows
517
        status |= run_trivial_explorers(prob, stgs);
15✔
518
        assert(!(status & REDUCED));
15✔
519
        RETURN_IF_NOT_UNCHANGED(status);
15✔
520

521
        // stats
522
        stats->nnz_removed_primal_propagation += nnz_before - *nnz;
14✔
523
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
14✔
524
        stats->time_primal_propagation += GET_ELAPSED_SECONDS(timer);
14✔
525
    }
526

527
    if (stgs->parallel_rows)
38✔
528
    {
529
        // stats
530
        nnz_before = *nnz;
38✔
531
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
38✔
532

533
        // the actual reduction (after removing parallel rows there can
534
        // be no new empty or ston rows so no need to run trivial presolvers)
535
        status |= remove_parallel_rows(prob->constraints);
38✔
536
        assert(prob->constraints->state->empty_rows->len == 0);
38✔
537
        assert(prob->constraints->state->ston_rows->len == 0);
38✔
538
        assert(!(status & REDUCED));
38✔
539
        RETURN_IF_NOT_UNCHANGED(status);
38✔
540

541
        // stats
542
        stats->nnz_removed_parallel_rows += nnz_before - *nnz;
38✔
543
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
38✔
544
        stats->time_parallel_rows += GET_ELAPSED_SECONDS(timer);
38✔
545
    }
546

547
    if (stgs->parallel_cols)
38✔
548
    {
549
        // stats
550
        nnz_before = *nnz;
12✔
551
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
12✔
552

553
        // the actual reduction (after removing parallel columns there can
554
        // be no new empty rows or empty cols but might be new stonrows, probably
555
        // although that's rare but it does happen
556
        status |= remove_parallel_cols(prob);
12✔
557
        assert(prob->constraints->state->empty_rows->len == 0);
12✔
558
        assert(prob->constraints->state->empty_cols->len == 0);
12✔
559
        status |= run_trivial_explorers(prob, stgs);
12✔
560
        assert(prob->constraints->state->ston_rows->len == 0);
12✔
561
        assert(!(status & REDUCED));
12✔
562
        RETURN_IF_NOT_UNCHANGED(status);
12✔
563

564
        // stats
565
        stats->nnz_removed_parallel_cols += nnz_before - *nnz;
12✔
566
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
12✔
567
        stats->time_parallel_cols += GET_ELAPSED_SECONDS(timer);
12✔
568
    }
569

570
    return status;
38✔
571
}
572

573
void populate_presolved_problem(Presolver *presolver)
20✔
574
{
575
    PresolvedProblem *reduced_prob = presolver->reduced_prob;
20✔
576
    Constraints *constraints = presolver->prob->constraints;
20✔
577
    Matrix *A = constraints->A;
20✔
578
    reduced_prob->m = A->m;
20✔
579
    reduced_prob->n = A->n;
20✔
580
    reduced_prob->nnz = A->nnz;
20✔
581
    reduced_prob->Ax = A->x;
20✔
582
    reduced_prob->Ai = A->i;
20✔
583
    reduced_prob->rhs = constraints->rhs;
20✔
584
    reduced_prob->lhs = constraints->lhs;
20✔
585
    reduced_prob->c = presolver->prob->obj->c;
20✔
586
    reduced_prob->obj_offset = presolver->prob->obj->offset;
20✔
587

588
    // create bounds arrays
589
    reduced_prob->lbs = (double *) malloc(A->n * sizeof(double));
20✔
590
    reduced_prob->ubs = (double *) malloc(A->n * sizeof(double));
20✔
591
    for (int i = 0; i < A->n; i++)
118✔
592
    {
593
        reduced_prob->lbs[i] = constraints->bounds[i].lb;
98✔
594
        reduced_prob->ubs[i] = constraints->bounds[i].ub;
98✔
595
    }
596

597
    // create row pointers
598
    reduced_prob->Ap = (int *) malloc((A->m + 1) * sizeof(int));
20✔
599
    for (int i = 0; i < A->m + 1; i++)
103✔
600
    {
601
        reduced_prob->Ap[i] = A->p[i].start;
83✔
602
    }
603
}
20✔
604

605
static inline void print_start_message(const PresolveStats *stats)
26✔
606
{
607
    printf("\n\t       PSLP v%s - LP presolver \n\t(c) Daniel "
26✔
608
           "Cederberg, Stanford University, 2025\n",
609
           PSLP_VERSION);
610
    printf("Original problem:  %zu rows, %zu columns, %zu nnz\n",
26✔
611
           stats->n_rows_original, stats->n_cols_original, stats->nnz_original);
26✔
612
}
26✔
613

614
static inline void print_end_message(const Matrix *A, const PresolveStats *stats)
20✔
615
{
616
    printf("Presolved problem: %zu rows, %zu columns, %zu nnz\n", A->m, A->n,
20✔
617
           A->nnz);
20✔
618
    printf("PSLP init & run time : %.3f seconds, %.3f \n", stats->time_init,
20✔
619
           stats->time_presolve);
20✔
620
}
20✔
621

622
static inline void print_infeas_or_unbnd_message(PresolveStatus status)
6✔
623
{
624
    if (status == INFEASIBLE)
6✔
625
    {
626
        printf("PSLP declares problem as infeasible.\n");
4✔
627
        return;
4✔
628
    }
629
    else if (status == UNBNDORINFEAS)
2✔
630
    {
631
        printf("PSLP declares problem as infeasible or unbounded.\n");
2✔
632
        return;
2✔
633
    }
634

635
    // should never reach here
636
    assert(false);
×
637
}
638

639
static inline PresolveStatus final_status(PresolveStats *stats)
20✔
640
{
641
    if (stats->n_cols_original == stats->n_cols_reduced &&
20✔
642
        stats->n_rows_original == stats->n_rows_reduced &&
2✔
643
        stats->nnz_original == stats->nnz_reduced)
×
644
    {
645
        return UNCHANGED;
×
646
    }
647
    else
648
    {
649
        return REDUCED;
20✔
650
    }
651
}
652

653
PresolveStatus run_presolver(Presolver *presolver)
26✔
654
{
655
    Timer inner_timer, outer_timer;
656
    size_t nnz_before_cycle, nnz_after_cycle;
657
    size_t nnz_before_phase, nnz_after_phase;
658
    size_t nnz_before_reduction;
659
    Problem *prob = presolver->prob;
26✔
660
    PresolveStats *stats = presolver->stats;
26✔
661
    Matrix *A = presolver->prob->constraints->A;
26✔
662
    const Settings *stgs = presolver->stgs;
26✔
663
    Complexity curr_complexity = FAST;
26✔
664
    bool terminate = false;
26✔
665
    PresolveStatus status = UNCHANGED;
26✔
666
    clock_gettime(CLOCK_MONOTONIC, &outer_timer.start);
26✔
667

668
    if (stgs->verbose)
26✔
669
    {
670
        print_start_message(stats);
26✔
671
    }
672

673
    DEBUG(run_debugger(prob->constraints, false));
26✔
674

675
    // ------------------------------------------------------------------------
676
    // Main loop for the presolver. The logic is organized into phases and
677
    // cycles. A phase is a sequence of presolvers belonging to a certain
678
    // complexity class. The cycle resets after the medium presolvers.
679
    // ------------------------------------------------------------------------
680
    nnz_before_cycle = A->nnz;
26✔
681
    nnz_after_cycle = A->nnz; /* to avoid unitialized variable warning */
26✔
682
    while (!terminate)
113✔
683
    {
684
        // before each phase we run the trivial presolvers
685
        nnz_before_reduction = A->nnz;
91✔
686
        status = run_trivial_explorers(prob, stgs);
91✔
687

688
        // run_trivial_explorers returns something else than UNCHANGED only
689
        // if the problem is detected to be infeasible or unbounded
690
        if (status != UNCHANGED)
91✔
691
        {
692
            break;
4✔
693
        }
694

695
        stats->nnz_removed_trivial += nnz_before_reduction - A->nnz;
87✔
696
        DEBUG(run_debugger(prob->constraints, false));
87✔
697
        nnz_before_phase = A->nnz;
87✔
698

699
        // apply presolvers belonging to a certain complexity class
700
        if (curr_complexity == FAST)
87✔
701
        {
702
            nnz_before_reduction = A->nnz;
48✔
703
            RUN_AND_TIME(run_fast_explorers, inner_timer,
48✔
704
                         stats->time_fast_reductions, status, prob, stgs);
705
            stats->nnz_removed_fast += nnz_before_reduction - A->nnz;
48✔
706
        }
707
        else if (curr_complexity == MEDIUM)
39✔
708
        {
709
            RUN_AND_TIME(run_medium_explorers, inner_timer,
39✔
710
                         stats->time_medium_reductions, status, prob, stgs, stats);
711
            nnz_after_cycle = A->nnz;
39✔
712
        }
713

714
        nnz_after_phase = A->nnz;
87✔
715

716
        terminate =
717
            update_termination(nnz_after_cycle, nnz_before_cycle, curr_complexity,
87✔
718
                               status, stgs->max_time, outer_timer);
87✔
719

720
        if (curr_complexity == MEDIUM)
87✔
721
        {
722
            // a new cycle starts after the medium presolvers
723
            nnz_before_cycle = nnz_after_cycle;
39✔
724
        }
725

726
        curr_complexity =
727
            update_complexity(curr_complexity, nnz_before_phase, nnz_after_phase);
87✔
728
    }
729

730
    if (status != UNCHANGED)
26✔
731
    {
732
        // problem detected to be infeasible or unbounded
733
        print_infeas_or_unbnd_message(status);
6✔
734
        return status;
6✔
735
    }
736

737
    if (stgs->relax_bounds)
20✔
738
    {
739
        remove_redundant_bounds(prob->constraints);
20✔
740
    }
741

742
    problem_clean(prob, true);
20✔
743
    DEBUG(run_debugger(prob->constraints, true));
20✔
744
    stats->n_rows_reduced = A->m;
20✔
745
    stats->n_cols_reduced = A->n;
20✔
746
    stats->nnz_reduced = A->nnz;
20✔
747
    DEBUG(run_debugger_stats_consistency_check(stats));
20✔
748
    populate_presolved_problem(presolver);
20✔
749
    clock_gettime(CLOCK_MONOTONIC, &outer_timer.end);
20✔
750
    stats->time_presolve = GET_ELAPSED_SECONDS(outer_timer);
20✔
751

752
    if (stgs->verbose)
20✔
753
    {
754
        print_end_message(A, stats);
20✔
755
    }
756

757
    return final_status(stats);
20✔
758
}
759

760
void postsolve(Presolver *presolver, const double *x, const double *y,
18✔
761
               const double *z)
762
{
763
    Timer timer;
764
    Solution *sol = presolver->sol;
18✔
765
    PresolveStats *stats = presolver->stats;
18✔
766
    State *data = presolver->prob->constraints->state;
18✔
767
    PostsolveInfo *postsolve_info = data->postsolve_info;
18✔
768
    clock_gettime(CLOCK_MONOTONIC, &timer.start);
18✔
769
    postsolver_update(postsolve_info, stats->n_cols_reduced, stats->n_rows_reduced,
18✔
770
                      data->work->mappings->cols, data->work->mappings->rows);
18✔
771
    postsolver_run(postsolve_info, sol, x, y, z);
18✔
772
    clock_gettime(CLOCK_MONOTONIC, &timer.end);
18✔
773
    stats->time_postsolve = GET_ELAPSED_SECONDS(timer);
18✔
774

775
    if (presolver->stgs->verbose)
18✔
776
    {
777
        printf("PSLP postsolve time: %.4f seconds\n", stats->time_postsolve);
18✔
778
    }
779
}
18✔
780

781
void free_presolver(Presolver *presolver)
99✔
782
{
783
    if (presolver == NULL)
99✔
784
    {
785
        return;
×
786
    }
787

788
    free_problem(presolver->prob);
99✔
789
    PS_FREE(presolver->stats);
99✔
790

791
    if (presolver->reduced_prob)
99✔
792
    {
793
        PS_FREE(presolver->reduced_prob->Ap);
99✔
794
        PS_FREE(presolver->reduced_prob->lbs);
99✔
795
        PS_FREE(presolver->reduced_prob->ubs);
99✔
796
        PS_FREE(presolver->reduced_prob);
99✔
797
    }
798

799
    if (presolver->sol)
99✔
800
    {
801
        PS_FREE(presolver->sol->x);
99✔
802
        PS_FREE(presolver->sol->y);
99✔
803
        PS_FREE(presolver->sol->z);
99✔
804
        PS_FREE(presolver->sol);
99✔
805
    }
806

807
    PS_FREE(presolver);
99✔
808
}
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc