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

dance858 / PSLP / 20194334020

13 Dec 2025 03:49PM UTC coverage: 87.707% (-0.5%) from 88.244%
20194334020

Pull #22

github

web-flow
Merge 106c7e353 into 90ab253fc
Pull Request #22: Separate thread for transpose

1217 of 1383 branches covered (88.0%)

Branch coverage included in aggregate %.

26 of 28 new or added lines in 1 file covered. (92.86%)

23 existing lines in 1 file now uncovered.

3763 of 4295 relevant lines covered (87.61%)

3398.25 hits per line

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

91.53
/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 <pthread.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(int n_rows, int n_cols, int nnz)
89✔
64
{
65
    PresolveStats *stats = (PresolveStats *) ps_calloc(1, sizeof(PresolveStats));
89✔
66
    stats->n_rows_original = n_rows;
89✔
67
    stats->n_cols_original = n_cols;
89✔
68
    stats->nnz_original = nnz;
89✔
69
    return stats;
89✔
70
}
71

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

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

96
// Helper struct for parallel initialization work
97
typedef struct
98
{
99
    Matrix *A;
100
    Work *work;
101
    int n_cols;
102
    int 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)
89✔
117
{
118
    ParallelInitData *data = (ParallelInitData *) arg;
89✔
119

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

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

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

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

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

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

143
    return NULL;
89✔
144
}
145

146
void set_settings_true(Settings *stgs)
16✔
147
{
148
    stgs->ston_cols = true;
16✔
149
    stgs->dton_eq = true;
16✔
150
    stgs->parallel_rows = true;
16✔
151
    stgs->parallel_cols = true;
16✔
152
    stgs->primal_propagation = true;
16✔
153
    stgs->dual_fix = true;
16✔
154
    // stgs->clean_small_coeff = true;
155
    stgs->finite_bound_tightening = true;
16✔
156
    stgs->relax_bounds = true;
16✔
157
    stgs->verbose = true;
16✔
158
}
16✔
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, int m,
89✔
213
                         int n, int 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);
89✔
219
    Matrix *A = NULL, *AT = NULL;
89✔
220
    double *lhs_copy = NULL, *rhs_copy = NULL, *c_copy = NULL;
89✔
221
    int *row_sizes = NULL, *col_sizes = NULL;
89✔
222
    RowTag *row_tags = NULL;
89✔
223
    Lock *locks = NULL;
89✔
224
    Activity *activities = NULL;
89✔
225
    State *data = NULL;
89✔
226
    Constraints *constraints = NULL;
89✔
227
    Presolver *presolver = NULL;
89✔
228
    Objective *obj = NULL;
89✔
229
    ColTag *col_tags = NULL;
89✔
230
    Bound *bounds = NULL;
89✔
231
    Work *work = NULL;
89✔
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
    int n_rows = m;
89✔
239
    int n_cols = n;
89✔
240
    lhs_copy = (double *) ps_malloc(n_rows, sizeof(double));
89✔
241
    rhs_copy = (double *) ps_malloc(n_rows, sizeof(double));
89✔
242
    c_copy = (double *) ps_malloc(n_cols, sizeof(double));
89✔
243
    col_tags = (ColTag *) ps_calloc(n_cols, sizeof(ColTag));
89✔
244
    bounds = (Bound *) ps_malloc(n_cols, sizeof(Bound));
89✔
245
    work = new_work(n_rows, n_cols);
89✔
246
    row_sizes = (int *) ps_malloc(n_rows, sizeof(int));
89✔
247
    col_sizes = (int *) ps_malloc(n_cols, sizeof(int));
89✔
248

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

255
    memcpy(lhs_copy, lhs, n_rows * sizeof(double));
89✔
256
    memcpy(rhs_copy, rhs, n_rows * sizeof(double));
89✔
257
    memcpy(c_copy, c, n_cols * sizeof(double));
89✔
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);
89✔
264
    if (!A) goto cleanup;
89✔
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
    pthread_t thread_id;
273
    ParallelInitData parallel_data = {A,    work,     n_cols,   n_rows,   lbs,
89✔
274
                                      ubs,  lhs_copy, rhs_copy, bounds,   col_tags,
275
                                      NULL, NULL,     NULL,     row_sizes};
276

277
    pthread_create(&thread_id, NULL, init_thread_func, &parallel_data);
89✔
278

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

288
    // sync threads
289
    pthread_join(thread_id, NULL);
89✔
290

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

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

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

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

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

338
    clock_gettime(CLOCK_MONOTONIC, &timer.end);
89✔
339
    presolver->stats->ps_time_init = GET_ELAPSED_SECONDS(timer);
89✔
340

341
    return presolver;
89✔
342

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

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

370
    if (complexity == MEDIUM &&
76✔
371
        nnz_after_cycle >= NNZ_REMOVED_CYCLE * nnz_before_cycle)
34✔
372
    {
373
        return true;
17✔
374
    }
375

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

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

384
    return false;
59✔
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,
76✔
396
                                           int nnz_before_phase, int nnz_after_phase)
397
{
398
    if (curr_complexity == FAST)
76✔
399
    {
400
        if (nnz_after_phase < NNZ_REMOVED_FAST * nnz_before_phase)
42✔
401
        {
402
            return FAST;
8✔
403
        }
404

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

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

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

426
    remove_variables_with_close_bounds(prob);
179✔
427

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

438
    do
439
    {
440
        status = remove_ston_rows(prob);
189✔
441
        RETURN_IF_INFEASIBLE(status);
189✔
442
    } while (status == REDUCED);
189✔
443

444
    status = remove_empty_rows(prob->constraints);
179✔
445
    RETURN_IF_INFEASIBLE(status);
179✔
446

447
    status = remove_empty_cols(prob);
179✔
448
    RETURN_IF_UNBNDORINFEAS(status);
179✔
449

450
    assert(prob->constraints->state->ston_rows->len == 0);
179✔
451
    assert(prob->constraints->state->empty_rows->len == 0);
179✔
452
    assert(prob->constraints->state->empty_cols->len == 0);
179✔
453
    return UNCHANGED;
179✔
454
}
455

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

467
    if (stgs->ston_cols)
42✔
468
    {
469
        status |= remove_ston_cols(prob);
39✔
470

471
        // after removing singleton columns, there can be new empty columns
472
        // and singleton rows, but no empty rows
473
        assert(prob->constraints->state->empty_rows->len == 0);
39✔
474
        status |= run_trivial_explorers(prob, stgs);
39✔
475
    }
476

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

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

502
    if (stgs->primal_propagation)
34✔
503
    {
504
        // stats
505
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
12✔
506
        nnz_before = *nnz;
12✔
507

508
        // the actual reduction
509
        status |= check_activities(prob);
12✔
510
        status |= propagate_primal(prob, stgs->finite_bound_tightening);
12✔
511

512
        // after dom prop propagation there can be new empty and ston rows
513
        status |= run_trivial_explorers(prob, stgs);
12✔
514

515
        // stats
516
        stats->nnz_removed_primal_propagation += nnz_before - *nnz;
12✔
517
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
12✔
518
        stats->ps_time_primal_propagation += GET_ELAPSED_SECONDS(timer);
12✔
519
    }
520

521
    if (stgs->parallel_rows)
34✔
522
    {
523
        // stats
524
        nnz_before = *nnz;
34✔
525
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
34✔
526

527
        // the actual reduction (after removing parallel rows there can
528
        // be no new empty or ston rows so no need to run trivial presolvers)
529
        status |= remove_parallel_rows(prob->constraints);
34✔
530
        assert(prob->constraints->state->empty_rows->len == 0);
34✔
531
        assert(prob->constraints->state->ston_rows->len == 0);
34✔
532

533
        // stats
534
        stats->nnz_removed_parallel_rows += nnz_before - *nnz;
34✔
535
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
34✔
536
        stats->ps_time_parallel_rows += GET_ELAPSED_SECONDS(timer);
34✔
537
    }
538

539
    if (stgs->parallel_cols)
34✔
540
    {
541
        // stats
542
        nnz_before = *nnz;
10✔
543
        clock_gettime(CLOCK_MONOTONIC, &timer.start);
10✔
544

545
        // the actual reduction (after removing parallel columns there can
546
        // be no new empty rows or empty cols but might be new stonrows, probably
547
        // although that's rare but it does happen
548
        status |= remove_parallel_cols(prob);
10✔
549
        assert(prob->constraints->state->empty_rows->len == 0);
10✔
550
        assert(prob->constraints->state->empty_cols->len == 0);
10✔
551
        status |= run_trivial_explorers(prob, stgs);
10✔
552
        assert(prob->constraints->state->ston_rows->len == 0);
10✔
553

554
        // stats
555
        stats->nnz_removed_parallel_cols += nnz_before - *nnz;
10✔
556
        clock_gettime(CLOCK_MONOTONIC, &timer.end);
10✔
557
        stats->ps_time_parallel_cols += GET_ELAPSED_SECONDS(timer);
10✔
558
    }
559

560
    return status;
34✔
561
}
562

563
void populate_presolved_problem(Presolver *presolver)
17✔
564
{
565
    PresolvedProblem *reduced_prob = presolver->reduced_prob;
17✔
566
    Constraints *constraints = presolver->prob->constraints;
17✔
567
    Matrix *A = constraints->A;
17✔
568
    reduced_prob->m = A->m;
17✔
569
    reduced_prob->n = A->n;
17✔
570
    reduced_prob->nnz = A->nnz;
17✔
571
    reduced_prob->Ax = A->x;
17✔
572
    reduced_prob->Ai = A->i;
17✔
573
    reduced_prob->rhs = constraints->rhs;
17✔
574
    reduced_prob->lhs = constraints->lhs;
17✔
575
    reduced_prob->c = presolver->prob->obj->c;
17✔
576
    reduced_prob->obj_offset = presolver->prob->obj->offset;
17✔
577

578
    // create bounds arrays
579
    reduced_prob->lbs = (double *) malloc(A->n * sizeof(double));
17✔
580
    reduced_prob->ubs = (double *) malloc(A->n * sizeof(double));
17✔
581
    for (int i = 0; i < A->n; i++)
113✔
582
    {
583
        reduced_prob->lbs[i] = constraints->bounds[i].lb;
96✔
584
        reduced_prob->ubs[i] = constraints->bounds[i].ub;
96✔
585
    }
586

587
    // create row pointers
588
    reduced_prob->Ap = (int *) malloc((A->m + 1) * sizeof(int));
17✔
589
    for (int i = 0; i < A->m + 1; i++)
95✔
590
    {
591
        reduced_prob->Ap[i] = A->p[i].start;
78✔
592
    }
593
}
17✔
594

595
static inline void print_start_message(const PresolveStats *stats)
17✔
596
{
597
    printf("\n\t       PSLP v%s - LP presolver \n\t(c) Daniel "
17✔
598
           "Cederberg, Stanford University, 2025\n",
599
           PSLP_VERSION);
600
    printf("Original problem:  %d rows, %d columns, %d nnz\n",
17✔
601
           stats->n_rows_original, stats->n_cols_original, stats->nnz_original);
17✔
602
}
17✔
603

604
static inline void print_end_message(const Matrix *A, const PresolveStats *stats)
17✔
605
{
606
    printf("Presolved problem: %d rows, %d columns, %d nnz\n", A->m, A->n, A->nnz);
17✔
607
    printf("PSLP init & run time : %.3f seconds, %.3f \n", stats->ps_time_init,
17✔
608
           stats->presolve_total_time);
17✔
609
}
17✔
610

611
PresolveStatus run_presolver(Presolver *presolver)
17✔
612
{
613
    Timer inner_timer, outer_timer;
614
    int nnz_before_cycle, nnz_after_cycle;
615
    int nnz_before_phase, nnz_after_phase;
616
    int nnz_before_reduction;
617
    Problem *prob = presolver->prob;
17✔
618
    PresolveStats *stats = presolver->stats;
17✔
619
    Matrix *A = presolver->prob->constraints->A;
17✔
620
    const Settings *stgs = presolver->stgs;
17✔
621
    Complexity curr_complexity = FAST;
17✔
622
    bool terminate = false;
17✔
623
    PresolveStatus status = UNCHANGED;
17✔
624
    clock_gettime(CLOCK_MONOTONIC, &outer_timer.start);
17✔
625

626
    if (stgs->verbose)
17✔
627
    {
628
        print_start_message(stats);
17✔
629
    }
630

631
    DEBUG(run_debugger(prob->constraints, false));
17✔
632

633
    // ------------------------------------------------------------------------
634
    // Main loop for the presolver. The logic is organized into phases and
635
    // cycles. A phase is a sequence of presolvers belonging to a certain
636
    // complexity class. The cycle resets after the medium presolvers.
637
    // ------------------------------------------------------------------------
638
    nnz_before_cycle = A->nnz;
17✔
639
    while (!terminate)
93✔
640
    {
641
        // before each phase we run the trivial presolvers
642
        nnz_before_reduction = A->nnz;
76✔
643
        status = run_trivial_explorers(prob, stgs);
76✔
644
        // run_trivial_explorers returns something else than UNCHANGED only
645
        // if the problem is detected to be infeasible or unbounded
646
        RETURN_IF_NOT_UNCHANGED(status);
76✔
647
        stats->nnz_removed_trivial += nnz_before_reduction - A->nnz;
76✔
648
        DEBUG(run_debugger(prob->constraints, false));
76✔
649
        nnz_before_phase = A->nnz;
76✔
650

651
        // apply presolvers belonging to a certain complexity class
652
        if (curr_complexity == FAST)
76✔
653
        {
654
            nnz_before_reduction = A->nnz;
42✔
655
            RUN_AND_TIME(run_fast_explorers, inner_timer, stats->ps_time_fast,
42✔
656
                         status, prob, stgs);
657
            stats->nnz_removed_fast += nnz_before_reduction - A->nnz;
42✔
658
        }
659
        else if (curr_complexity == MEDIUM)
34✔
660
        {
661
            RUN_AND_TIME(run_medium_explorers, inner_timer, stats->ps_time_medium,
34✔
662
                         status, prob, stgs, stats);
663
            nnz_after_cycle = A->nnz;
34✔
664
        }
665

666
        nnz_after_phase = A->nnz;
76✔
667

668
        terminate =
669
            update_termination(nnz_after_cycle, nnz_before_cycle, curr_complexity,
76✔
670
                               status, stgs->max_time, outer_timer);
76✔
671

672
        if (curr_complexity == MEDIUM)
76✔
673
        {
674
            // a new cycle starts after the medium presolvers
675
            nnz_before_cycle = nnz_after_cycle;
34✔
676
        }
677

678
        curr_complexity =
679
            update_complexity(curr_complexity, nnz_before_phase, nnz_after_phase);
76✔
680
    }
681

682
    if (stgs->relax_bounds)
17✔
683
    {
684
        remove_redundant_bounds(prob->constraints);
17✔
685
    }
686

687
    problem_clean(prob, true);
17✔
688
    DEBUG(run_debugger(prob->constraints, true));
17✔
689

690
    clock_gettime(CLOCK_MONOTONIC, &outer_timer.end);
17✔
691
    stats->n_rows_reduced = A->m;
17✔
692
    stats->n_cols_reduced = A->n;
17✔
693
    stats->nnz_reduced = A->nnz;
17✔
694
    DEBUG(run_debugger_stats_consistency_check(stats));
17✔
695
    populate_presolved_problem(presolver);
17✔
696
    stats->presolve_total_time = GET_ELAPSED_SECONDS(outer_timer);
17✔
697

698
    if (stgs->verbose)
17✔
699
    {
700
        print_end_message(A, stats);
17✔
701
    }
702

703
    return status;
17✔
704
}
705

706
void postsolve(Presolver *presolver, const double *x, const double *y,
16✔
707
               const double *z, double obj)
708
{
709
    Timer timer;
710
    Solution *sol = presolver->sol;
16✔
711
    PresolveStats *stats = presolver->stats;
16✔
712
    State *data = presolver->prob->constraints->state;
16✔
713
    PostsolveInfo *postsolve_info = data->postsolve_info;
16✔
714
    clock_gettime(CLOCK_MONOTONIC, &timer.start);
16✔
715
    postsolver_update(postsolve_info, stats->n_cols_reduced, stats->n_rows_reduced,
16✔
716
                      data->work->mappings->cols, data->work->mappings->rows);
16✔
717
    postsolver_run(postsolve_info, sol, x, y, z);
16✔
718
    sol->obj = obj + presolver->prob->obj->offset;
16✔
719
    clock_gettime(CLOCK_MONOTONIC, &timer.end);
16✔
720
    stats->ps_time_post_solve = GET_ELAPSED_SECONDS(timer);
16✔
721

722
    if (presolver->stgs->verbose)
16✔
723
    {
724
        printf("PSLP postsolve time: %.4f seconds\n", stats->ps_time_post_solve);
16✔
725
    }
726
}
16✔
727

728
void free_presolver(Presolver *presolver)
89✔
729
{
730
    if (presolver == NULL)
89✔
731
    {
732
        return;
×
733
    }
734

735
    free_problem(presolver->prob);
89✔
736
    PS_FREE(presolver->stats);
89✔
737

738
    if (presolver->reduced_prob)
89✔
739
    {
740
        PS_FREE(presolver->reduced_prob->Ap);
89✔
741
        PS_FREE(presolver->reduced_prob->lbs);
89✔
742
        PS_FREE(presolver->reduced_prob->ubs);
89✔
743
        PS_FREE(presolver->reduced_prob);
89✔
744
    }
745

746
    if (presolver->sol)
89✔
747
    {
748
        PS_FREE(presolver->sol->x);
89✔
749
        PS_FREE(presolver->sol->y);
89✔
750
        PS_FREE(presolver->sol->z);
89✔
751
        PS_FREE(presolver->sol);
89✔
752
    }
753

754
    PS_FREE(presolver);
89✔
755
}
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