1 /* Solve linear programming problems by either vertex/point enumeration
2 or the primal simplex algorithm.
3 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
4 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
5
6 This file is part of the Parma Polyhedra Library (PPL).
7
8 The PPL is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The PPL is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software Foundation,
20 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
21
22 For the most up-to-date information see the Parma Polyhedra Library
23 site: http://bugseng.com/products/ppl/ . */
24
25 #include "ppl-config.h"
26 #include "ppl_c.h"
27 #include <gmp.h>
28 #include <stdio.h>
29 #include <assert.h>
30 #include <limits.h>
31 #include <time.h>
32 #include <stdarg.h>
33 #include <stdlib.h>
34 #include <errno.h>
35 #include <string.h>
36 #include <math.h>
37
38 #if defined(PPL_HAVE_GLPK_GLPK_H)
39 #include <glpk/glpk.h>
40 #elif defined(PPL_HAVE_GLPK_H)
41 #include <glpk.h>
42 #endif
43
44 #ifdef PPL_HAVE_GETOPT_H
45 # include <getopt.h>
46
47 /* Try to accommodate non-GNU implementations of `getopt()'. */
48 #if !defined(no_argument) && defined(NO_ARG)
49 #define no_argument NO_ARG
50 #endif
51
52 #if !defined(required_argument) && defined(REQUIRED_ARG)
53 #define required_argument REQUIRED_ARG
54 #endif
55
56 #if !defined(optional_argument) && defined(OPTIONAL_ARG)
57 #define optional_argument OPTIONAL_ARG
58 #endif
59
60 #endif /* defined(PPL_HAVE_GETOPT_H) */
61
62 #ifdef PPL_HAVE_UNISTD_H
63 /* Include this for `getopt()': especially important if we do not have
64 <getopt.h>. */
65 # include <unistd.h>
66 #endif
67
68 #ifdef PPL_HAVE_SIGNAL_H
69 # include <signal.h>
70 #endif
71
72 #ifdef PPL_HAVE_SYS_TIME_H
73 # include <sys/time.h>
74 #endif
75
76 #ifdef PPL_HAVE_SYS_RESOURCE_H
77 /* This should be included after <time.h> and <sys/time.h> so as to make
78 sure we have the definitions for, e.g., `ru_utime'. */
79 # include <sys/resource.h>
80 #endif
81
82 #if PPL_VERSION_MAJOR == 0 && PPL_VERSION_MINOR < 10
83 # error "PPL version 0.10 or following is required"
84 #endif
85
86 static const char* ppl_source_version = PPL_VERSION;
87
88 #ifdef __GNUC__
89 # define ATTRIBUTE_UNUSED __attribute__ ((__unused__))
90 #else
91 # define ATTRIBUTE_UNUSED
92 #endif
93
94 #if PPL_HAVE_DECL_GETRUSAGE
95 # define PPL_LPSOL_SUPPORTS_TIMINGS
96 #endif
97
98 #if defined(PPL_HAVE_SYS_RESOURCE_H) \
99 && PPL_CXX_SUPPORTS_LIMITING_MEMORY \
100 && (defined(SA_ONESHOT) || defined(SA_RESETHAND))
101 # define PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
102 #endif
103
104 #ifdef PPL_HAVE_GETOPT_H
105 static struct option long_options[] = {
106 {"check", optional_argument, 0, 'c'},
107 {"help", no_argument, 0, 'h'},
108 {"incremental", no_argument, 0, 'i'},
109 {"min", no_argument, 0, 'm'},
110 {"max", no_argument, 0, 'M'},
111 {"no-optimization", no_argument, 0, 'n'},
112 {"no-mip", no_argument, 0, 'r'},
113 {"max-cpu", required_argument, 0, 'C'},
114 {"max-memory", required_argument, 0, 'R'},
115 {"output", required_argument, 0, 'o'},
116 {"pricing", required_argument, 0, 'p'},
117 {"enumerate", no_argument, 0, 'e'},
118 {"simplex", no_argument, 0, 's'},
119 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
120 {"timings", no_argument, 0, 't'},
121 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
122 {"verbosity", required_argument, 0, 'v'},
123 {"version", no_argument, 0, 'V'},
124 {0, 0, 0, 0}
125 };
126 #endif
127
128 #define USAGE_STRING0 \
129 "Usage: %s [OPTION]... [FILE]\n" \
130 "Reads a file in MPS format and attempts solution using the optimization\n" \
131 "algorithms provided by the PPL.\n\n" \
132 "Options:\n" \
133 " -c, --check[=THRESHOLD] checks the obtained results using GLPK;\n" \
134 " optima are checked with a tolerance of\n" \
135 " THRESHOLD (default %.10g); input data\n" \
136 " are also perturbed the same way as GLPK does\n" \
137 " -i, --incremental solves the problem incrementally\n"
138 #define USAGE_STRING1 \
139 " -m, --min minimizes the objective function\n" \
140 " -M, --max maximizes the objective function (default)\n" \
141 " -n, --no-optimization checks for satisfiability only\n" \
142 " -r, --no-mip consider integer variables as real variables\n" \
143 " -CSECS, --max-cpu=SECS limits CPU usage to SECS seconds\n" \
144 " -RMB, --max-memory=MB limits memory usage to MB megabytes\n" \
145 " -h, --help prints this help text to stdout\n" \
146 " -oPATH, --output=PATH appends output to PATH\n"
147 #define USAGE_STRING2 \
148 " -e, --enumerate use the (expensive!) enumeration method\n" \
149 " -pM, --pricing=M use pricing method M for simplex (assumes -s);\n" \
150 " M is an int from 0 to 2, default 0:\n" \
151 " 0 --> steepest-edge using floating point\n" \
152 " 1 --> steepest-edge using exact arithmetic\n" \
153 " 2 --> textbook\n" \
154 " -s, --simplex use the simplex method\n"
155 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
156 #define USAGE_STRING3 \
157 " -t, --timings prints timings to stderr\n"
158 #else /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
159 #define USAGE_STRING3 \
160 ""
161 #endif /* !defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
162 #define USAGE_STRING4 \
163 " -v, --verbosity=LEVEL sets verbosity level (from 0 to 4, default 3):\n" \
164 " 0 --> quiet: no output except for errors and\n" \
165 " explicitly required notifications\n" \
166 " 1 --> solver state only\n" \
167 " 2 --> state + optimal value\n" \
168 " 3 --> state + optimal value + optimum location\n" \
169 " 4 --> lots of output\n" \
170 " -V, --version prints version information to stdout\n"
171 #ifndef PPL_HAVE_GETOPT_H
172 #define USAGE_STRING5 \
173 "\n" \
174 "NOTE: this version does not support long options.\n"
175 #else /* defined(PPL_HAVE_GETOPT_H) */
176 #define USAGE_STRING5 \
177 ""
178 #endif /* !defined(PPL_HAVE_GETOPT_H) */
179 #define USAGE_STRING6 \
180 "\n" \
181 "Report bugs to <ppl-devel@cs.unipr.it>.\n"
182
183
184 #define OPTION_LETTERS "bc::eimnMC:R:ho:p:rstVv:"
185
186 static const char* program_name = 0;
187 static unsigned long max_bytes_of_virtual_memory = 0;
188 static const char* output_argument = 0;
189 FILE* output_file = NULL;
190 static int check_results = 0;
191 static int use_simplex = 0;
192 static int pricing_method = 0;
193 static int verbosity = 3;
194 static int maximize = 1;
195 static int incremental = 0;
196 static int no_optimization = 0;
197 static int no_mip = 0;
198 static int check_results_failed = 0;
199 static double check_threshold = 0.0;
200 static const double default_check_threshold = 0.000000001;
201
202 #ifdef PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
203 static unsigned long max_seconds_of_cpu_time = 0;
204 #endif /* defined (PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME) */
205
206 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
207 static int print_timings = 0;
208 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
209
210 static void
my_exit(int status)211 my_exit(int status) {
212 (void) ppl_finalize();
213 exit(status);
214 }
215
216 void
fatal(const char * format,...)217 fatal(const char* format, ...) {
218 va_list ap;
219 fprintf(stderr, "%s: ", program_name);
220 va_start(ap, format);
221 vfprintf(stderr, format, ap);
222 va_end(ap);
223 fprintf(stderr, "\n");
224 my_exit(1);
225 }
226
227 #if 0
228 static void
229 warning(const char* format, ...) {
230 va_list ap;
231 va_start(ap, format);
232 fprintf(stderr, "%s: warning: ", program_name);
233 vfprintf(stderr, format, ap);
234 fprintf(stderr, "\n");
235 va_end(ap);
236 }
237 #endif
238
239 static void
error(const char * format,...)240 error(const char* format, ...) {
241 va_list ap;
242 fprintf(stderr, "%s: ", program_name);
243 va_start(ap, format);
244 vfprintf(stderr, format, ap);
245 va_end(ap);
246 fprintf(stderr, "\n");
247 if (output_argument) {
248 va_start(ap, format);
249 vfprintf(output_file, format, ap);
250 va_end(ap);
251 fprintf(output_file, "\n");
252 }
253 }
254
255 static const char*
get_ppl_version()256 get_ppl_version() {
257 const char* p;
258 (void) ppl_version(&p);
259 return p;
260 }
261
262 static const char*
get_ppl_banner()263 get_ppl_banner() {
264 const char* p;
265 (void) ppl_banner(&p);
266 return p;
267 }
268
269 static void
process_options(int argc,char * argv[])270 process_options(int argc, char* argv[]) {
271 #ifdef PPL_HAVE_GETOPT_H
272 int option_index;
273 #endif
274 int enumerate_required = 0;
275 int simplex_required = 0;
276 int incremental_required = 0;
277 int no_optimization_required = 0;
278 int no_mip_required = 0;
279 int c;
280 char* endptr;
281 long l;
282 double d;
283
284 while (1) {
285 #ifdef PPL_HAVE_GETOPT_H
286 option_index = 0;
287 c = getopt_long(argc, argv, OPTION_LETTERS, long_options, &option_index);
288 #else
289 c = getopt(argc, argv, OPTION_LETTERS);
290 #endif
291 if (c == EOF)
292 break;
293
294 switch (c) {
295 case 0:
296 break;
297
298 case 'c':
299 check_results = 1;
300 if (optarg) {
301 d = strtod(optarg, &endptr);
302 if (*endptr || errno == ERANGE || d < 0.0)
303 fatal("only a non-negative floating point number can `-c'");
304 else
305 check_threshold = d;
306 }
307 else
308 check_threshold = default_check_threshold;
309 break;
310
311 case 'm':
312 maximize = 0;
313 break;
314
315 case 'M':
316 maximize = 1;
317 break;
318
319 case '?':
320 case 'h':
321 fprintf(stdout, USAGE_STRING0, argv[0], default_check_threshold);
322 fputs(USAGE_STRING1, stdout);
323 fputs(USAGE_STRING2, stdout);
324 fputs(USAGE_STRING3, stdout);
325 fputs(USAGE_STRING4, stdout);
326 my_exit(0);
327 break;
328
329 #ifdef PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
330
331 case 'C':
332 l = strtol(optarg, &endptr, 10);
333 if (*endptr || l < 0)
334 fatal("a non-negative integer must follow `-C'");
335 else
336 max_seconds_of_cpu_time = l;
337 break;
338
339 #endif /* defined (PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME) */
340
341 case 'R':
342 l = strtol(optarg, &endptr, 10);
343 if (*endptr || l < 0)
344 fatal("a non-negative integer must follow `-R'");
345 else if (((unsigned long) l) > ULONG_MAX/(1024*1024))
346 max_bytes_of_virtual_memory = ULONG_MAX;
347 else
348 max_bytes_of_virtual_memory = l*1024*1024;
349 break;
350
351 case 'o':
352 output_argument = optarg;
353 break;
354
355 case 'p':
356 l = strtol(optarg, &endptr, 10);
357 if (*endptr || l < 0 || l > 2)
358 fatal("0 or 1 or 2 must follow `-p'");
359 else
360 pricing_method = l;
361 break;
362
363 case 'e':
364 enumerate_required = 1;
365 break;
366
367 case 's':
368 simplex_required = 1;
369 break;
370
371 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
372
373 case 't':
374 print_timings = 1;
375 break;
376
377 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
378
379 case 'v':
380 l = strtol(optarg, &endptr, 10);
381 if (*endptr || l < 0 || l > 4)
382 fatal("verbosity must be an integer between 0 and 4");
383 else
384 verbosity = l;
385 break;
386
387 case 'V':
388 fprintf(stdout, "%s\n", PPL_VERSION);
389 my_exit(0);
390 break;
391
392 case 'i':
393 incremental_required = 1;
394 break;
395
396 case 'n':
397 no_optimization_required = 1;
398 break;
399
400 case 'r':
401 no_mip_required = 1;
402 break;
403
404 default:
405 abort();
406 }
407 }
408
409 if (enumerate_required
410 && (simplex_required
411 || incremental_required))
412 fatal("-e option is incompatible with -i and -s");
413
414 if (enumerate_required)
415 use_simplex = 0;
416 else if (simplex_required)
417 use_simplex = 1;
418
419 if (incremental_required)
420 incremental = 1;
421
422 if (no_optimization_required)
423 no_optimization = 1;
424
425 if (no_mip_required)
426 no_mip = 1;
427
428 if (optind >= argc) {
429 if (verbosity >= 4)
430 fprintf(stderr,
431 "Parma Polyhedra Library version:\n%s\n\n"
432 "Parma Polyhedra Library banner:\n%s\n",
433 get_ppl_version(),
434 get_ppl_banner());
435 else
436 fatal("no input files");
437 }
438
439 if (argc - optind > 1)
440 /* We have multiple input files. */
441 fatal("only one input file is accepted");
442
443 if (output_argument) {
444 output_file = fopen(output_argument, "a");
445 if (output_file == NULL)
446 fatal("cannot open output file `%s'", output_argument);
447 }
448 else
449 output_file = stdout;
450 }
451
452 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
453
454 /* To save the time when start_clock is called. */
455 static struct timeval saved_ru_utime;
456
457 static void
start_clock()458 start_clock() {
459 struct rusage rsg;
460 if (getrusage(RUSAGE_SELF, &rsg) != 0)
461 fatal("getrusage failed: %s", strerror(errno));
462 else
463 saved_ru_utime = rsg.ru_utime;
464 }
465
466 static void
print_clock(FILE * f)467 print_clock(FILE* f) {
468 struct rusage rsg;
469 if (getrusage(RUSAGE_SELF, &rsg) != 0)
470 fatal("getrusage failed: %s", strerror(errno));
471 else {
472 time_t current_secs = rsg.ru_utime.tv_sec;
473 time_t current_usecs = rsg.ru_utime.tv_usec;
474 time_t saved_secs = saved_ru_utime.tv_sec;
475 time_t saved_usecs = saved_ru_utime.tv_usec;
476 int secs;
477 int csecs;
478 if (current_usecs < saved_usecs) {
479 csecs = (((1000000 + current_usecs) - saved_usecs) + 5000) / 10000;
480 secs = (current_secs - saved_secs) -1;
481 }
482 else {
483 csecs = ((current_usecs - saved_usecs) + 5000) / 10000;
484 secs = current_secs - saved_secs;
485 }
486 assert(csecs >= 0 && csecs < 100 && secs >= 0);
487 fprintf(f, "%d.%.2d", secs, csecs);
488 }
489 }
490
491 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
492
493 #ifdef PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
494
495 void
set_alarm_on_cpu_time(unsigned seconds,void (* handler)(int))496 set_alarm_on_cpu_time(unsigned seconds, void (*handler)(int)) {
497 sigset_t mask;
498 struct sigaction s;
499 struct rlimit t;
500
501 sigemptyset(&mask);
502
503 s.sa_handler = handler;
504 s.sa_mask = mask;
505 #if defined(SA_ONESHOT)
506 s.sa_flags = SA_ONESHOT;
507 #elif defined(SA_RESETHAND)
508 s.sa_flags = SA_RESETHAND;
509 #else
510 # error "Either SA_ONESHOT or SA_RESETHAND must be defined."
511 #endif
512
513 if (sigaction(SIGXCPU, &s, 0) != 0)
514 fatal("sigaction failed: %s", strerror(errno));
515
516 if (getrlimit(RLIMIT_CPU, &t) != 0)
517 fatal("getrlimit failed: %s", strerror(errno));
518
519 if (seconds < t.rlim_cur) {
520 t.rlim_cur = seconds;
521 if (setrlimit(RLIMIT_CPU, &t) != 0)
522 fatal("setrlimit failed: %s", strerror(errno));
523 }
524 }
525
526 #endif /* defined(PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME) */
527
528 #if PPL_CXX_SUPPORTS_LIMITING_MEMORY && PPL_HAVE_DECL_RLIMIT_AS
529
530 void
limit_virtual_memory(unsigned long bytes)531 limit_virtual_memory(unsigned long bytes) {
532 struct rlimit t;
533
534 if (getrlimit(RLIMIT_AS, &t) != 0)
535 fatal("getrlimit failed: %s", strerror(errno));
536
537 if (bytes < t.rlim_cur) {
538 t.rlim_cur = bytes;
539 if (setrlimit(RLIMIT_AS, &t) != 0)
540 fatal("setrlimit failed: %s", strerror(errno));
541 }
542 }
543
544 #else
545
546 void
limit_virtual_memory(unsigned long bytes ATTRIBUTE_UNUSED)547 limit_virtual_memory(unsigned long bytes ATTRIBUTE_UNUSED) {
548 }
549
550 #endif /* !PPL_HAVE_DECL_RLIMIT_AS */
551
552 #ifdef PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
553
554 static void
my_timeout(int dummy ATTRIBUTE_UNUSED)555 my_timeout(int dummy ATTRIBUTE_UNUSED) {
556 fprintf(stderr, "TIMEOUT\n");
557 if (output_argument)
558 fprintf(output_file, "TIMEOUT\n");
559 my_exit(0);
560 }
561
562 #endif /* defined(PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME) */
563
564 static mpz_t tmp_z;
565 static mpq_t tmp1_q;
566 static mpq_t tmp2_q;
567 static ppl_Coefficient_t ppl_coeff;
568 static glp_prob* glpk_lp;
569 static int glpk_lp_num_int;
570 static ppl_dimension_type* integer_variables;
571
572 static void
maybe_check_results(const int ppl_status,const double ppl_optimum_value)573 maybe_check_results(const int ppl_status, const double ppl_optimum_value) {
574 const char* ppl_status_string;
575 const char* glpk_status_string;
576 int glpk_status;
577 int treat_as_lp = 0;
578 glp_smcp glpk_smcp;
579
580 if (!check_results)
581 return;
582
583 if (no_mip || glpk_lp_num_int == 0)
584 treat_as_lp = 1;
585
586 glp_set_obj_dir(glpk_lp, (maximize ? GLP_MAX : GLP_MIN));
587
588 glp_init_smcp(&glpk_smcp);
589 /* Disable GLPK output. */
590 glpk_smcp.msg_lev = GLP_MSG_OFF;
591
592 if (treat_as_lp) {
593 /* Set the problem class to LP: MIP problems are thus treated as
594 LP ones. */
595 glp_exact(glpk_lp, &glpk_smcp);
596 glpk_status = glp_get_status(glpk_lp);
597 }
598 else {
599 /* MIP case. */
600 glp_simplex(glpk_lp, &glpk_smcp);
601 glpk_status = glp_get_status(glpk_lp);
602 if (glpk_status != GLP_NOFEAS && glpk_status != GLP_UNBND) {
603 glp_iocp glpk_iocp;
604 glp_init_iocp(&glpk_iocp);
605 /* Disable GLPK output. */
606 glpk_iocp.msg_lev = GLP_MSG_OFF;
607 glp_intopt(glpk_lp, &glpk_iocp);
608 glpk_status = glp_mip_status(glpk_lp);
609 }
610 }
611 /* If no_optimization is enabled, the second case is not possibile. */
612 if (!((ppl_status == PPL_MIP_PROBLEM_STATUS_UNFEASIBLE
613 && glpk_status == GLP_NOFEAS)
614 || (ppl_status == PPL_MIP_PROBLEM_STATUS_UNBOUNDED
615 && glpk_status == GLP_UNBND)
616 || (ppl_status == PPL_MIP_PROBLEM_STATUS_OPTIMIZED
617 && (glpk_status == GLP_OPT
618 /* If no_optimization is enabled, check if the problem is
619 unbounded for GLPK. */
620 || (no_optimization && (glpk_status == GLP_UNBND
621 || glpk_status == GLP_UNDEF)))))) {
622
623 if (ppl_status == PPL_MIP_PROBLEM_STATUS_UNFEASIBLE)
624 ppl_status_string = "unfeasible";
625 else if (ppl_status == PPL_MIP_PROBLEM_STATUS_UNBOUNDED)
626 ppl_status_string = "unbounded";
627 else if (ppl_status == PPL_MIP_PROBLEM_STATUS_OPTIMIZED)
628 ppl_status_string = "optimizable";
629 else
630 ppl_status_string = "<?>";
631
632 switch (glpk_status) {
633 case GLP_NOFEAS:
634 glpk_status_string = "unfeasible";
635 break;
636 case GLP_UNBND:
637 glpk_status_string = "unbounded";
638 break;
639 case GLP_OPT:
640 glpk_status_string = "optimizable";
641 break;
642 case GLP_UNDEF:
643 glpk_status_string = "undefined";
644 break;
645 default:
646 glpk_status_string = "<?>";
647 break;
648 }
649
650 error("check failed: for GLPK the problem is %s, not %s",
651 glpk_status_string, ppl_status_string);
652
653 check_results_failed = 1;
654 }
655 else if (!no_optimization
656 && ppl_status == PPL_MIP_PROBLEM_STATUS_OPTIMIZED) {
657
658 double glpk_optimum_value
659 = (treat_as_lp ? glp_get_obj_val(glpk_lp) : glp_mip_obj_val(glpk_lp));
660
661 if (fabs(ppl_optimum_value - glpk_optimum_value) > check_threshold) {
662 error("check failed: for GLPK the problem's optimum is %.20g,"
663 " not %.20g", glpk_optimum_value, ppl_optimum_value);
664 check_results_failed = 1;
665 }
666 }
667 return;
668 }
669
670
671 static const char*
variable_output_function(ppl_dimension_type var)672 variable_output_function(ppl_dimension_type var) {
673 const char* name = glp_get_col_name(glpk_lp, var+1);
674 if (name != NULL)
675 return name;
676 else
677 return 0;
678 }
679
680 static void
add_constraints(ppl_Linear_Expression_t ppl_le,int type,mpq_t rational_lb,mpq_t rational_ub,mpz_t den_lcm,ppl_Constraint_System_t ppl_cs)681 add_constraints(ppl_Linear_Expression_t ppl_le,
682 int type, mpq_t rational_lb, mpq_t rational_ub, mpz_t den_lcm,
683 ppl_Constraint_System_t ppl_cs) {
684 ppl_Constraint_t ppl_c;
685 ppl_Linear_Expression_t ppl_le2;
686 switch (type) {
687 case GLP_FR:
688 break;
689
690 case GLP_LO:
691 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_lb));
692 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_lb));
693 mpz_neg(tmp_z, tmp_z);
694 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
695 ppl_Linear_Expression_add_to_inhomogeneous(ppl_le, ppl_coeff);
696 ppl_new_Constraint(&ppl_c, ppl_le, PPL_CONSTRAINT_TYPE_GREATER_OR_EQUAL);
697 if (verbosity >= 4) {
698 ppl_io_fprint_Constraint(output_file, ppl_c);
699 fprintf(output_file, "\n");
700 }
701 ppl_Constraint_System_insert_Constraint(ppl_cs, ppl_c);
702 ppl_delete_Constraint(ppl_c);
703 break;
704
705 case GLP_UP:
706 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_ub));
707 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_ub));
708 mpz_neg(tmp_z, tmp_z);
709 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
710 ppl_Linear_Expression_add_to_inhomogeneous(ppl_le, ppl_coeff);
711 ppl_new_Constraint(&ppl_c, ppl_le,
712 PPL_CONSTRAINT_TYPE_LESS_OR_EQUAL);
713 if (verbosity >= 4) {
714 ppl_io_fprint_Constraint(output_file, ppl_c);
715 fprintf(output_file, "\n");
716 }
717 ppl_Constraint_System_insert_Constraint(ppl_cs, ppl_c);
718 ppl_delete_Constraint(ppl_c);
719 break;
720
721 case GLP_DB:
722 ppl_new_Linear_Expression_from_Linear_Expression(&ppl_le2, ppl_le);
723
724 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_lb));
725 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_lb));
726 mpz_neg(tmp_z, tmp_z);
727 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
728 ppl_Linear_Expression_add_to_inhomogeneous(ppl_le, ppl_coeff);
729 ppl_new_Constraint(&ppl_c, ppl_le, PPL_CONSTRAINT_TYPE_GREATER_OR_EQUAL);
730 if (verbosity >= 4) {
731 ppl_io_fprint_Constraint(output_file, ppl_c);
732 fprintf(output_file, "\n");
733 }
734 ppl_Constraint_System_insert_Constraint(ppl_cs, ppl_c);
735 ppl_delete_Constraint(ppl_c);
736
737 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_ub));
738 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_ub));
739 mpz_neg(tmp_z, tmp_z);
740 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
741 ppl_Linear_Expression_add_to_inhomogeneous(ppl_le2, ppl_coeff);
742 ppl_new_Constraint(&ppl_c, ppl_le2, PPL_CONSTRAINT_TYPE_LESS_OR_EQUAL);
743 ppl_delete_Linear_Expression(ppl_le2);
744 if (verbosity >= 4) {
745 ppl_io_fprint_Constraint(output_file, ppl_c);
746 fprintf(output_file, "\n");
747 }
748 ppl_Constraint_System_insert_Constraint(ppl_cs, ppl_c);
749 ppl_delete_Constraint(ppl_c);
750 break;
751
752 case GLP_FX:
753 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_lb));
754 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_lb));
755 mpz_neg(tmp_z, tmp_z);
756 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
757 ppl_Linear_Expression_add_to_inhomogeneous(ppl_le, ppl_coeff);
758 ppl_new_Constraint(&ppl_c, ppl_le,
759 PPL_CONSTRAINT_TYPE_EQUAL);
760 if (verbosity >= 4) {
761 ppl_io_fprint_Constraint(output_file, ppl_c);
762 fprintf(output_file, "\n");
763 }
764 ppl_Constraint_System_insert_Constraint(ppl_cs, ppl_c);
765 ppl_delete_Constraint(ppl_c);
766 break;
767
768 default:
769 fatal("internal error");
770 break;
771 }
772 }
773
774 static int
solve_with_generators(ppl_Constraint_System_t ppl_cs,ppl_const_Linear_Expression_t ppl_objective_le,ppl_Coefficient_t optimum_n,ppl_Coefficient_t optimum_d,ppl_Generator_t point)775 solve_with_generators(ppl_Constraint_System_t ppl_cs,
776 ppl_const_Linear_Expression_t ppl_objective_le,
777 ppl_Coefficient_t optimum_n,
778 ppl_Coefficient_t optimum_d,
779 ppl_Generator_t point) {
780 ppl_Polyhedron_t ppl_ph;
781 int optimum_found = 0;
782 int empty;
783 int unbounded;
784 int included;
785
786 /* Create the polyhedron (recycling the data structures of ppl_cs). */
787 ppl_new_C_Polyhedron_recycle_Constraint_System(&ppl_ph, ppl_cs);
788
789 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
790
791 if (print_timings) {
792 fprintf(stderr, "Time to create a PPL polyhedron: ");
793 print_clock(stderr);
794 fprintf(stderr, " s\n");
795 start_clock();
796 }
797
798 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
799
800 empty = ppl_Polyhedron_is_empty(ppl_ph);
801
802 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
803
804 if (print_timings) {
805 fprintf(stderr, "Time to check for emptiness: ");
806 print_clock(stderr);
807 fprintf(stderr, " s\n");
808 start_clock();
809 }
810
811 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
812
813 if (empty) {
814 if (verbosity >= 1)
815 fprintf(output_file, "Unfeasible problem.\n");
816 maybe_check_results(PPL_MIP_PROBLEM_STATUS_UNFEASIBLE, 0.0);
817 goto exit;
818 }
819
820 if (!empty && no_optimization) {
821 if (verbosity >= 1)
822 fprintf(output_file, "Feasible problem.\n");
823 /* Kludge: let's pass PPL_MIP_PROBLEM_STATUS_OPTIMIZED,
824 to let work `maybe_check_results'. */
825 maybe_check_results(PPL_MIP_PROBLEM_STATUS_OPTIMIZED, 0.0);
826 goto exit;
827 }
828
829 /* Check whether the problem is unbounded. */
830 unbounded = maximize
831 ? !ppl_Polyhedron_bounds_from_above(ppl_ph, ppl_objective_le)
832 : !ppl_Polyhedron_bounds_from_below(ppl_ph, ppl_objective_le);
833
834 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
835
836 if (print_timings) {
837 fprintf(stderr, "Time to check for unboundedness: ");
838 print_clock(stderr);
839 fprintf(stderr, " s\n");
840 start_clock();
841 }
842
843 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
844
845 if (unbounded) {
846 if (verbosity >= 1)
847 fprintf(output_file, "Unbounded problem.\n");
848 maybe_check_results(PPL_MIP_PROBLEM_STATUS_UNBOUNDED, 0.0);
849 goto exit;
850 }
851
852 optimum_found = maximize
853 ? ppl_Polyhedron_maximize_with_point(ppl_ph, ppl_objective_le,
854 optimum_n, optimum_d, &included,
855 point)
856 : ppl_Polyhedron_minimize_with_point(ppl_ph, ppl_objective_le,
857 optimum_n, optimum_d, &included,
858 point);
859
860 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
861
862 if (print_timings) {
863 fprintf(stderr, "Time to find the optimum: ");
864 print_clock(stderr);
865 fprintf(stderr, " s\n");
866 start_clock();
867 }
868
869 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
870
871 if (!optimum_found)
872 fatal("internal error");
873
874 if (!included)
875 fatal("internal error");
876
877 exit:
878 ppl_delete_Polyhedron(ppl_ph);
879 return optimum_found;
880 }
881
882 static int
solve_with_simplex(ppl_const_Constraint_System_t cs,ppl_const_Linear_Expression_t objective,ppl_Coefficient_t optimum_n,ppl_Coefficient_t optimum_d,ppl_Generator_t point)883 solve_with_simplex(ppl_const_Constraint_System_t cs,
884 ppl_const_Linear_Expression_t objective,
885 ppl_Coefficient_t optimum_n,
886 ppl_Coefficient_t optimum_d,
887 ppl_Generator_t point) {
888 ppl_MIP_Problem_t ppl_mip;
889 int optimum_found = 0;
890 int pricing = 0;
891 int status = 0;
892 int satisfiable = 0;
893 ppl_dimension_type space_dim;
894 ppl_const_Constraint_t c;
895 ppl_const_Generator_t g;
896 ppl_Constraint_System_const_iterator_t i;
897 ppl_Constraint_System_const_iterator_t iend;
898 int counter;
899 int mode = maximize
900 ? PPL_OPTIMIZATION_MODE_MAXIMIZATION
901 : PPL_OPTIMIZATION_MODE_MINIMIZATION;
902
903 ppl_Constraint_System_space_dimension(cs, &space_dim);
904 ppl_new_MIP_Problem_from_space_dimension(&ppl_mip, space_dim);
905 switch (pricing_method) {
906 case 0:
907 pricing = PPL_MIP_PROBLEM_CONTROL_PARAMETER_PRICING_STEEPEST_EDGE_FLOAT;
908 break;
909 case 1:
910 pricing = PPL_MIP_PROBLEM_CONTROL_PARAMETER_PRICING_STEEPEST_EDGE_EXACT;
911 break;
912 case 2:
913 pricing = PPL_MIP_PROBLEM_CONTROL_PARAMETER_PRICING_TEXTBOOK;
914 break;
915 default:
916 fatal("ppl_lpsol internal error");
917 }
918 ppl_MIP_Problem_set_control_parameter(ppl_mip, pricing);
919 ppl_MIP_Problem_set_objective_function(ppl_mip, objective);
920 ppl_MIP_Problem_set_optimization_mode(ppl_mip, mode);
921 if (!no_mip)
922 ppl_MIP_Problem_add_to_integer_space_dimensions(ppl_mip, integer_variables,
923 glpk_lp_num_int);
924 if (incremental) {
925 /* Add the constraints of `cs' one at a time. */
926 ppl_new_Constraint_System_const_iterator(&i);
927 ppl_new_Constraint_System_const_iterator(&iend);
928 ppl_Constraint_System_begin(cs, i);
929 ppl_Constraint_System_end(cs, iend);
930
931 counter = 0;
932 while (!ppl_Constraint_System_const_iterator_equal_test(i, iend)) {
933 ++counter;
934 if (verbosity >= 4)
935 fprintf(output_file, "\nSolving constraint %d\n", counter);
936 ppl_Constraint_System_const_iterator_dereference(i, &c);
937 ppl_MIP_Problem_add_constraint(ppl_mip, c);
938
939 if (no_optimization) {
940 satisfiable = ppl_MIP_Problem_is_satisfiable(ppl_mip);
941 if (!satisfiable)
942 break;
943 }
944 else
945 status = ppl_MIP_Problem_solve(ppl_mip);
946 ppl_Constraint_System_const_iterator_increment(i);
947 }
948 ppl_delete_Constraint_System_const_iterator(i);
949 ppl_delete_Constraint_System_const_iterator(iend);
950 }
951
952 else {
953 ppl_MIP_Problem_add_constraints(ppl_mip, cs);
954 if (no_optimization)
955 satisfiable = ppl_MIP_Problem_is_satisfiable(ppl_mip);
956 else
957 status = ppl_MIP_Problem_solve(ppl_mip);
958 }
959
960 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
961
962 if (print_timings) {
963 fprintf(stderr, "Time to solve the problem: ");
964 print_clock(stderr);
965 fprintf(stderr, " s\n");
966 start_clock();
967 }
968
969 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
970
971 if ((no_optimization && !satisfiable)
972 || (!no_optimization && status == PPL_MIP_PROBLEM_STATUS_UNFEASIBLE)) {
973 if (verbosity >= 1)
974 fprintf(output_file, "Unfeasible problem.\n");
975 maybe_check_results(status, 0.0);
976 goto exit;
977 }
978 else if (no_optimization && satisfiable) {
979 if (verbosity >= 1)
980 fprintf(output_file, "Feasible problem.\n");
981 /* Kludge: let's pass PPL_MIP_PROBLEM_STATUS_OPTIMIZED,
982 to let work `maybe_check_results'. */
983 maybe_check_results(PPL_MIP_PROBLEM_STATUS_OPTIMIZED, 0.0);
984 goto exit;
985 }
986 else if (status == PPL_MIP_PROBLEM_STATUS_UNBOUNDED) {
987 if (verbosity >= 1)
988 fprintf(output_file, "Unbounded problem.\n");
989 maybe_check_results(status, 0.0);
990 goto exit;
991 }
992 else if (status == PPL_MIP_PROBLEM_STATUS_OPTIMIZED) {
993 ppl_MIP_Problem_optimal_value(ppl_mip, optimum_n, optimum_d);
994 ppl_MIP_Problem_optimizing_point(ppl_mip, &g);
995 ppl_assign_Generator_from_Generator(point, g);
996 optimum_found = 1;
997 goto exit;
998 }
999 else
1000 fatal("internal error");
1001
1002 exit:
1003 ppl_delete_MIP_Problem(ppl_mip);
1004 return optimum_found;
1005 }
1006
1007 extern void set_d_eps(mpq_t x, double val);
1008
1009 static void
set_mpq_t_from_double(mpq_t q,double d)1010 set_mpq_t_from_double(mpq_t q, double d) {
1011 if (check_results)
1012 set_d_eps(q, d);
1013 else
1014 mpq_set_d(q, d);
1015 }
1016
1017 static void
solve(char * file_name)1018 solve(char* file_name) {
1019 ppl_Constraint_System_t ppl_cs;
1020 #ifndef NDEBUG
1021 ppl_Constraint_System_t ppl_cs_copy;
1022 #endif
1023 ppl_Generator_t optimum_location;
1024 ppl_Linear_Expression_t ppl_le;
1025 int dimension, row, num_rows, column, nz, i, j, type;
1026 int* coefficient_index;
1027 double lb, ub;
1028 double* coefficient_value;
1029 mpq_t rational_lb, rational_ub;
1030 mpq_t* rational_coefficient;
1031 mpq_t* objective;
1032 ppl_Linear_Expression_t ppl_objective_le;
1033 ppl_Coefficient_t optimum_n;
1034 ppl_Coefficient_t optimum_d;
1035 mpq_t optimum;
1036 mpz_t den_lcm;
1037 int optimum_found;
1038 glp_mpscp glpk_mpscp;
1039
1040 glpk_lp = glp_create_prob();
1041 glp_init_mpscp(&glpk_mpscp);
1042
1043 if (verbosity == 0) {
1044 /* FIXME: find a way to suppress output from glp_read_mps. */
1045 }
1046
1047 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
1048
1049 if (print_timings)
1050 start_clock();
1051
1052 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
1053
1054 if (glp_read_mps(glpk_lp, GLP_MPS_FILE, &glpk_mpscp, file_name) != 0)
1055 fatal("cannot read MPS file `%s'", file_name);
1056
1057 #ifdef PPL_LPSOL_SUPPORTS_TIMINGS
1058
1059 if (print_timings) {
1060 fprintf(stderr, "Time to read the input file: ");
1061 print_clock(stderr);
1062 fprintf(stderr, " s\n");
1063 start_clock();
1064 }
1065
1066 #endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
1067
1068 glpk_lp_num_int = glp_get_num_int(glpk_lp);
1069
1070 if (glpk_lp_num_int > 0 && !no_mip && !use_simplex)
1071 fatal("the enumeration solving method can not handle MIP problems");
1072
1073 dimension = glp_get_num_cols(glpk_lp);
1074
1075 /* Read variables constrained to be integer. */
1076 if (glpk_lp_num_int > 0 && !no_mip && use_simplex) {
1077 if (verbosity >= 4)
1078 fprintf(output_file, "Integer variables:\n");
1079 integer_variables = (ppl_dimension_type*)
1080 malloc((glpk_lp_num_int + 1)*sizeof(ppl_dimension_type));
1081 for (i = 0, j = 0; i < dimension; ++i) {
1082 int col_kind = glp_get_col_kind(glpk_lp, i+1);
1083 if (col_kind == GLP_IV || col_kind == GLP_BV) {
1084 integer_variables[j] = i;
1085 if (verbosity >= 4) {
1086 ppl_io_fprint_variable(output_file, i);
1087 fprintf(output_file, " ");
1088 }
1089 ++j;
1090 }
1091 }
1092 }
1093 coefficient_index = (int*) malloc((dimension+1)*sizeof(int));
1094 coefficient_value = (double*) malloc((dimension+1)*sizeof(double));
1095 rational_coefficient = (mpq_t*) malloc((dimension+1)*sizeof(mpq_t));
1096
1097
1098 ppl_new_Constraint_System(&ppl_cs);
1099
1100 mpq_init(rational_lb);
1101 mpq_init(rational_ub);
1102 for (i = 1; i <= dimension; ++i)
1103 mpq_init(rational_coefficient[i]);
1104
1105 mpz_init(den_lcm);
1106
1107 if (verbosity >= 4)
1108 fprintf(output_file, "\nConstraints:\n");
1109
1110 /* Set up the row (ordinary) constraints. */
1111 num_rows = glp_get_num_rows(glpk_lp);
1112 for (row = 1; row <= num_rows; ++row) {
1113 /* Initialize the least common multiple computation. */
1114 mpz_set_si(den_lcm, 1);
1115 /* Set `nz' to the number of non-zero coefficients. */
1116 nz = glp_get_mat_row(glpk_lp, row, coefficient_index, coefficient_value);
1117 for (i = 1; i <= nz; ++i) {
1118 set_mpq_t_from_double(rational_coefficient[i], coefficient_value[i]);
1119 /* Update den_lcm. */
1120 mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_coefficient[i]));
1121 }
1122
1123 lb = glp_get_row_lb(glpk_lp, row);
1124 ub = glp_get_row_ub(glpk_lp, row);
1125
1126 set_mpq_t_from_double(rational_lb, lb);
1127 set_mpq_t_from_double(rational_ub, ub);
1128
1129 mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_lb));
1130 mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_ub));
1131
1132 ppl_new_Linear_Expression_with_dimension(&ppl_le, dimension);
1133
1134 for (i = 1; i <= nz; ++i) {
1135 mpz_mul(tmp_z, den_lcm, mpq_numref(rational_coefficient[i]));
1136 mpz_divexact(tmp_z, tmp_z, mpq_denref(rational_coefficient[i]));
1137 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
1138 ppl_Linear_Expression_add_to_coefficient(ppl_le, coefficient_index[i]-1,
1139 ppl_coeff);
1140 }
1141
1142 type = glp_get_row_type(glpk_lp, row);
1143 add_constraints(ppl_le, type, rational_lb, rational_ub, den_lcm, ppl_cs);
1144
1145 ppl_delete_Linear_Expression(ppl_le);
1146 }
1147
1148 free(coefficient_value);
1149 for (i = 1; i <= dimension; ++i)
1150 mpq_clear(rational_coefficient[i]);
1151 free(rational_coefficient);
1152 free(coefficient_index);
1153
1154 #ifndef NDEBUG
1155 ppl_new_Constraint_System_from_Constraint_System(&ppl_cs_copy, ppl_cs);
1156 #endif
1157
1158 /*
1159 FIXME: here we could build the polyhedron and minimize it before
1160 adding the variable bounds.
1161 */
1162
1163 /* Set up the columns constraints, i.e., variable bounds. */
1164 for (column = 1; column <= dimension; ++column) {
1165
1166 lb = glp_get_col_lb(glpk_lp, column);
1167 ub = glp_get_col_ub(glpk_lp, column);
1168
1169 set_mpq_t_from_double(rational_lb, lb);
1170 set_mpq_t_from_double(rational_ub, ub);
1171
1172 /* Initialize the least common multiple computation. */
1173 mpz_set_si(den_lcm, 1);
1174 mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_lb));
1175 mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_ub));
1176
1177 ppl_new_Linear_Expression_with_dimension(&ppl_le, dimension);
1178 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, den_lcm);
1179 ppl_Linear_Expression_add_to_coefficient(ppl_le, column-1, ppl_coeff);
1180
1181 type = glp_get_col_type(glpk_lp, column);
1182 add_constraints(ppl_le, type, rational_lb, rational_ub, den_lcm, ppl_cs);
1183
1184 ppl_delete_Linear_Expression(ppl_le);
1185 }
1186
1187 mpq_clear(rational_ub);
1188 mpq_clear(rational_lb);
1189
1190 /* Deal with the objective function. */
1191 objective = (mpq_t*) malloc((dimension+1)*sizeof(mpq_t));
1192
1193 /* Initialize the least common multiple computation. */
1194 mpz_set_si(den_lcm, 1);
1195
1196 mpq_init(objective[0]);
1197 set_mpq_t_from_double(objective[0], glp_get_obj_coef(glpk_lp, 0));
1198 for (i = 1; i <= dimension; ++i) {
1199 mpq_init(objective[i]);
1200 set_mpq_t_from_double(objective[i], glp_get_obj_coef(glpk_lp, i));
1201 /* Update den_lcm. */
1202 mpz_lcm(den_lcm, den_lcm, mpq_denref(objective[i]));
1203 }
1204
1205 /* Set the ppl_objective_le to be the objective function. */
1206 ppl_new_Linear_Expression_with_dimension(&ppl_objective_le, dimension);
1207 /* Set value for objective function's inhomogeneous term. */
1208 mpz_mul(tmp_z, den_lcm, mpq_numref(objective[0]));
1209 mpz_divexact(tmp_z, tmp_z, mpq_denref(objective[0]));
1210 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
1211 ppl_Linear_Expression_add_to_inhomogeneous(ppl_objective_le, ppl_coeff);
1212 /* Set values for objective function's variable coefficients. */
1213 for (i = 1; i <= dimension; ++i) {
1214 mpz_mul(tmp_z, den_lcm, mpq_numref(objective[i]));
1215 mpz_divexact(tmp_z, tmp_z, mpq_denref(objective[i]));
1216 ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
1217 ppl_Linear_Expression_add_to_coefficient(ppl_objective_le, i-1, ppl_coeff);
1218 }
1219
1220 if (verbosity >= 4) {
1221 fprintf(output_file, "Objective function:\n");
1222 if (mpz_cmp_si(den_lcm, 1) != 0)
1223 fprintf(output_file, "(");
1224 ppl_io_fprint_Linear_Expression(output_file, ppl_objective_le);
1225 }
1226
1227 for (i = 0; i <= dimension; ++i)
1228 mpq_clear(objective[i]);
1229 free(objective);
1230
1231 if (verbosity >= 4) {
1232 if (mpz_cmp_si(den_lcm, 1) != 0) {
1233 fprintf(output_file, ")/");
1234 mpz_out_str(output_file, 10, den_lcm);
1235 }
1236 fprintf(output_file, "\n%s\n",
1237 (maximize ? "Maximizing." : "Minimizing."));
1238 }
1239
1240 ppl_new_Coefficient(&optimum_n);
1241 ppl_new_Coefficient(&optimum_d);
1242 ppl_new_Generator_zero_dim_point(&optimum_location);
1243
1244 optimum_found = use_simplex
1245 ? solve_with_simplex(ppl_cs,
1246 ppl_objective_le,
1247 optimum_n,
1248 optimum_d,
1249 optimum_location)
1250 : solve_with_generators(ppl_cs,
1251 ppl_objective_le,
1252 optimum_n,
1253 optimum_d,
1254 optimum_location);
1255
1256 ppl_delete_Linear_Expression(ppl_objective_le);
1257
1258 if (glpk_lp_num_int > 0)
1259 free(integer_variables);
1260
1261 if (optimum_found) {
1262 mpq_init(optimum);
1263 ppl_Coefficient_to_mpz_t(optimum_n, tmp_z);
1264 mpq_set_num(optimum, tmp_z);
1265 ppl_Coefficient_to_mpz_t(optimum_d, tmp_z);
1266 mpz_mul(tmp_z, tmp_z, den_lcm);
1267 mpq_set_den(optimum, tmp_z);
1268 if (verbosity == 1)
1269 fprintf(output_file, "Optimized problem.\n");
1270 if (verbosity >= 2)
1271 fprintf(output_file, "Optimum value: %.10g\n", mpq_get_d(optimum));
1272 if (verbosity >= 3) {
1273 fprintf(output_file, "Optimum location:\n");
1274 ppl_Generator_divisor(optimum_location, ppl_coeff);
1275 ppl_Coefficient_to_mpz_t(ppl_coeff, tmp_z);
1276 for (i = 0; i < dimension; ++i) {
1277 mpz_set(mpq_denref(tmp1_q), tmp_z);
1278 ppl_Generator_coefficient(optimum_location, i, ppl_coeff);
1279 ppl_Coefficient_to_mpz_t(ppl_coeff, mpq_numref(tmp1_q));
1280 ppl_io_fprint_variable(output_file, i);
1281 fprintf(output_file, " = %.10g\n", mpq_get_d(tmp1_q));
1282 }
1283 }
1284 #ifndef NDEBUG
1285 {
1286 ppl_Polyhedron_t ph;
1287 unsigned int relation;
1288 ppl_new_C_Polyhedron_recycle_Constraint_System(&ph, ppl_cs_copy);
1289 ppl_delete_Constraint_System(ppl_cs_copy);
1290 relation = ppl_Polyhedron_relation_with_Generator(ph, optimum_location);
1291 ppl_delete_Polyhedron(ph);
1292 assert(relation == PPL_POLY_GEN_RELATION_SUBSUMES);
1293 }
1294 #endif
1295 maybe_check_results(PPL_MIP_PROBLEM_STATUS_OPTIMIZED,
1296 mpq_get_d(optimum));
1297 mpq_clear(optimum);
1298 }
1299
1300 ppl_delete_Constraint_System(ppl_cs);
1301 ppl_delete_Coefficient(optimum_d);
1302 ppl_delete_Coefficient(optimum_n);
1303 ppl_delete_Generator(optimum_location);
1304
1305 glp_delete_prob(glpk_lp);
1306 }
1307
1308 static void
error_handler(enum ppl_enum_error_code code,const char * description)1309 error_handler(enum ppl_enum_error_code code,
1310 const char* description) {
1311 if (output_argument)
1312 fprintf(output_file, "PPL error code %d: %s\n", code, description);
1313 fatal("PPL error code %d: %s", code, description);
1314 }
1315
1316 #if defined(NDEBUG)
1317
1318 #if !(defined(PPL_GLPK_HAS_GLP_TERM_OUT) && defined(GLP_OFF))
1319
1320 #if defined(PPL_GLPK_HAS_GLP_TERM_HOOK) \
1321 || defined(PPL_GLPK_HAS__GLP_LIB_PRINT_HOOK)
1322
1323 static int
glpk_message_interceptor(void * info,const char * msg)1324 glpk_message_interceptor(void* info, const char* msg) {
1325 (void) info;
1326 (void) msg;
1327 return 1;
1328 }
1329
1330 #elif defined(PPL_GLPK_HAS_LIB_SET_PRINT_HOOK)
1331
1332 static int
glpk_message_interceptor(void * info,char * msg)1333 glpk_message_interceptor(void* info, char* msg) {
1334 (void) info;
1335 (void) msg;
1336 return 1;
1337 }
1338
1339 #endif /* !(defined(PPL_GLPK_HAS_GLP_TERM_HOOK)
1340 || defined(PPL_GLPK_HAS__GLP_LIB_PRINT_HOOK))
1341 && defined(PPL_GLPK_HAS_LIB_SET_PRINT_HOOK) */
1342
1343 #endif /* !(defined(PPL_GLPK_HAS_GLP_TERM_OUT) && defined(GLP_OFF)) */
1344
1345 #endif /* defined(NDEBUG) */
1346
1347 int
main(int argc,char * argv[])1348 main(int argc, char* argv[]) {
1349 #if defined(PPL_GLPK_HAS__GLP_LIB_PRINT_HOOK)
1350 extern void _glp_lib_print_hook(int (*func)(void *info, const char *buf),
1351 void *info);
1352 #endif
1353 program_name = argv[0];
1354 if (ppl_initialize() < 0)
1355 fatal("cannot initialize the Parma Polyhedra Library");
1356
1357 /* The PPL solver does not use floating point numbers, except
1358 perhaps for the steepest edge heuristics. In contrast, GLPK does
1359 use them, so it is best to restore the rounding mode as it was
1360 prior to the PPL initialization. */
1361 if (ppl_restore_pre_PPL_rounding() < 0)
1362 fatal("cannot restore the rounding mode");
1363
1364 if (ppl_set_error_handler(error_handler) < 0)
1365 fatal("cannot install the custom error handler");
1366
1367 if (strcmp(ppl_source_version, get_ppl_version()) != 0)
1368 fatal("was compiled with PPL version %s, but linked with version %s",
1369 ppl_source_version, get_ppl_version());
1370
1371 if (ppl_io_set_variable_output_function(variable_output_function) < 0)
1372 fatal("cannot install the custom variable output function");
1373
1374 #if defined(NDEBUG)
1375 #if defined(PPL_GLPK_HAS_GLP_TERM_OUT) && defined(GLP_OFF)
1376 glp_term_out(GLP_OFF);
1377 #elif defined(PPL_GLPK_HAS_GLP_TERM_HOOK)
1378 glp_term_hook(glpk_message_interceptor, 0);
1379 #elif defined(PPL_GLPK_HAS__GLP_LIB_PRINT_HOOK)
1380 _glp_lib_print_hook(glpk_message_interceptor, 0);
1381 #elif defined(PPL_GLPK_HAS_LIB_SET_PRINT_HOOK)
1382 lib_set_print_hook(0, glpk_message_interceptor);
1383 #endif
1384 #endif
1385
1386 /* Process command line options. */
1387 process_options(argc, argv);
1388
1389 /* Initialize globals. */
1390 mpz_init(tmp_z);
1391 mpq_init(tmp1_q);
1392 mpq_init(tmp2_q);
1393 ppl_new_Coefficient(&ppl_coeff);
1394
1395 #ifdef PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME
1396
1397 if (max_seconds_of_cpu_time > 0)
1398 set_alarm_on_cpu_time(max_seconds_of_cpu_time, my_timeout);
1399
1400 #endif /* defined (PPL_LPSOL_SUPPORTS_LIMIT_ON_CPU_TIME) */
1401
1402 if (max_bytes_of_virtual_memory > 0)
1403 limit_virtual_memory(max_bytes_of_virtual_memory);
1404
1405 while (optind < argc) {
1406 if (check_results)
1407 check_results_failed = 0;
1408
1409 solve(argv[optind++]);
1410
1411 if (check_results && check_results_failed)
1412 break;
1413 }
1414
1415 /* Finalize globals. */
1416 ppl_delete_Coefficient(ppl_coeff);
1417 mpq_clear(tmp2_q);
1418 mpq_clear(tmp1_q);
1419 mpz_clear(tmp_z);
1420
1421 /* Close output file, if any. */
1422 if (output_argument)
1423 fclose(output_file);
1424
1425 my_exit((check_results && check_results_failed) ? 1 : 0);
1426
1427 /* This is just to avoid a compiler warning. */
1428 return 0;
1429 }
1430