1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /*                                                                           */
3 /*   File....: lpstore.c                                                     */
4 /*   Name....: Store Linear Programm                                         */
5 /*   Author..: Thorsten Koch                                                 */
6 /*   Copyright by Author, All rights reserved                                */
7 /*                                                                           */
8 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
9 /*
10  * Copyright (C) 2003-2018 by Thorsten Koch <koch@zib.de>
11  *
12  * This program is free software; you can redistribute it and/or
13  * modify it under the terms of the GNU Lesser General Public License
14  * as published by the Free Software Foundation; either version 3
15  * of the License, or (at your option) any later version.
16  *
17  * This program is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20  * GNU Lesser General Public License for more details.
21  *
22  * You should have received a copy of the GNU Lesser General Public License
23  * along with this program; if not, write to the Free Software
24  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
25  */
26 #if 0 /* Not used anymore ??? */
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <string.h>
30 #include <ctype.h>
31 #include <assert.h>
32 
33 #include <gmp.h>
34 
35 #include "zimpl/lint.h"
36 #include "zimpl/mshell.h"
37 #include <stdbool.h>
38 #include "zimpl/gmpmisc.h"
39 #include "zimpl/ratlptypes.h"
40 #include "zimpl/numb.h"
41 #include "zimpl/bound.h"
42 #include "zimpl/mme.h"
43 #include "zimpl/mono.h"
44 #include "zimpl/term.h"
45 #include "zimpl/ratlp.h"
46 #include "zimpl/ratlpstore.h"
47 
48 /*lint -e{818} supress "Pointer parameter 'var' could be declared as pointing to const" */
49 static void remove_fixed_var(
50    Lps* lp,
51    Var* var,
52    int  verbose_level)
53 {
54    Nzo*  nzo;
55    mpq_t x;
56    mpq_t temp;
57    bool  is_zero;
58 
59    assert(lp  != NULL);
60    assert(var != NULL);
61 
62    assert(var->type == VAR_FIXED && mpq_equal(var->lower, var->upper));
63 
64    if (verbose_level > 0)
65       printf("Removing variable %s fixed to %g\n", var->name, mpq_get_d(var->lower));
66 
67    mpq_init(x);
68    mpq_init(temp);
69 
70    is_zero = mpq_equal(var->lower, x /* zero */);
71 
72    while(var->first != NULL)
73    {
74       nzo = var->first;
75 
76       /* Do we have to ajust the lhs/rhs ?
77        */
78       if (!is_zero)
79       {
80          mpq_mul(x, var->lower, nzo->value);
81 
82          /* do we have a valid rhs ?
83           */
84          if (HAS_RHS(nzo->con))
85          {
86             mpq_sub(temp, nzo->con->rhs, x);
87             lps_setrhs(nzo->con, temp);
88          }
89          /* do we have a valid lhs ?
90           */
91          if (HAS_LHS(nzo->con))
92          {
93             mpq_sub(temp, nzo->con->lhs, x);
94             lps_setlhs(nzo->con, temp);
95          }
96       }
97       lps_delnzo(lp, nzo);
98    }
99    mpq_clear(temp);
100    mpq_clear(x);
101 }
102 
103 static PSResult simple_rows(
104    Lps*  lp,
105    Bool* again,
106    int   verbose_level)
107 {
108    bool  have_up;
109    bool  have_lo;
110    mpq_t up;
111    mpq_t lo;
112    Nzo*  nzo;
113    Var*  var;
114    Con*  con;
115    Con*  con_next;
116    int   rem_rows = 0;
117    int   rem_nzos = 0;
118 
119    mpq_init(up);
120    mpq_init(lo);
121 
122    for(con = lp->con_root; con != NULL; con = con_next)
123    {
124       con_next = con->next;
125 
126       /* infeasible range row
127        */
128       if (con->type == CON_RANGE && mpq_cmp(con->lhs, con->rhs) > 0)
129       {
130          printf("Infeasible range row %s lhs=%g rhs=%g\n",
131             con->name, mpq_get_d(con->lhs), mpq_get_d(con->rhs));
132 
133          return PRESOLVE_INFEASIBLE;
134       }
135       /* empty row ?
136        */
137       if (con->size == 0)
138       {
139          if (  (HAS_RHS(con) && mpq_cmp(con->rhs, const_zero) < 0)
140             || (HAS_LHS(con) && mpq_cmp(con->lhs, const_zero) > 0))
141          {
142             printf("Infeasible empty row %s lhs=%g rhs=%g\n",
143                con->name, mpq_get_d(con->lhs), mpq_get_d(con->rhs));
144 
145             return PRESOLVE_INFEASIBLE;
146          }
147          if (verbose_level > 0)
148             printf("Empty row %s removed\n", con->name);
149 
150          lps_delcon(lp, con);
151 
152          rem_rows++;
153          continue;
154       }
155       /* unconstraint constraint
156        */
157       if (con->type == CON_FREE)
158       {
159          if (verbose_level > 0)
160             printf("Unconstraint row %s removed\n", con->name);
161 
162          rem_rows++;
163          rem_nzos += con->size;
164 
165          lps_delcon(lp, con);
166          continue;
167       }
168       /* row singleton
169        */
170       if (con->size == 1)
171       {
172          have_up = false;
173          have_lo = false;
174          nzo     = con->first;
175          var     = nzo->var;
176 
177          if (mpq_cmp(nzo->value, const_zero) > 0) /* x > 0 */
178          {
179             if (HAS_RHS(con))
180             {
181                mpq_div(up, con->rhs, nzo->value);
182                have_up = true;
183             }
184             if (HAS_LHS(con))
185             {
186                mpq_div(lo, con->lhs, nzo->value);
187                have_lo = true;
188             }
189          }
190          else if (mpq_cmp(nzo->value, const_zero) < 0) /* x < 0 */
191          {
192             if (HAS_RHS(con))
193             {
194                mpq_div(lo, con->rhs, nzo->value);
195                have_lo = true;
196             }
197             if (HAS_LHS(con))
198             {
199                mpq_div(up, con->lhs, nzo->value);
200                have_up = true;
201             }
202          }
203          else if ((HAS_RHS(con) && !mpq_equal(con->rhs, const_zero))
204             ||    (HAS_LHS(con) && !mpq_equal(con->lhs, const_zero)))
205          {
206             /* x == 0 rhs/lhs != 0 */
207             printf("Infeasibel row %s Zero row singleton with non zero lhs or rhs\n",
208                con->name);
209 
210             return PRESOLVE_INFEASIBLE;
211          }
212          assert(!HAS_LOWER(var) || !HAS_UPPER(var) || mpq_cmp(var->lower, var->upper) <= 0);
213 
214          if (have_up && (!HAS_UPPER(var) || mpq_cmp(up, var->upper) < 0))
215             lps_setupper(var, up);
216 
217          if (have_lo && (!HAS_LOWER(var) || mpq_cmp(lo, var->lower) > 0))
218             lps_setlower(var, lo);
219 
220          if (HAS_LOWER(var) && HAS_UPPER(var) && mpq_cmp(var->lower, var->upper) > 0)
221          {
222             printf("Row %s implise infeasible bounds on var %s, lower=%g upper=%g\n",
223                con->name, var->name, mpq_get_d(var->lower), mpq_get_d(var->upper));
224 
225             return PRESOLVE_INFEASIBLE;
226          }
227 
228          if (verbose_level > 1)
229             printf("Row %s singleton var %s set to lower=%g upper=%g\n",
230                con->name, var->name, mpq_get_d(var->lower), mpq_get_d(var->upper));
231 
232          rem_rows++;
233          rem_nzos++;
234 
235          lps_delcon(lp, con);
236          continue;
237       }
238    }
239    assert(rem_rows != 0 || rem_nzos == 0);
240 
241    if (rem_rows > 0)
242    {
243       *again = true;
244 
245       if (verbose_level > 0)
246          printf("Simple row presolve removed %d rows and %d non-zeros\n",
247             rem_rows, rem_nzos);
248    }
249 
250    mpq_clear(up);
251    mpq_clear(lo);
252 
253    return PRESOLVE_OKAY;
254 }
255 
256 static PSResult handle_col_singleton(
257    Lps* lp,
258    Var* var,
259    int  verbose_level,
260    int* rem_cols,
261    int* rem_nzos)
262 {
263    Nzo*  nzo;
264    Con*  con;
265    int   cmpres;
266    mpq_t maxobj;
267 
268    assert(lp            != NULL);
269    assert(var           != NULL);
270    assert(verbose_level >= 0);
271    assert(rem_cols      != NULL);
272    assert(rem_nzos      != NULL);
273 
274    nzo = var->first;
275    con = nzo->con;
276 
277    assert(!mpq_equal(nzo->value, const_zero));
278 
279    mpq_init(maxobj);
280    mpq_set(maxobj, var->cost);
281 
282    if (lp->direct == LP_MIN)
283       mpq_neg(maxobj, maxobj);
284 
285    /* Value is positive
286     */
287    if (mpq_cmp(nzo->value, const_zero) > 0.0)
288    {
289       /* max -3 x
290        * s.t. 5 x (<= 8)
291        * l <= x
292        *
293        * and we have NO lhs, (rhs does not matter)
294        */
295       if (!HAS_LHS(con))
296       {
297          cmpres = mpq_cmp(maxobj, const_zero);
298 
299          /* The objective is either zero or negative
300           */
301          if (cmpres <= 0)
302          {
303             /* If we have no lower bound
304              */
305             if (!HAS_LOWER(var))
306             {
307                /* ... but a negative objective, so we get unbounded
308                 */
309                if (cmpres < 0)
310                {
311                   printf("Unbounded var %s\n", var->name);
312                   return PRESOLVE_UNBOUNDED;
313                }
314                /* With a zero objective and no bounds, there is not
315                 * much we can do.
316                 */
317             }
318             else
319             {
320                /* now we know we want to go to the lower bound.
321                 */
322                lps_setupper(var, var->lower);
323 
324                remove_fixed_var(lp, var, verbose_level);
325 
326                (*rem_cols)++;
327                (*rem_nzos)++;
328 
329                return PRESOLVE_OKAY;
330             }
331          }
332       }
333       /* max  3 x
334        * s.t. 5 x (>= 8)
335        * x <= u
336        *
337        * and we have NO rhs, (lhs does not matter)
338        */
339       if (!HAS_RHS(con))
340       {
341          cmpres = mpq_cmp(maxobj, const_zero);
342 
343          /* The objective is either zero or positive
344           */
345          if (cmpres >= 0)
346          {
347             /* If we have no upper bound
348              */
349             if (!HAS_UPPER(var))
350             {
351                /* ... but a positive objective, so we get unbounded
352                 */
353                if (cmpres > 0)
354                {
355                   printf("Unbounded var %s\n", var->name);
356                   return PRESOLVE_UNBOUNDED;
357                }
358                /* With a zero objective and no bounds, there is not
359                 * much we can do.
360                 */
361             }
362             else
363             {
364                /* now we know we want to go to the upper bound.
365                 */
366                lps_setlower(var, var->upper);
367 
368                remove_fixed_var(lp, var, verbose_level);
369 
370                (*rem_cols)++;
371                (*rem_nzos)++;
372 
373                return PRESOLVE_OKAY;
374             }
375          }
376       }
377    }
378    else
379    {
380       /* Value is negative
381        */
382       assert(mpq_cmp(nzo->value, const_zero) < 0.0);
383 
384       /* max  -3 x
385        * s.t. -5 x (>= 8)
386        * l <= x
387        */
388       if (!HAS_RHS(con))
389       {
390          cmpres = mpq_cmp(maxobj, const_zero);
391 
392          if (cmpres <= 0)
393          {
394             if (!HAS_LOWER(var))
395             {
396                if (cmpres < 0)
397                {
398                   printf("Unbounded var %s\n", var->name);
399                   return PRESOLVE_UNBOUNDED;
400                }
401             }
402             else
403             {
404                lps_setupper(var, var->lower);
405 
406                remove_fixed_var(lp, var, verbose_level);
407 
408                (*rem_cols)++;
409                (*rem_nzos)++;
410 
411                return PRESOLVE_OKAY;
412             }
413          }
414       }
415       /* max  3 x
416        * s.t. -5 x (<= 8)
417        * x <= u
418        */
419       if (!HAS_LHS(con))
420       {
421          cmpres = mpq_cmp(maxobj, const_zero);
422 
423          if (cmpres >= 0)
424          {
425             if (!HAS_UPPER(var))
426             {
427                if (cmpres > 0)
428                {
429                   printf("Unbounded var %s\n", var->name);
430                   return PRESOLVE_UNBOUNDED;
431                }
432             }
433             else
434             {
435                lps_setlower(var, var->upper);
436 
437                remove_fixed_var(lp, var, verbose_level);
438 
439                (*rem_cols)++;
440                (*rem_nzos)++;
441 
442                return PRESOLVE_OKAY;
443             }
444          }
445       }
446    }
447    mpq_clear(maxobj);
448 
449    return PRESOLVE_OKAY;
450 }
451 
452 static PSResult simple_cols(
453    Lps*  lp,
454    Bool* again,
455    int   verbose_level)
456 {
457    PSResult res;
458    mpq_t    maxobj;
459    int      rem_cols = 0;
460    int      rem_nzos = 0;
461    Var*     var;
462 
463    mpq_init(maxobj);
464 
465    for(var = lp->var_root; var != NULL; var = var->next)
466    {
467       if (var->type == VAR_FIXED && var->size == 0)
468          continue;
469 
470       /* Empty column ?
471        */
472       if (var->size == 0)
473       {
474          mpq_set(maxobj, var->cost);
475 
476          if (lp->direct == LP_MIN)
477             mpq_neg(maxobj, maxobj);
478 
479          if (mpq_cmp(maxobj, const_zero) > 0)
480          {
481             /* Do we not have an upper bound ?
482              */
483             if (!HAS_UPPER(var))
484             {
485                printf("Var %s unbounded\n", var->name);
486 
487                return PRESOLVE_UNBOUNDED;
488             }
489             lps_setlower(var, var->upper);
490          }
491          else if (mpq_cmp(maxobj, const_zero) < 0)
492          {
493             /* Do we not have an lower bound ?
494              */
495             if (!HAS_LOWER(var))
496             {
497                printf("Var %s unbounded\n", var->name);
498 
499                return PRESOLVE_UNBOUNDED;
500             }
501             lps_setupper(var, var->lower);
502          }
503          else
504          {
505             assert(mpq_equal(maxobj, const_zero));
506 
507             /* any value within the bounds is ok
508              */
509             if (HAS_LOWER(var))
510                lps_setupper(var, var->lower);
511             else if (HAS_UPPER(var))
512                lps_setlower(var, var->upper);
513             else
514             {
515                lps_setlower(var, const_zero);
516                lps_setupper(var, const_zero);
517             }
518          }
519          if (verbose_level > 1)
520             printf("Empty var %s fixed\n", var->name);
521 
522          rem_cols++;
523          continue;
524       }
525 
526       /* infeasible bounds ?
527        */
528       if (var->type == VAR_BOXED && mpq_cmp(var->lower, var->upper) > 0)
529       {
530          printf("Var %s infeasible bounds lower=%g upper=%g\n",
531             var->name, mpq_get_d(var->lower), mpq_get_d(var->upper));
532 
533          return PRESOLVE_INFEASIBLE;
534       }
535       /* Fixed column ?
536        */
537       if (var->type == VAR_FIXED)
538       {
539          assert(var->size > 0);
540 
541          rem_cols++;
542          rem_nzos += var->size;
543 
544          remove_fixed_var(lp, var, verbose_level);
545 
546          continue;
547       }
548       /* Column singleton
549        */
550       if (var->size == 1)
551       {
552          res = handle_col_singleton(lp, var, verbose_level, &rem_cols, &rem_nzos);
553 
554          if (res != PRESOLVE_OKAY)
555             return res;
556       }
557    }
558    assert(rem_cols != 0 || rem_nzos == 0);
559 
560    if (rem_cols > 0)
561    {
562       *again = true;
563 
564       if (verbose_level > 0)
565          printf("Simple col presolve removed %d cols and %d non-zeros\n",
566             rem_cols, rem_nzos);
567    }
568 
569    mpq_clear(maxobj);
570 
571    return PRESOLVE_OKAY;
572 }
573 
574 
575 
576 PSResult lps_presolve(Lps* lp, int verbose_level)
577 {
578    PSResult ret = PRESOLVE_OKAY;
579    bool     again;
580    /*
581    bool     rcagain;
582    bool     rragain;
583    */
584    do
585    {
586       again = false;
587 
588       if (ret == PRESOLVE_OKAY)
589          ret = simple_rows(lp, &again, verbose_level);
590 
591       if (ret == PRESOLVE_OKAY)
592          ret = simple_cols(lp, &again, verbose_level);
593 
594       assert(ret == PRESOLVE_OKAY || again == false);
595    }
596    while(again);
597 
598 #if 0
599    if (ret == OKAY)
600       ret = redundantCols(lp, rcagain);
601 
602    if (ret == OKAY)
603       ret = redundantRows(lp, rragain);
604 
605    again = (ret == OKAY) && (rcagain || rragain);
606 
607    /* This has to be a loop, otherwise we could end up with
608     * empty rows.
609     */
610    while(again)
611    {
612       again = false;
613 
614       if (ret == OKAY)
615          ret = simpleRows(lp, again);
616 
617       if (ret == OKAY)
618          ret = simpleCols(lp, again);
619 
620       assert(ret == OKAY || again == false);
621    }
622    VERBOSE1({ std::cout << "IREDSM25 redundant simplifier removed "
623                         << m_remRows << " rows, "
624                         << m_remNzos << " nzos, changed "
625                         << m_chgBnds << " col bounds "
626                         << m_chgLRhs << " row bounds,"
627                         << std::endl; });
628 
629 #endif
630    /* ??? does not work, because vars are not deleted */
631    if (lp->vars == 0 || lp->cons == 0)
632    {
633       /* ??? Check if this is ok.
634        */
635       assert(lp->vars == 0 && lp->cons == 0);
636 
637       printf("Simplifier removed all variables and constraints\n");
638       ret = PRESOLVE_VANISHED;
639    }
640    return ret;
641 }
642 
643 #endif
644 
645