← Back to branch summary

~ubuntu-branches/ubuntu/quantal/ppl/quantal

« back to all changes in this revision

Viewing changes to demos/ppl_lpsol/ppl_lpsol.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2011-04-21 09:08:39 UTC
  • mfrom: (3.2.3 sid)
  • Revision ID: james.westby@ubuntu.com-20110421090839-hiiobptchgzoqiik
Tags: 0.11.2-3ubuntu1
https://launchpad.net/bugs/697305
* Merge with Debian; remaining changes:
  - Ignore results running the testsuite on amd64. One failing test.
    See LP: #697305.
  - Disable building the SWI-Prolog bindings as its in universe.

Show diffs side-by-side

added added

removed removed

Lines of Context:
demos/ppl_lpsol/ppl_lpsol.c
1
1
/* Solve linear programming problems by either vertex/point enumeration
2
2
   or the primal simplex algorithm.
3
3
   Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
 
4
   Copyright (C) 2010-2011 BUGSENG srl (http://bugseng.com)
4
5
 
5
6
This file is part of the Parma Polyhedra Library (PPL).
6
7
 
26
27
#include <gmp.h>
27
28
#include <stdio.h>
28
29
#include <assert.h>
 
30
#include <limits.h>
29
31
#include <time.h>
30
32
#include <stdarg.h>
31
33
#include <stdlib.h>
127
129
  "Reads a file in MPS format and attempts solution using the optimization\n" \
128
130
  "algorithms provided by the PPL.\n\n"                                 \
129
131
  "Options:\n"                                                          \
130
 
  "  -c, --check[=THRESHOLD] checks the obtained results;  optima are checked\n" \
131
 
  "                          with a tolerance of THRESHOLD (default %.10g)\n" \
132
 
  "  -i, --incremental       solves the problem incrementally\n"        \
 
132
  "  -c, --check[=THRESHOLD] checks the obtained results using GLPK;\n" \
 
133
  "                          optima are checked with a tolerance of\n"  \
 
134
  "                          THRESHOLD (default %.10g);  input data\n"  \
 
135
  "                          are also perturbed the same way as GLPK does\n" \
 
136
  "  -i, --incremental       solves the problem incrementally\n"
 
137
#define USAGE_STRING1                                                   \
133
138
  "  -m, --min               minimizes the objective function\n"        \
134
 
  "  -M, --max               maximizes the objective function (default)\n"
135
 
#define USAGE_STRING1                                                   \
 
139
  "  -M, --max               maximizes the objective function (default)\n" \
136
140
  "  -n, --no-optimization   checks for satisfiability only\n"          \
137
141
  "  -r, --no-mip            consider integer variables as real variables\n" \
138
142
  "  -CSECS, --max-cpu=SECS  limits CPU usage to SECS seconds\n"        \
139
143
  "  -RMB, --max-memory=MB   limits memory usage to MB megabytes\n"     \
140
144
  "  -h, --help              prints this help text to stdout\n"         \
141
 
  "  -oPATH, --output=PATH   appends output to PATH\n"                  \
142
 
  "  -e, --enumerate         use the (expensive!) enumeration method\n"
 
145
  "  -oPATH, --output=PATH   appends output to PATH\n"
143
146
#define USAGE_STRING2                                                   \
 
147
  "  -e, --enumerate         use the (expensive!) enumeration method\n" \
144
148
  "  -pM, --pricing=M        use pricing method M for simplex (assumes -s);\n" \
145
149
  "                          M is an int from 0 to 2, default 0:\n"     \
146
150
  "                          0 --> steepest-edge using floating point\n" \
208
212
  exit(status);
209
213
}
210
214
 
211
 
static void
 
215
void
212
216
fatal(const char* format, ...) {
213
217
  va_list ap;
214
218
  fprintf(stderr, "%s: ", program_name);
337
341
      l = strtol(optarg, &endptr, 10);
338
342
      if (*endptr || l < 0)
339
343
        fatal("a non-negative integer must follow `-R'");
 
344
      else if (((unsigned long) l) > ULONG_MAX/(1024*1024))
 
345
        max_bytes_of_virtual_memory = ULONG_MAX;
340
346
      else
341
347
        max_bytes_of_virtual_memory = l*1024*1024;
342
348
      break;
521
527
#if PPL_HAVE_DECL_RLIMIT_AS
522
528
 
523
529
void
524
 
limit_virtual_memory(unsigned bytes) {
 
530
limit_virtual_memory(unsigned long bytes) {
525
531
  struct rlimit t;
526
532
 
527
533
  if (getrlimit(RLIMIT_AS, &t) != 0)
537
543
#else
538
544
 
539
545
void
540
 
limit_virtual_memory(unsigned bytes ATTRIBUTE_UNUSED) {
 
546
limit_virtual_memory(unsigned long bytes ATTRIBUTE_UNUSED) {
541
547
}
542
548
 
543
549
#endif /* !PPL_HAVE_DECL_RLIMIT_AS */
584
590
    /* Set the problem class to LP: MIP problems are thus treated as
585
591
       LP ones. */
586
592
    lpx_set_class(glpk_lp, LPX_LP);
587
 
    lpx_simplex(glpk_lp);
 
593
    lpx_exact(glpk_lp);
588
594
    glpk_status = lpx_get_status(glpk_lp);
589
595
  }
590
596
  else {
649
655
      : lpx_mip_obj_val(glpk_lp);
650
656
 
651
657
    if (fabs(ppl_optimum_value - glpk_optimum_value) > check_threshold) {
652
 
      error("check failed: for GLPK the problem's optimum is %.10g,"
653
 
            " not %.10g", glpk_optimum_value, ppl_optimum_value);
 
658
      error("check failed: for GLPK the problem's optimum is %.20g,"
 
659
            " not %.20g", glpk_optimum_value, ppl_optimum_value);
654
660
      check_results_failed = 1;
655
661
    }
656
662
  }
768
774
                      ppl_Coefficient_t optimum_d,
769
775
                      ppl_Generator_t point) {
770
776
  ppl_Polyhedron_t ppl_ph;
 
777
  int optimum_found = 0;
771
778
  int empty;
772
779
  int unbounded;
773
780
  int included;
774
 
  int ok;
775
781
 
776
782
  /* Create the polyhedron (recycling the data structures of ppl_cs). */
777
783
  ppl_new_C_Polyhedron_recycle_Constraint_System(&ppl_ph, ppl_cs);
804
810
    if (verbosity >= 1)
805
811
      fprintf(output_file, "Unfeasible problem.\n");
806
812
    maybe_check_results(PPL_MIP_PROBLEM_STATUS_UNFEASIBLE, 0.0);
807
 
    return 0;
 
813
    goto exit;
808
814
  }
809
815
 
810
816
  if (!empty && no_optimization) {
813
819
    /* Kludge: let's pass PPL_MIP_PROBLEM_STATUS_OPTIMIZED,
814
820
       to let work `maybe_check_results'. */
815
821
    maybe_check_results(PPL_MIP_PROBLEM_STATUS_OPTIMIZED, 0.0);
816
 
    return 0;
 
822
    goto exit;
817
823
  }
818
824
 
819
825
  /* Check whether the problem is unbounded. */
836
842
    if (verbosity >= 1)
837
843
      fprintf(output_file, "Unbounded problem.\n");
838
844
    maybe_check_results(PPL_MIP_PROBLEM_STATUS_UNBOUNDED, 0.0);
839
 
    return 0;
 
845
    goto exit;
840
846
  }
841
847
 
842
 
  ok = maximize
 
848
  optimum_found = maximize
843
849
    ? ppl_Polyhedron_maximize_with_point(ppl_ph, ppl_objective_le,
844
 
                              optimum_n, optimum_d, &included,
845
 
                              point)
 
850
                                         optimum_n, optimum_d, &included,
 
851
                                         point)
846
852
    : ppl_Polyhedron_minimize_with_point(ppl_ph, ppl_objective_le,
847
 
                              optimum_n, optimum_d, &included,
848
 
                              point);
849
 
 
850
 
  if (!ok)
851
 
    fatal("internal error");
852
 
 
853
 
  ppl_delete_Polyhedron(ppl_ph);
 
853
                                         optimum_n, optimum_d, &included,
 
854
                                         point);
854
855
 
855
856
#ifdef PPL_LPSOL_SUPPORTS_TIMINGS
856
857
 
863
864
 
864
865
#endif /* defined(PPL_LPSOL_SUPPORTS_TIMINGS) */
865
866
 
 
867
  if (!optimum_found)
 
868
    fatal("internal error");
 
869
 
866
870
  if (!included)
867
871
    fatal("internal error");
868
872
 
869
 
  return 1;
 
873
 exit:
 
874
  ppl_delete_Polyhedron(ppl_ph);
 
875
  return optimum_found;
870
876
}
871
877
 
872
878
static int
876
882
                   ppl_Coefficient_t optimum_d,
877
883
                   ppl_Generator_t point) {
878
884
  ppl_MIP_Problem_t ppl_mip;
 
885
  int optimum_found = 0;
879
886
  int pricing = 0;
880
887
  int status = 0;
881
888
  int satisfiable = 0;
962
969
    if (verbosity >= 1)
963
970
      fprintf(output_file, "Unfeasible problem.\n");
964
971
    maybe_check_results(status, 0.0);
965
 
    return 0;
 
972
    goto exit;
966
973
  }
967
974
  else if (no_optimization && satisfiable) {
968
975
    if (verbosity >= 1)
970
977
    /* Kludge: let's pass PPL_MIP_PROBLEM_STATUS_OPTIMIZED,
971
978
       to let work `maybe_check_results'. */
972
979
    maybe_check_results(PPL_MIP_PROBLEM_STATUS_OPTIMIZED, 0.0);
973
 
    return 0;
 
980
    goto exit;
974
981
  }
975
982
  else if (status == PPL_MIP_PROBLEM_STATUS_UNBOUNDED) {
976
983
    if (verbosity >= 1)
977
984
      fprintf(output_file, "Unbounded problem.\n");
978
985
    maybe_check_results(status, 0.0);
979
 
    return 0;
 
986
    goto exit;
980
987
  }
981
988
  else if (status == PPL_MIP_PROBLEM_STATUS_OPTIMIZED) {
982
989
    ppl_MIP_Problem_optimal_value(ppl_mip, optimum_n, optimum_d);
983
990
    ppl_MIP_Problem_optimizing_point(ppl_mip, &g);
984
991
    ppl_assign_Generator_from_Generator(point, g);
985
 
    return 1;
 
992
    optimum_found = 1;
 
993
    goto exit;
986
994
  }
987
995
  else
988
996
    fatal("internal error");
989
997
 
990
 
  /* This is just to avoid a compiler warning. */
991
 
  return 0;
 
998
 exit:
 
999
  ppl_delete_MIP_Problem(ppl_mip);
 
1000
  return optimum_found;
 
1001
}
 
1002
 
 
1003
static void
 
1004
set_mpq_t_from_double(mpq_t q, double d) {
 
1005
  void set_d_eps(mpq_t x, double val);
 
1006
  if (check_results)
 
1007
    set_d_eps(q, d);
 
1008
  else
 
1009
    mpq_set_d(q, d);
992
1010
}
993
1011
 
994
1012
static void
1088
1106
    /* Set `nz' to the number of non-zero coefficients. */
1089
1107
    nz = lpx_get_mat_row(glpk_lp, row, coefficient_index, coefficient_value);
1090
1108
    for (i = 1; i <= nz; ++i) {
1091
 
      mpq_set_d(rational_coefficient[i], coefficient_value[i]);
 
1109
      set_mpq_t_from_double(rational_coefficient[i], coefficient_value[i]);
1092
1110
      /* Update den_lcm. */
1093
1111
      mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_coefficient[i]));
1094
1112
    }
1095
1113
    lpx_get_row_bnds(glpk_lp, row, &type, &lb, &ub);
1096
 
    mpq_set_d(rational_lb, lb);
 
1114
    set_mpq_t_from_double(rational_lb, lb);
1097
1115
    mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_lb));
1098
 
    mpq_set_d(rational_ub, ub);
 
1116
    set_mpq_t_from_double(rational_ub, ub);
1099
1117
    mpz_lcm(den_lcm, den_lcm, mpq_denref(rational_ub));
1100
1118
 
1101
1119
    ppl_new_Linear_Expression_with_dimension(&ppl_le, dimension);
1132
1150
  for (column = 1; column <= dimension; ++column) {
1133
1151
    lpx_get_col_bnds(glpk_lp, column, &type, &lb, &ub);
1134
1152
 
1135
 
    mpq_set_d(rational_lb, lb);
1136
 
    mpq_set_d(rational_ub, ub);
 
1153
    set_mpq_t_from_double(rational_lb, lb);
 
1154
    set_mpq_t_from_double(rational_ub, ub);
1137
1155
 
1138
1156
    /* Initialize the least common multiple computation. */
1139
1157
    mpz_set_si(den_lcm, 1);
1159
1177
  mpz_set_si(den_lcm, 1);
1160
1178
 
1161
1179
  mpq_init(objective[0]);
1162
 
  mpq_set_d(objective[0], lpx_get_obj_coef(glpk_lp, 0));
 
1180
  set_mpq_t_from_double(objective[0], lpx_get_obj_coef(glpk_lp, 0));
1163
1181
  for (i = 1; i <= dimension; ++i) {
1164
1182
    mpq_init(objective[i]);
1165
 
    mpq_set_d(objective[i], lpx_get_obj_coef(glpk_lp, i));
 
1183
    set_mpq_t_from_double(objective[i], lpx_get_obj_coef(glpk_lp, i));
1166
1184
    /* Update den_lcm. */
1167
1185
    mpz_lcm(den_lcm, den_lcm, mpq_denref(objective[i]));
1168
1186
  }
1169
1187
 
1170
1188
  /* Set the ppl_objective_le to be the objective function. */
1171
1189
  ppl_new_Linear_Expression_with_dimension(&ppl_objective_le, dimension);
1172
 
  /* The inhomogeneous term is completely useless for our purpose. */
 
1190
  /* Set value for objective function's inhomogeneous term. */
 
1191
  mpz_mul(tmp_z, den_lcm, mpq_numref(objective[0]));
 
1192
  mpz_divexact(tmp_z, tmp_z, mpq_denref(objective[0]));
 
1193
  ppl_assign_Coefficient_from_mpz_t(ppl_coeff, tmp_z);
 
1194
  ppl_Linear_Expression_add_to_inhomogeneous(ppl_objective_le, ppl_coeff);
 
1195
  /* Set values for objective function's variable coefficients. */
1173
1196
  for (i = 1; i <= dimension; ++i) {
1174
1197
    mpz_mul(tmp_z, den_lcm, mpq_numref(objective[i]));
1175
1198
    mpz_divexact(tmp_z, tmp_z, mpq_denref(objective[i]));
1179
1202
 
1180
1203
  if (verbosity >= 4) {
1181
1204
    fprintf(output_file, "Objective function:\n");
 
1205
    if (mpz_cmp_si(den_lcm, 1) != 0)
 
1206
      fprintf(output_file, "(");
1182
1207
    ppl_io_fprint_Linear_Expression(output_file, ppl_objective_le);
1183
1208
  }
1184
1209