diff --git a/README.md b/README.md index c7140e7..8d7fb1e 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ This repository contains a stripped-down version of the [mccs solver](http://www.i3s.unice.fr/~cpjm/misc/mccs.html), taken from snapshot 1.1, with a binding as an OCaml library, and building with `dune`. The [GLPK](https://www.gnu.org/software/glpk/glpk.html) source it links against is -also included within src/glpk, at version 4.63 (unmodified, apart from many +also included within src/glpk, at version 4.65 (unmodified, apart from many removed modules, corresponding to the parts that we don't use). The binding enables interoperation with binary CUDF data from diff --git a/src/glpk/api/cpxbas.c b/src/glpk/api/cpxbas.c index 7fd1528..e1c656a 100644 --- a/src/glpk/api/cpxbas.c +++ b/src/glpk/api/cpxbas.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2008-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2008-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "env.h" #include "prob.h" @@ -36,11 +32,7 @@ struct var /* penalty value */ }; -#ifndef _MSC_VER -static int fcmp(const void *ptr1, const void *ptr2) -#else -static int __cdecl fcmp(const void *ptr1, const void *ptr2) -#endif +static int CDECL fcmp(const void *ptr1, const void *ptr2) { /* this routine is passed to the qsort() function */ struct var *col1 = (void *)ptr1, *col2 = (void *)ptr2; if (col1->q < col2->q) return -1; diff --git a/src/glpk/api/prob1.c b/src/glpk/api/prob1.c index ff49524..6afad44 100644 --- a/src/glpk/api/prob1.c +++ b/src/glpk/api/prob1.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000-2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2000-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -22,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /* CAUTION: DO NOT CHANGE THE LIMITS BELOW */ diff --git a/src/glpk/api/wrmip.c b/src/glpk/api/wrmip.c new file mode 100644 index 0000000..407a5fe --- /dev/null +++ b/src/glpk/api/wrmip.c @@ -0,0 +1,122 @@ +/* wrmip.c (write MIP solution in GLPK format) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2010-2016 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "env.h" +#include "prob.h" + +/*********************************************************************** +* NAME +* +* glp_write_mip - write MIP solution in GLPK format +* +* SYNOPSIS +* +* int glp_write_mip(glp_prob *P, const char *fname); +* +* DESCRIPTION +* +* The routine glp_write_mip writes MIP solution to a text file in GLPK +* format. +* +* RETURNS +* +* If the operation was successful, the routine returns zero. Otherwise +* it prints an error message and returns non-zero. */ + +int glp_write_mip(glp_prob *P, const char *fname) +{ glp_file *fp; + GLPROW *row; + GLPCOL *col; + int i, j, count, ret = 1; + char *s; +#if 0 /* 04/IV-2016 */ + if (P == NULL || P->magic != GLP_PROB_MAGIC) + xerror("glp_write_mip: P = %p; invalid problem object\n", P); +#endif + if (fname == NULL) + xerror("glp_write_mip: fname = %d; invalid parameter\n", fname) + ; + xprintf("Writing MIP solution to '%s'...\n", fname); + fp = glp_open(fname, "w"), count = 0; + if (fp == NULL) + { xprintf("Unable to create '%s' - %s\n", fname, get_err_msg()); + goto done; + } + /* write comment lines */ + glp_format(fp, "c %-12s%s\n", "Problem:", + P->name == NULL ? "" : P->name), count++; + glp_format(fp, "c %-12s%d\n", "Rows:", P->m), count++; + glp_format(fp, "c %-12s%d\n", "Columns:", P->n), count++; + glp_format(fp, "c %-12s%d\n", "Non-zeros:", P->nnz), count++; + switch (P->mip_stat) + { case GLP_OPT: s = "INTEGER OPTIMAL"; break; + case GLP_FEAS: s = "INTEGER NON-OPTIMAL"; break; + case GLP_NOFEAS: s = "INTEGER EMPTY"; break; + case GLP_UNDEF: s = "INTEGER UNDEFINED"; break; + default: s = "???"; break; + } + glp_format(fp, "c %-12s%s\n", "Status:", s), count++; + switch (P->dir) + { case GLP_MIN: s = "MINimum"; break; + case GLP_MAX: s = "MAXimum"; break; + default: s = "???"; break; + } + glp_format(fp, "c %-12s%s%s%.10g (%s)\n", "Objective:", + P->obj == NULL ? "" : P->obj, + P->obj == NULL ? "" : " = ", P->mip_obj, s), count++; + glp_format(fp, "c\n"), count++; + /* write solution line */ + glp_format(fp, "s mip %d %d ", P->m, P->n), count++; + switch (P->mip_stat) + { case GLP_OPT: glp_format(fp, "o"); break; + case GLP_FEAS: glp_format(fp, "f"); break; + case GLP_NOFEAS: glp_format(fp, "n"); break; + case GLP_UNDEF: glp_format(fp, "u"); break; + default: glp_format(fp, "?"); break; + } + glp_format(fp, " %.*g\n", DBL_DIG, P->mip_obj); + /* write row solution descriptor lines */ + for (i = 1; i <= P->m; i++) + { row = P->row[i]; + glp_format(fp, "i %d %.*g\n", i, DBL_DIG, row->mipx), count++; + } + /* write column solution descriptor lines */ + for (j = 1; j <= P->n; j++) + { col = P->col[j]; + glp_format(fp, "j %d %.*g\n", j, DBL_DIG, col->mipx), count++; + } + /* write end line */ + glp_format(fp, "e o f\n"), count++; + if (glp_ioerr(fp)) + { xprintf("Write error on '%s' - %s\n", fname, get_err_msg()); + goto done; + } + /* MIP solution has been successfully written */ + xprintf("%d lines were written\n", count); + ret = 0; +done: if (fp != NULL) + glp_close(fp); + return ret; +} + +/* eof */ diff --git a/src/glpk/context_flags.ml b/src/glpk/context_flags.ml index d119501..db89072 100644 --- a/src/glpk/context_flags.ml +++ b/src/glpk/context_flags.ml @@ -4,7 +4,7 @@ match Sys.argv.(1) with | "clibs" -> if Sys.win32 && Config.ccomp_type = "msvc" then - print_string "(glpk_4_63.lib)" + print_string "(glpk_4_65.lib)" else let glpk = if Array.length Sys.argv > 2 && Sys.argv.(2) = "static" then diff --git a/src/glpk/bfd.c b/src/glpk/draft/bfd.c similarity index 100% rename from src/glpk/bfd.c rename to src/glpk/draft/bfd.c diff --git a/src/glpk/bfd.h b/src/glpk/draft/bfd.h similarity index 100% rename from src/glpk/bfd.h rename to src/glpk/draft/bfd.h diff --git a/src/glpk/draft.h b/src/glpk/draft/draft.h similarity index 100% rename from src/glpk/draft.h rename to src/glpk/draft/draft.h diff --git a/src/glpk/glpapi06.c b/src/glpk/draft/glpapi06.c similarity index 99% rename from src/glpk/glpapi06.c rename to src/glpk/draft/glpapi06.c index acadb81..a31e396 100644 --- a/src/glpk/glpapi06.c +++ b/src/glpk/draft/glpapi06.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,8 +23,8 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" -#include "glpnpp.h" +#include "ios.h" +#include "npp.h" #if 0 /* 07/XI-2015 */ #include "glpspx.h" #else diff --git a/src/glpk/glpapi09.c b/src/glpk/draft/glpapi09.c similarity index 99% rename from src/glpk/glpapi09.c rename to src/glpk/draft/glpapi09.c index 7b8fd18..0d3ab57 100644 --- a/src/glpk/glpapi09.c +++ b/src/glpk/draft/glpapi09.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -24,8 +24,8 @@ #include "draft.h" #include "env.h" -#include "glpios.h" -#include "glpnpp.h" +#include "ios.h" +#include "npp.h" /*********************************************************************** * NAME @@ -686,8 +686,12 @@ void glp_init_iocp(glp_iocp *parm) parm->save_sol = NULL; parm->alien = GLP_OFF; #endif +#if 0 /* 20/I-2018 */ #if 1 /* 16/III-2016; not documented--should not be used */ parm->flip = GLP_OFF; +#endif +#else + parm->flip = GLP_ON; #endif return; } diff --git a/src/glpk/glpapi10.c b/src/glpk/draft/glpapi10.c similarity index 100% rename from src/glpk/glpapi10.c rename to src/glpk/draft/glpapi10.c diff --git a/src/glpk/glpapi12.c b/src/glpk/draft/glpapi12.c similarity index 100% rename from src/glpk/glpapi12.c rename to src/glpk/draft/glpapi12.c diff --git a/src/glpk/glpapi13.c b/src/glpk/draft/glpapi13.c similarity index 99% rename from src/glpk/glpapi13.c rename to src/glpk/draft/glpapi13.c index a4ad10b..1181b39 100644 --- a/src/glpk/glpapi13.c +++ b/src/glpk/draft/glpapi13.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * NAME @@ -457,7 +457,11 @@ int glp_ios_pool_size(glp_tree *tree) if (tree->reason != GLP_ICUTGEN) xerror("glp_ios_pool_size: operation not allowed\n"); xassert(tree->local != NULL); +#ifdef NEW_LOCAL /* 02/II-2018 */ + return tree->local->m; +#else return tree->local->size; +#endif } /**********************************************************************/ diff --git a/src/glpk/glpios01.c b/src/glpk/draft/glpios01.c similarity index 95% rename from src/glpk/glpios01.c rename to src/glpk/draft/glpios01.c index 13d04f8..cb1a0da 100644 --- a/src/glpk/glpios01.c +++ b/src/glpk/draft/glpios01.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" #include "misc.h" static int lpx_eval_tab_row(glp_prob *lp, int k, int ind[], @@ -125,13 +125,16 @@ glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm) tree->pred_type = NULL; tree->pred_lb = tree->pred_ub = NULL; tree->pred_stat = NULL; - /* cut generator */ + /* cut generators */ tree->local = ios_create_pool(tree); /*tree->first_attempt = 1;*/ /*tree->max_added_cuts = 0;*/ /*tree->min_eff = 0.0;*/ /*tree->miss = 0;*/ /*tree->just_selected = 0;*/ +#ifdef NEW_COVER /* 13/II-2018 */ + tree->cov_gen = NULL; +#endif tree->mir_gen = NULL; tree->clq_gen = NULL; /*tree->round = 0;*/ @@ -1384,6 +1387,15 @@ int ios_solve_node(glp_tree *tree) /**********************************************************************/ +#ifdef NEW_LOCAL /* 02/II-2018 */ +IOSPOOL *ios_create_pool(glp_tree *tree) +{ /* create cut pool */ + IOSPOOL *pool; + pool = glp_create_prob(); + glp_add_cols(pool, tree->mip->n); + return pool; +} +#else IOSPOOL *ios_create_pool(glp_tree *tree) { /* create cut pool */ IOSPOOL *pool; @@ -1398,7 +1410,23 @@ IOSPOOL *ios_create_pool(glp_tree *tree) pool->ord = 0, pool->curr = NULL; return pool; } +#endif +#ifdef NEW_LOCAL /* 02/II-2018 */ +int ios_add_row(glp_tree *tree, IOSPOOL *pool, + const char *name, int klass, int flags, int len, const int ind[], + const double val[], int type, double rhs) +{ /* add row (constraint) to the cut pool */ + int i; + i = glp_add_rows(pool, 1); + glp_set_row_name(pool, i, name); + pool->row[i]->klass = klass; + xassert(flags == 0); + glp_set_mat_row(pool, i, len, ind, val); + glp_set_row_bnds(pool, i, type, rhs, rhs); + return i; +} +#else int ios_add_row(glp_tree *tree, IOSPOOL *pool, const char *name, int klass, int flags, int len, const int ind[], const double val[], int type, double rhs) @@ -1457,7 +1485,14 @@ int ios_add_row(glp_tree *tree, IOSPOOL *pool, pool->size++; return pool->size; } +#endif +#ifdef NEW_LOCAL /* 02/II-2018 */ +IOSCUT *ios_find_row(IOSPOOL *pool, int i) +{ /* find row (constraint) in the cut pool */ + xassert(0); +} +#else IOSCUT *ios_find_row(IOSPOOL *pool, int i) { /* find row (constraint) in the cut pool */ /* (smart linear search) */ @@ -1509,7 +1544,14 @@ IOSCUT *ios_find_row(IOSPOOL *pool, int i) xassert(pool->curr != NULL); return pool->curr; } +#endif +#ifdef NEW_LOCAL /* 02/II-2018 */ +void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i) +{ /* remove row (constraint) from the cut pool */ + xassert(0); +} +#else void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i) { /* remove row (constraint) from the cut pool */ IOSCUT *cut; @@ -1553,7 +1595,22 @@ void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i) pool->size--; return; } +#endif +#ifdef NEW_LOCAL /* 02/II-2018 */ +void ios_clear_pool(glp_tree *tree, IOSPOOL *pool) +{ /* remove all rows (constraints) from the cut pool */ + if (pool->m > 0) + { int i, *num; + num = talloc(1+pool->m, int); + for (i = 1; i <= pool->m; i++) + num[i] = i; + glp_del_rows(pool, pool->m, num); + tfree(num); + } + return; +} +#else void ios_clear_pool(glp_tree *tree, IOSPOOL *pool) { /* remove all rows (constraints) from the cut pool */ xassert(pool != NULL); @@ -1574,7 +1631,16 @@ void ios_clear_pool(glp_tree *tree, IOSPOOL *pool) pool->ord = 0, pool->curr = NULL; return; } +#endif +#ifdef NEW_LOCAL /* 02/II-2018 */ +void ios_delete_pool(glp_tree *tree, IOSPOOL *pool) +{ /* delete cut pool */ + xassert(pool != NULL); + glp_delete_prob(pool); + return; +} +#else void ios_delete_pool(glp_tree *tree, IOSPOOL *pool) { /* delete cut pool */ xassert(pool != NULL); @@ -1582,9 +1648,10 @@ void ios_delete_pool(glp_tree *tree, IOSPOOL *pool) xfree(pool); return; } +#endif #if 1 /* 11/VII-2013 */ -#include "glpnpp.h" +#include "npp.h" void ios_process_sol(glp_tree *T) { /* process integer feasible solution just found */ @@ -1596,6 +1663,21 @@ void ios_process_sol(glp_tree *T) } xassert(T->P != NULL); /* save solution to text file, if requested */ + if (T->save_sol != NULL) + { char *fn, *mark; + fn = talloc(strlen(T->save_sol) + 50, char); + mark = strrchr(T->save_sol, '*'); + if (mark == NULL) + strcpy(fn, T->save_sol); + else + { memcpy(fn, T->save_sol, mark - T->save_sol); + fn[mark - T->save_sol] = '\0'; + sprintf(fn + strlen(fn), "%03d", ++(T->save_cnt)); + strcat(fn, &mark[1]); + } + glp_write_mip(T->P, fn); + tfree(fn); + } return; } #endif diff --git a/src/glpk/glpios02.c b/src/glpk/draft/glpios02.c similarity index 99% rename from src/glpk/glpios02.c rename to src/glpk/draft/glpios02.c index 43cff57..a73458a 100644 --- a/src/glpk/glpios02.c +++ b/src/glpk/draft/glpios02.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * prepare_row_info - prepare row info to determine implied bounds diff --git a/src/glpk/glpios03.c b/src/glpk/draft/glpios03.c similarity index 96% rename from src/glpk/glpios03.c rename to src/glpk/draft/glpios03.c index eeb67dd..21d6a00 100644 --- a/src/glpk/glpios03.c +++ b/src/glpk/draft/glpios03.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * show_progress - display current progress of the search @@ -72,7 +72,10 @@ static void show_progress(glp_tree *T, int bingo) else if (temp == +DBL_MAX) sprintf(best_bound, "%17s", "+inf"); else + { if (fabs(temp) < 1e-9) + temp = 0; sprintf(best_bound, "%17.9e", temp); + } } /* choose the relation sign between global bounds */ if (T->mip->dir == GLP_MIN) @@ -633,60 +636,6 @@ done: tfree(x); return ret; } -#if 0 -#define round_heur round_heur2 -static int round_heur(glp_tree *T) -{ glp_prob *lp; - int *ind, ret, i, j, len; - double *val; - lp = glp_create_prob(); - ind = talloc(1+T->mip->n, int); - val = talloc(1+T->mip->n, double); - glp_add_rows(lp, T->orig_m); - glp_add_cols(lp, T->n); - for (i = 1; i <= T->orig_m; i++) - { glp_set_row_bnds(lp, i, - T->orig_type[i], T->orig_lb[i], T->orig_ub[i]); - len = glp_get_mat_row(T->mip, i, ind, val); - glp_set_mat_row(lp, i, len, ind, val); - } - for (j = 1; j <= T->n; j++) - { GLPCOL *col = T->mip->col[j]; - glp_set_obj_coef(lp, j, col->coef); - if (col->kind == GLP_IV) - { /* integer variable */ - glp_set_col_bnds(lp, j, GLP_FX, floor(col->prim + .5), 0); - } - else - { glp_set_col_bnds(lp, j, T->orig_type[T->orig_m+j], - T->orig_lb[T->orig_m+j], T->orig_ub[T->orig_m+j]); - } - } -glp_term_out(GLP_OFF); - glp_adv_basis(lp, 0); - ret = glp_simplex(lp, NULL); -glp_term_out(GLP_ON); - if (ret != 0) - { ret = 1; - goto done; - } - if (glp_get_status(lp) != GLP_OPT) - { ret = 2; - goto done; - } - for (j = 1; j <= lp->n; j++) - val[j] = lp->col[j]->prim; - if (glp_ios_heur_sol(T, val) == 0) - ret = 0; - else - ret = 3; -done: glp_delete_prob(lp); - tfree(ind); - tfree(val); - return ret; -} -#endif - /**********************************************************************/ #if 1 /* 08/III-2016 */ @@ -715,6 +664,34 @@ static void gmi_gen(glp_tree *T) } #endif +#ifdef NEW_COVER /* 13/II-2018 */ +static void cov_gen(glp_tree *T) +{ /* generate cover cuts */ + glp_prob *P, *pool; + if (T->cov_gen == NULL) + return; + P = T->mip; + pool = glp_create_prob(); + glp_add_cols(pool, P->n); + glp_cov_gen1(P, T->cov_gen, pool); + if (pool->m > 0) + { int i, len, *ind; + double *val; + ind = xcalloc(1+P->n, sizeof(int)); + val = xcalloc(1+P->n, sizeof(double)); + for (i = 1; i <= pool->m; i++) + { len = glp_get_mat_row(pool, i, ind, val); + glp_ios_add_row(T, NULL, GLP_RF_COV, 0, len, ind, val, + GLP_UP, pool->row[i]->ub); + } + xfree(ind); + xfree(val); + } + glp_delete_prob(pool); + return; +} +#endif + #if 1 /* 08/III-2016 */ static void mir_gen(glp_tree *T) { /* generate mixed integer rounding cuts */ @@ -798,8 +775,11 @@ static void generate_cuts(glp_tree *T) } if (T->parm->cov_cuts == GLP_ON) { /* cover cuts works well along with mir cuts */ - /*if (T->round <= 5)*/ - ios_cov_gen(T); +#ifdef NEW_COVER /* 13/II-2018 */ + cov_gen(T); +#else + ios_cov_gen(T); +#endif } if (T->parm->clq_cuts == GLP_ON) { if (T->clq_gen != NULL) @@ -942,7 +922,11 @@ int ios_driver(glp_tree *T) #endif #if 1 /* 16/III-2016 */ if (((glp_iocp *)T->parm)->flip) +#if 0 /* 20/I-2018 */ xprintf("WARNING: LONG-STEP DUAL SIMPLEX WILL BE USED\n"); +#else + xprintf("Long-step dual simplex will be used\n"); +#endif #endif /* on entry to the B&B driver it is assumed that the active list contains the only active (i.e. root) subproblem, which is the @@ -1038,6 +1022,10 @@ int ios_driver(glp_tree *T) if (T->parm->cov_cuts == GLP_ON) { if (T->parm->msg_lev >= GLP_MSG_ALL) xprintf("Cover cuts enabled\n"); +#ifdef NEW_COVER /* 13/II-2018 */ + xassert(T->cov_gen == NULL); + T->cov_gen = glp_cov_init(T->mip); +#endif } if (T->parm->clq_cuts == GLP_ON) { xassert(T->clq_gen == NULL); @@ -1356,7 +1344,11 @@ int ios_driver(glp_tree *T) #endif /* it's time to generate cutting planes */ xassert(T->local != NULL); +#ifdef NEW_LOCAL /* 02/II-2018 */ + xassert(T->local->m == 0); +#else xassert(T->local->size == 0); +#endif /* let the application program generate some cuts; note that it can add cuts either to the local cut pool or directly to the current subproblem */ @@ -1402,7 +1394,11 @@ int ios_driver(glp_tree *T) } /* if the local cut pool is not empty, select useful cuts and add them to the current subproblem */ +#ifdef NEW_LOCAL /* 02/II-2018 */ + if (T->local->m > 0) +#else if (T->local->size > 0) +#endif { xassert(T->reason == 0); T->reason = GLP_ICUTGEN; ios_process_cuts(T); @@ -1498,6 +1494,10 @@ int ios_driver(glp_tree *T) ios_mir_term(T->mir_gen), T->mir_gen = NULL; #else glp_mir_free(T->mir_gen), T->mir_gen = NULL; +#endif +#ifdef NEW_COVER /* 13/II-2018 */ + if (T->cov_gen != NULL) + glp_cov_free(T->cov_gen), T->cov_gen = NULL; #endif if (T->clq_gen != NULL) #if 0 /* 08/III-2016 */ diff --git a/src/glpk/glpios07.c b/src/glpk/draft/glpios07.c similarity index 99% rename from src/glpk/glpios07.c rename to src/glpk/draft/glpios07.c index 3ea515b..f750e57 100644 --- a/src/glpk/glpios07.c +++ b/src/glpk/draft/glpios07.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*---------------------------------------------------------------------- -- COVER INEQUALITIES diff --git a/src/glpk/glpios09.c b/src/glpk/draft/glpios09.c similarity index 99% rename from src/glpk/glpios09.c rename to src/glpk/draft/glpios09.c index d87578c..d80ed9a 100644 --- a/src/glpk/glpios09.c +++ b/src/glpk/draft/glpios09.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * NAME diff --git a/src/glpk/glpios11.c b/src/glpk/draft/glpios11.c similarity index 63% rename from src/glpk/glpios11.c rename to src/glpk/draft/glpios11.c index 78b3187..09fccef 100644 --- a/src/glpk/glpios11.c +++ b/src/glpk/draft/glpios11.c @@ -4,9 +4,9 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013, 2017 Andrew Makhorin, Department for Applied -* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights -* reserved. E-mail: . +* 2009, 2010, 2011, 2013, 2017, 2018 Andrew Makhorin, Department for +* Applied Informatics, Moscow Aviation Institute, Moscow, Russia. All +* rights reserved. E-mail: . * * GLPK is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by @@ -22,13 +22,9 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "draft.h" #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * NAME @@ -70,11 +66,7 @@ struct info /* lower bound to objective degradation */ }; -#ifndef _MSC_VER -static int fcmp(const void *arg1, const void *arg2) -#else -static int __cdecl fcmp(const void *arg1, const void *arg2) -#endif +static int CDECL fcmp(const void *arg1, const void *arg2) { const struct info *info1 = arg1, *info2 = arg2; if (info1->deg == 0.0 && info2->deg == 0.0) { if (info1->eff > info2->eff) return -1; @@ -89,6 +81,138 @@ static int __cdecl fcmp(const void *arg1, const void *arg2) static double parallel(IOSCUT *a, IOSCUT *b, double work[]); +#ifdef NEW_LOCAL /* 02/II-2018 */ +void ios_process_cuts(glp_tree *T) +{ IOSPOOL *pool; + IOSCUT *cut; + GLPAIJ *aij; + struct info *info; + int k, kk, max_cuts, len, ret, *ind; + double *val, *work, rhs; + /* the current subproblem must exist */ + xassert(T->curr != NULL); + /* the pool must exist and be non-empty */ + pool = T->local; + xassert(pool != NULL); + xassert(pool->m > 0); + /* allocate working arrays */ + info = xcalloc(1+pool->m, sizeof(struct info)); + ind = xcalloc(1+T->n, sizeof(int)); + val = xcalloc(1+T->n, sizeof(double)); + work = xcalloc(1+T->n, sizeof(double)); + for (k = 1; k <= T->n; k++) work[k] = 0.0; + /* build the list of cuts stored in the cut pool */ + for (k = 1; k <= pool->m; k++) + info[k].cut = pool->row[k], info[k].flag = 0; + /* estimate efficiency of all cuts in the cut pool */ + for (k = 1; k <= pool->m; k++) + { double temp, dy, dz; + cut = info[k].cut; + /* build the vector of cut coefficients and compute its + Euclidean norm */ + len = 0; temp = 0.0; + for (aij = cut->ptr; aij != NULL; aij = aij->r_next) + { xassert(1 <= aij->col->j && aij->col->j <= T->n); + len++, ind[len] = aij->col->j, val[len] = aij->val; + temp += aij->val * aij->val; + } + if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON; + /* transform the cut to express it only through non-basic + (auxiliary and structural) variables */ + len = glp_transform_row(T->mip, len, ind, val); + /* determine change in the cut value and in the objective + value for the adjacent basis by simulating one step of the + dual simplex */ + switch (cut->type) + { case GLP_LO: rhs = cut->lb; break; + case GLP_UP: rhs = cut->ub; break; + default: xassert(cut != cut); + } + ret = _glp_analyze_row(T->mip, len, ind, val, cut->type, + rhs, 1e-9, NULL, NULL, NULL, NULL, &dy, &dz); + /* determine normalized residual and lower bound to objective + degradation */ + if (ret == 0) + { info[k].eff = fabs(dy) / sqrt(temp); + /* if some reduced costs violates (slightly) their zero + bounds (i.e. have wrong signs) due to round-off errors, + dz also may have wrong sign being close to zero */ + if (T->mip->dir == GLP_MIN) + { if (dz < 0.0) dz = 0.0; + info[k].deg = + dz; + } + else /* GLP_MAX */ + { if (dz > 0.0) dz = 0.0; + info[k].deg = - dz; + } + } + else if (ret == 1) + { /* the constraint is not violated at the current point */ + info[k].eff = info[k].deg = 0.0; + } + else if (ret == 2) + { /* no dual feasible adjacent basis exists */ + info[k].eff = 1.0; + info[k].deg = DBL_MAX; + } + else + xassert(ret != ret); + /* if the degradation is too small, just ignore it */ + if (info[k].deg < 0.01) info[k].deg = 0.0; + } + /* sort the list of cuts by decreasing objective degradation and + then by decreasing efficacy */ + qsort(&info[1], pool->m, sizeof(struct info), fcmp); + /* only first (most efficient) max_cuts in the list are qualified + as candidates to be added to the current subproblem */ + max_cuts = (T->curr->level == 0 ? 90 : 10); + if (max_cuts > pool->m) max_cuts = pool->m; + /* add cuts to the current subproblem */ +#if 0 + xprintf("*** adding cuts ***\n"); +#endif + for (k = 1; k <= max_cuts; k++) + { int i, len; + /* if this cut seems to be inefficient, skip it */ + if (info[k].deg < 0.01 && info[k].eff < 0.01) continue; + /* if the angle between this cut and every other cut included + in the current subproblem is small, skip this cut */ + for (kk = 1; kk < k; kk++) + { if (info[kk].flag) + { if (parallel(info[k].cut, info[kk].cut, work) > 0.90) + break; + } + } + if (kk < k) continue; + /* add this cut to the current subproblem */ +#if 0 + xprintf("eff = %g; deg = %g\n", info[k].eff, info[k].deg); +#endif + cut = info[k].cut, info[k].flag = 1; + i = glp_add_rows(T->mip, 1); + if (cut->name != NULL) + glp_set_row_name(T->mip, i, cut->name); + xassert(T->mip->row[i]->origin == GLP_RF_CUT); + T->mip->row[i]->klass = cut->klass; + len = 0; + for (aij = cut->ptr; aij != NULL; aij = aij->r_next) + len++, ind[len] = aij->col->j, val[len] = aij->val; + glp_set_mat_row(T->mip, i, len, ind, val); + switch (cut->type) + { case GLP_LO: rhs = cut->lb; break; + case GLP_UP: rhs = cut->ub; break; + default: xassert(cut != cut); + } + glp_set_row_bnds(T->mip, i, cut->type, rhs, rhs); + } + /* free working arrays */ + xfree(info); + xfree(ind); + xfree(val); + xfree(work); + return; +} +#else void ios_process_cuts(glp_tree *T) { IOSPOOL *pool; IOSCUT *cut; @@ -211,6 +335,7 @@ void ios_process_cuts(glp_tree *T) xfree(work); return; } +#endif #if 0 /*********************************************************************** @@ -269,6 +394,25 @@ static double efficacy(glp_tree *T, IOSCUT *cut) * i.e. with disjoint support, while requirement cos phi <= 0.999 means * only avoiding duplicate (parallel) cuts [1]. */ +#ifdef NEW_LOCAL /* 02/II-2018 */ +static double parallel(IOSCUT *a, IOSCUT *b, double work[]) +{ GLPAIJ *aij; + double s = 0.0, sa = 0.0, sb = 0.0, temp; + for (aij = a->ptr; aij != NULL; aij = aij->r_next) + { work[aij->col->j] = aij->val; + sa += aij->val * aij->val; + } + for (aij = b->ptr; aij != NULL; aij = aij->r_next) + { s += work[aij->col->j] * aij->val; + sb += aij->val * aij->val; + } + for (aij = a->ptr; aij != NULL; aij = aij->r_next) + work[aij->col->j] = 0.0; + temp = sqrt(sa) * sqrt(sb); + if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON; + return s / temp; +} +#else static double parallel(IOSCUT *a, IOSCUT *b, double work[]) { IOSAIJ *aij; double s = 0.0, sa = 0.0, sb = 0.0, temp; @@ -286,5 +430,6 @@ static double parallel(IOSCUT *a, IOSCUT *b, double work[]) if (temp < DBL_EPSILON * DBL_EPSILON) temp = DBL_EPSILON; return s / temp; } +#endif /* eof */ diff --git a/src/glpk/glpios12.c b/src/glpk/draft/glpios12.c similarity index 98% rename from src/glpk/glpios12.c rename to src/glpk/draft/glpios12.c index d5cf302..bec6fa2 100644 --- a/src/glpk/glpios12.c +++ b/src/glpk/draft/glpios12.c @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +23,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" /*********************************************************************** * NAME diff --git a/src/glpk/glpscl.c b/src/glpk/draft/glpscl.c similarity index 100% rename from src/glpk/glpscl.c rename to src/glpk/draft/glpscl.c diff --git a/src/glpk/glpssx.h b/src/glpk/draft/glpssx.h similarity index 97% rename from src/glpk/glpssx.h rename to src/glpk/draft/glpssx.h index 1b55d7a..3b52b3c 100644 --- a/src/glpk/glpssx.h +++ b/src/glpk/draft/glpssx.h @@ -4,7 +4,7 @@ * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -27,6 +27,9 @@ #include "bfx.h" #include "env.h" +#if 1 /* 25/XI-2017 */ +#include "glpk.h" +#endif typedef struct SSX SSX; @@ -305,6 +308,14 @@ struct SSX /* actual change of xN[q] in the adjacent basis (it has the same sign as q_dir) */ /*--------------------------------------------------------------------*/ +#if 1 /* 25/XI-2017 */ + int msg_lev; + /* verbosity level: + GLP_MSG_OFF no output + GLP_MSG_ERR report errors and warnings + GLP_MSG_ON normal output + GLP_MSG_ALL highest verbosity */ +#endif int it_lim; /* simplex iterations limit; if this value is positive, it is decreased by one each time when one simplex iteration has been diff --git a/src/glpk/glpios.h b/src/glpk/draft/ios.h similarity index 97% rename from src/glpk/glpios.h rename to src/glpk/draft/ios.h index 8db97df..1cb07ee 100644 --- a/src/glpk/glpios.h +++ b/src/glpk/draft/ios.h @@ -1,10 +1,10 @@ -/* glpios.h (integer optimization suite) */ +/* ios.h (integer optimization suite) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* 2009, 2010, 2011, 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -22,19 +22,32 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifndef GLPIOS_H -#define GLPIOS_H +#ifndef IOS_H +#define IOS_H #include "prob.h" +#if 1 /* 02/II-2018 */ +#define NEW_LOCAL 1 +#endif + +#if 1 /* 15/II-2018 */ +#define NEW_COVER 1 +#endif + typedef struct IOSLOT IOSLOT; typedef struct IOSNPD IOSNPD; typedef struct IOSBND IOSBND; typedef struct IOSTAT IOSTAT; typedef struct IOSROW IOSROW; typedef struct IOSAIJ IOSAIJ; +#ifdef NEW_LOCAL /* 02/II-2018 */ +typedef glp_prob IOSPOOL; +typedef GLPROW IOSCUT; +#else typedef struct IOSPOOL IOSPOOL; typedef struct IOSCUT IOSCUT; +#endif struct glp_tree { /* branch-and-bound tree */ @@ -148,19 +161,14 @@ struct glp_tree /* built-in cut generators segment */ IOSPOOL *local; /* local cut pool */ -#if 0 /* 06/III-2016 */ - void *mir_gen; -#else - glp_mir *mir_gen; +#if 1 /* 13/II-2018 */ + glp_cov *cov_gen; + /* pointer to working area used by the cover cut generator */ #endif + glp_mir *mir_gen; /* pointer to working area used by the MIR cut generator */ -#if 0 /* 08/III-2016 */ - void *clq_gen; - /* pointer to working area used by the clique cut generator */ -#else glp_cfg *clq_gen; /* pointer to conflict graph used by the clique cut generator */ -#endif /*--------------------------------------------------------------*/ void *pcost; /* pointer to working area used on pseudocost branching */ @@ -172,18 +180,10 @@ struct glp_tree /* control parameters and statistics */ const glp_iocp *parm; /* copy of control parameters passed to the solver */ -#if 0 /* 10/VI-2013 */ - glp_long tm_beg; -#else double tm_beg; -#endif /* starting time of the search, in seconds; the total time of the search is the difference between xtime() and tm_beg */ -#if 0 /* 10/VI-2013 */ - glp_long tm_lag; -#else double tm_lag; -#endif /* the most recent time, in seconds, at which the progress of the the search was displayed */ int sol_cnt; @@ -373,6 +373,7 @@ struct IOSAIJ /* pointer to next coefficient for the same row */ }; +#ifndef NEW_LOCAL /* 02/II-2018 */ struct IOSPOOL { /* cut pool */ int size; @@ -386,7 +387,9 @@ struct IOSPOOL IOSCUT *curr; /* pointer to the current cut */ }; +#endif +#ifndef NEW_LOCAL /* 02/II-2018 */ struct IOSCUT { /* cut (cutting plane constraint) */ char *name; @@ -407,6 +410,7 @@ struct IOSCUT IOSCUT *next; /* pointer to next cut */ }; +#endif #define ios_create_tree _glp_ios_create_tree glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm); diff --git a/src/glpk/dune b/src/glpk/dune index 400ef0b..b2ec7ba 100644 --- a/src/glpk/dune +++ b/src/glpk/dune @@ -3,6 +3,7 @@ (public_name mccs.glpk.internal) (modules ocaml_mccs_glpk) (c_names ; Core GLPK + ; draft bfd glpapi06 glpapi09 @@ -14,14 +15,8 @@ glpios03 glpios07 glpios09 - glpios10 glpios11 glpios12 - glpnpp01 - glpnpp02 - glpnpp03 - glpnpp04 - glpnpp05 glpscl ; api advbas @@ -30,6 +25,7 @@ prob2 prob4 prob5 + wrmip ; bflib btf btfint @@ -43,11 +39,13 @@ scfint sgf sva - ; cglib + ; intopt cfg cfg1 cfg2 clqcut + covgen + fpump gmicut gmigen mirgen @@ -59,6 +57,7 @@ error stdc stdout + stream time tls ; misc @@ -66,14 +65,22 @@ dmp gcd jd + ks mc13d mc21a + mt1 rng rng1 round2n triang wclique wclique1 + ; npp + npp1 + npp2 + npp3 + npp4 + npp5 ; proxy proxy proxy1 @@ -97,12 +104,16 @@ (copy_files# bflib/*.{c,h}) -(copy_files# cglib/*.{c,h}) +(copy_files# draft/*.{c,h}) (copy_files# env/*.{c,h}) +(copy_files# intopt/*.{c,h}) + (copy_files# misc/*.{c,h}) +(copy_files# npp/*.{c,h}) + (copy_files# proxy/*.{c,h}) (copy_files# simplex/*.{c,h}) diff --git a/src/glpk/env/env.h b/src/glpk/env/env.h index 4b839c8..67214ef 100644 --- a/src/glpk/env/env.h +++ b/src/glpk/env/env.h @@ -29,10 +29,7 @@ typedef struct ENV ENV; typedef struct MBD MBD; -#ifndef SIZE_T_MAX #define SIZE_T_MAX (~(size_t)0) -#endif - /* largest value of size_t type */ #define TBUF_SIZE 4096 diff --git a/src/glpk/env/stdc.h b/src/glpk/env/stdc.h index 49988ef..793ddc6 100644 --- a/src/glpk/env/stdc.h +++ b/src/glpk/env/stdc.h @@ -57,6 +57,17 @@ char *xstrerr(int); char *xstrtok(char *, const char *); #endif +#if 1 /* 06/II-2018 */ +#ifdef HAVE_CONFIG_H +#include +#endif +#ifndef _MSC_VER +#define CDECL +#else +#define CDECL __cdecl +#endif +#endif + #endif /* eof */ diff --git a/src/glpk/env/stream.c b/src/glpk/env/stream.c new file mode 100644 index 0000000..906e5b0 --- /dev/null +++ b/src/glpk/env/stream.c @@ -0,0 +1,517 @@ +/* stream.c (stream input/output) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2008-2017 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "env.h" +#include "zlib.h" + +struct glp_file +{ /* sequential stream descriptor */ + char *base; + /* pointer to buffer */ + int size; + /* size of buffer, in bytes */ + char *ptr; + /* pointer to next byte in buffer */ + int cnt; + /* count of bytes in buffer */ + int flag; + /* stream flags: */ +#define IONULL 0x01 /* null file */ +#define IOSTD 0x02 /* standard stream */ +#define IOGZIP 0x04 /* gzipped file */ +#define IOWRT 0x08 /* output stream */ +#define IOEOF 0x10 /* end of file */ +#define IOERR 0x20 /* input/output error */ + void *file; + /* pointer to underlying control object */ +}; + +/*********************************************************************** +* NAME +* +* glp_open - open stream +* +* SYNOPSIS +* +* glp_file *glp_open(const char *name, const char *mode); +* +* DESCRIPTION +* +* The routine glp_open opens a file whose name is a string pointed to +* by name and associates a stream with it. +* +* The following special filenames are recognized by the routine (this +* feature is platform independent): +* +* "/dev/null" empty (null) file; +* "/dev/stdin" standard input stream; +* "/dev/stdout" standard output stream; +* "/dev/stderr" standard error stream. +* +* If the specified filename is ended with ".gz", it is assumed that +* the file is in gzipped format. In this case the file is compressed +* or decompressed by the I/O routines "on the fly". +* +* The parameter mode points to a string, which indicates the open mode +* and should be one of the following: +* +* "r" open text file for reading; +* "w" truncate to zero length or create text file for writing; +* "a" append, open or create text file for writing at end-of-file; +* "rb" open binary file for reading; +* "wb" truncate to zero length or create binary file for writing; +* "ab" append, open or create binary file for writing at end-of-file. +* +* RETURNS +* +* The routine glp_open returns a pointer to the object controlling the +* stream. If the operation fails, the routine returns NULL. */ + +glp_file *glp_open(const char *name, const char *mode) +{ glp_file *f; + int flag; + void *file; + if (strcmp(mode, "r") == 0 || strcmp(mode, "rb") == 0) + flag = 0; + else if (strcmp(mode, "w") == 0 || strcmp(mode, "wb") == 0) + flag = IOWRT; +#if 1 /* 08/V-2014 */ + else if (strcmp(mode, "a") == 0 || strcmp(mode, "ab") == 0) + flag = IOWRT; +#endif + else + xerror("glp_open: invalid mode string\n"); + if (strcmp(name, "/dev/null") == 0) + { flag |= IONULL; + file = NULL; + } + else if (strcmp(name, "/dev/stdin") == 0) + { flag |= IOSTD; + file = stdin; + } + else if (strcmp(name, "/dev/stdout") == 0) + { flag |= IOSTD; + file = stdout; + } + else if (strcmp(name, "/dev/stderr") == 0) + { flag |= IOSTD; + file = stderr; + } + else + { char *ext = strrchr(name, '.'); + if (ext == NULL || strcmp(ext, ".gz") != 0) + { file = fopen(name, mode); + if (file == NULL) +#if 0 /* 29/I-2017 */ + { put_err_msg(strerror(errno)); +#else + { put_err_msg(xstrerr(errno)); +#endif + return NULL; + } + } + else + { flag |= IOGZIP; + if (strcmp(mode, "r") == 0) + mode = "rb"; + else if (strcmp(mode, "w") == 0) + mode = "wb"; +#if 1 /* 08/V-2014; this mode seems not to work */ + else if (strcmp(mode, "a") == 0) + mode = "ab"; +#endif + file = gzopen(name, mode); + if (file == NULL) +#if 0 /* 29/I-2017 */ + { put_err_msg(strerror(errno)); +#else + { put_err_msg(xstrerr(errno)); +#endif + return NULL; + } + } + } + f = talloc(1, glp_file); + f->base = talloc(BUFSIZ, char); + f->size = BUFSIZ; + f->ptr = f->base; + f->cnt = 0; + f->flag = flag; + f->file = file; + return f; +} + +/*********************************************************************** +* NAME +* +* glp_eof - test end-of-file indicator +* +* SYNOPSIS +* +* int glp_eof(glp_file *f); +* +* DESCRIPTION +* +* The routine glp_eof tests the end-of-file indicator for the stream +* pointed to by f. +* +* RETURNS +* +* The routine glp_eof returns non-zero if and only if the end-of-file +* indicator is set for the specified stream. */ + +int glp_eof(glp_file *f) +{ return + f->flag & IOEOF; +} + +/*********************************************************************** +* NAME +* +* glp_ioerr - test I/O error indicator +* +* SYNOPSIS +* +* int glp_ioerr(glp_file *f); +* +* DESCRIPTION +* +* The routine glp_ioerr tests the I/O error indicator for the stream +* pointed to by f. +* +* RETURNS +* +* The routine glp_ioerr returns non-zero if and only if the I/O error +* indicator is set for the specified stream. */ + +int glp_ioerr(glp_file *f) +{ return + f->flag & IOERR; +} + +/*********************************************************************** +* NAME +* +* glp_read - read data from stream +* +* SYNOPSIS +* +* int glp_read(glp_file *f, void *buf, int nnn); +* +* DESCRIPTION +* +* The routine glp_read reads, into the buffer pointed to by buf, up to +* nnn bytes, from the stream pointed to by f. +* +* RETURNS +* +* The routine glp_read returns the number of bytes successfully read +* (which may be less than nnn). If an end-of-file is encountered, the +* end-of-file indicator for the stream is set and glp_read returns +* zero. If a read error occurs, the error indicator for the stream is +* set and glp_read returns a negative value. */ + +int glp_read(glp_file *f, void *buf, int nnn) +{ int nrd, cnt; + if (f->flag & IOWRT) + xerror("glp_read: attempt to read from output stream\n"); + if (nnn < 1) + xerror("glp_read: nnn = %d; invalid parameter\n", nnn); + for (nrd = 0; nrd < nnn; nrd += cnt) + { if (f->cnt == 0) + { /* buffer is empty; fill it */ + if (f->flag & IONULL) + cnt = 0; + else if (!(f->flag & IOGZIP)) + { cnt = fread(f->base, 1, f->size, (FILE *)(f->file)); + if (ferror((FILE *)(f->file))) + { f->flag |= IOERR; +#if 0 /* 29/I-2017 */ + put_err_msg(strerror(errno)); +#else + put_err_msg(xstrerr(errno)); +#endif + return EOF; + } + } + else + { int errnum; + const char *msg; + cnt = gzread((gzFile)(f->file), f->base, f->size); + if (cnt < 0) + { f->flag |= IOERR; + msg = gzerror((gzFile)(f->file), &errnum); + if (errnum == Z_ERRNO) +#if 0 /* 29/I-2017 */ + put_err_msg(strerror(errno)); +#else + put_err_msg(xstrerr(errno)); +#endif + else + put_err_msg(msg); + return EOF; + } + } + if (cnt == 0) + { if (nrd == 0) + f->flag |= IOEOF; + break; + } + f->ptr = f->base; + f->cnt = cnt; + } + cnt = nnn - nrd; + if (cnt > f->cnt) + cnt = f->cnt; + memcpy((char *)buf + nrd, f->ptr, cnt); + f->ptr += cnt; + f->cnt -= cnt; + } + return nrd; +} + +/*********************************************************************** +* NAME +* +* glp_getc - read character from stream +* +* SYNOPSIS +* +* int glp_getc(glp_file *f); +* +* DESCRIPTION +* +* The routine glp_getc obtains a next character as an unsigned char +* converted to an int from the input stream pointed to by f. +* +* RETURNS +* +* The routine glp_getc returns the next character obtained. However, +* if an end-of-file is encountered or a read error occurs, the routine +* returns EOF. (An end-of-file and a read error can be distinguished +* by use of the routines glp_eof and glp_ioerr.) */ + +int glp_getc(glp_file *f) +{ unsigned char buf[1]; + if (f->flag & IOWRT) + xerror("glp_getc: attempt to read from output stream\n"); + if (glp_read(f, buf, 1) != 1) + return EOF; + return buf[0]; +} + +/*********************************************************************** +* do_flush - flush output stream +* +* This routine causes buffered data for the specified output stream to +* be written to the associated file. +* +* If the operation was successful, the routine returns zero, otherwise +* non-zero. */ + +static int do_flush(glp_file *f) +{ xassert(f->flag & IOWRT); + if (f->cnt > 0) + { if (f->flag & IONULL) + ; + else if (!(f->flag & IOGZIP)) + { if ((int)fwrite(f->base, 1, f->cnt, (FILE *)(f->file)) + != f->cnt) + { f->flag |= IOERR; +#if 0 /* 29/I-2017 */ + put_err_msg(strerror(errno)); +#else + put_err_msg(xstrerr(errno)); +#endif + return EOF; + } + } + else + { int errnum; + const char *msg; + if (gzwrite((gzFile)(f->file), f->base, f->cnt) != f->cnt) + { f->flag |= IOERR; + msg = gzerror((gzFile)(f->file), &errnum); + if (errnum == Z_ERRNO) +#if 0 /* 29/I-2017 */ + put_err_msg(strerror(errno)); +#else + put_err_msg(xstrerr(errno)); +#endif + else + put_err_msg(msg); + return EOF; + } + } + } + f->ptr = f->base; + f->cnt = 0; + return 0; +} + +/*********************************************************************** +* NAME +* +* glp_write - write data to stream +* +* SYNOPSIS +* +* int glp_write(glp_file *f, const void *buf, int nnn); +* +* DESCRIPTION +* +* The routine glp_write writes, from the buffer pointed to by buf, up +* to nnn bytes, to the stream pointed to by f. +* +* RETURNS +* +* The routine glp_write returns the number of bytes successfully +* written (which is equal to nnn). If a write error occurs, the error +* indicator for the stream is set and glp_write returns a negative +* value. */ + +int glp_write(glp_file *f, const void *buf, int nnn) +{ int nwr, cnt; + if (!(f->flag & IOWRT)) + xerror("glp_write: attempt to write to input stream\n"); + if (nnn < 1) + xerror("glp_write: nnn = %d; invalid parameter\n", nnn); + for (nwr = 0; nwr < nnn; nwr += cnt) + { cnt = nnn - nwr; + if (cnt > f->size - f->cnt) + cnt = f->size - f->cnt; + memcpy(f->ptr, (const char *)buf + nwr, cnt); + f->ptr += cnt; + f->cnt += cnt; + if (f->cnt == f->size) + { /* buffer is full; flush it */ + if (do_flush(f) != 0) + return EOF; + } + } + return nwr; +} + +/*********************************************************************** +* NAME +* +* glp_format - write formatted data to stream +* +* SYNOPSIS +* +* int glp_format(glp_file *f, const char *fmt, ...); +* +* DESCRIPTION +* +* The routine glp_format writes formatted data to the stream pointed +* to by f. The format control string pointed to by fmt specifies how +* subsequent arguments are converted for output. +* +* RETURNS +* +* The routine glp_format returns the number of characters written, or +* a negative value if an output error occurs. */ + +int glp_format(glp_file *f, const char *fmt, ...) +{ ENV *env = get_env_ptr(); + va_list arg; + int nnn; + if (!(f->flag & IOWRT)) + xerror("glp_format: attempt to write to input stream\n"); + va_start(arg, fmt); + nnn = vsprintf(env->term_buf, fmt, arg); + xassert(0 <= nnn && nnn < TBUF_SIZE); + va_end(arg); + return nnn == 0 ? 0 : glp_write(f, env->term_buf, nnn); +} + +/*********************************************************************** +* NAME +* +* glp_close - close stream +* +* SYNOPSIS +* +* int glp_close(glp_file *f); +* +* DESCRIPTION +* +* The routine glp_close closes the stream pointed to by f. +* +* RETURNS +* +* If the operation was successful, the routine returns zero, otherwise +* non-zero. */ + +int glp_close(glp_file *f) +{ int ret = 0; + if (f->flag & IOWRT) + { if (do_flush(f) != 0) + ret = EOF; + } + if (f->flag & (IONULL | IOSTD)) + ; + else if (!(f->flag & IOGZIP)) + { if (fclose((FILE *)(f->file)) != 0) + { if (ret == 0) +#if 0 /* 29/I-2017 */ + { put_err_msg(strerror(errno)); +#else + { put_err_msg(xstrerr(errno)); +#endif + ret = EOF; + } + } + } + else + { int errnum; + errnum = gzclose((gzFile)(f->file)); + if (errnum == Z_OK) + ; + else if (errnum == Z_ERRNO) + { if (ret == 0) +#if 0 /* 29/I-2017 */ + { put_err_msg(strerror(errno)); +#else + { put_err_msg(xstrerr(errno)); +#endif + ret = EOF; + } + } +#if 1 /* FIXME */ + else + { if (ret == 0) + { ENV *env = get_env_ptr(); + sprintf(env->term_buf, "gzclose returned %d", errnum); + put_err_msg(env->term_buf); + ret = EOF; + } + } +#endif + } + tfree(f->base); + tfree(f); + return ret; +} + +/* eof */ diff --git a/src/glpk/glpk.h b/src/glpk/glpk.h index 83fa3f1..f4e250f 100644 --- a/src/glpk/glpk.h +++ b/src/glpk/glpk.h @@ -1,12 +1,11 @@ -/* glpk.h (GLPK API) */ +/* glpk.h */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013, 2014, 2015, 2016, 2017 Andrew Makhorin, -* Department for Applied Informatics, Moscow Aviation Institute, -* Moscow, Russia. All rights reserved. E-mail: . +* Copyright (C) 2000-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . * * GLPK is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by @@ -34,7 +33,7 @@ extern "C" { /* library version numbers: */ #define GLP_MAJOR_VERSION 4 -#define GLP_MINOR_VERSION 63 +#define GLP_MINOR_VERSION 65 typedef struct glp_prob glp_prob; /* LP/MIP problem object */ @@ -107,7 +106,7 @@ typedef struct } glp_bfcp; typedef struct -{ /* simplex method control parameters */ +{ /* simplex solver control parameters */ int msg_lev; /* message level: */ #define GLP_MSG_OFF 0 /* no output */ #define GLP_MSG_ERR 1 /* warning and error messages only */ @@ -297,6 +296,11 @@ typedef struct /* (reserved for use in the future) */ } glp_cpxcp; +#if 1 /* 10/XII-2017 */ +typedef struct glp_prep glp_prep; +/* LP/MIP preprocessor workspace */ +#endif + typedef struct glp_tran glp_tran; /* MathProg translator workspace */ @@ -666,6 +670,30 @@ void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1, double *value1, double *coef2, int *var2, double *value2); /* analyze objective coefficient at basic variable */ +#if 1 /* 10/XII-2017 */ +glp_prep *glp_npp_alloc_wksp(void); +/* allocate the preprocessor workspace */ + +void glp_npp_load_prob(glp_prep *prep, glp_prob *P, int sol, + int names); +/* load original problem instance */ + +int glp_npp_preprocess1(glp_prep *prep, int hard); +/* perform basic LP/MIP preprocessing */ + +void glp_npp_build_prob(glp_prep *prep, glp_prob *Q); +/* build resultant problem instance */ + +void glp_npp_postprocess(glp_prep *prep, glp_prob *Q); +/* postprocess solution to resultant problem */ + +void glp_npp_obtain_sol(glp_prep *prep, glp_prob *P); +/* obtain solution to original problem */ + +void glp_npp_free_wksp(glp_prep *prep); +/* free the preprocessor workspace */ +#endif + int glp_ios_reason(glp_tree *T); /* determine reason for calling the callback routine */ @@ -739,52 +767,46 @@ void glp_ios_terminate(glp_tree *T); int glp_gmi_cut(glp_prob *P, int j, int ind[], double val[], double phi[]); /* generate Gomory's mixed integer cut (core routine) */ -#endif -#ifdef GLP_UNDOC int glp_gmi_gen(glp_prob *P, glp_prob *pool, int max_cuts); /* generate Gomory's mixed integer cuts */ -#endif -#ifdef GLP_UNDOC +typedef struct glp_cov glp_cov; +/* cover cur generator workspace */ + +glp_cov *glp_cov_init(glp_prob *P); +/* create and initialize cover cut generator */ + +void glp_cov_gen1(glp_prob *P, glp_cov *cov, glp_prob *pool); +/* generate locally valid simple cover cuts */ + +void glp_cov_free(glp_cov *cov); +/* delete cover cut generator workspace */ + typedef struct glp_mir glp_mir; /* MIR cut generator workspace */ -#endif -#ifdef GLP_UNDOC glp_mir *glp_mir_init(glp_prob *P); /* create and initialize MIR cut generator */ -#endif -#ifdef GLP_UNDOC int glp_mir_gen(glp_prob *P, glp_mir *mir, glp_prob *pool); /* generate mixed integer rounding (MIR) cuts */ -#endif -#ifdef GLP_UNDOC void glp_mir_free(glp_mir *mir); /* delete MIR cut generator workspace */ -#endif -#ifdef GLP_UNDOC typedef struct glp_cfg glp_cfg; /* conflict graph descriptor */ -#endif -#ifdef GLP_UNDOC glp_cfg *glp_cfg_init(glp_prob *P); /* create and initialize conflict graph */ -#endif -#ifdef GLP_UNDOC void glp_cfg_free(glp_cfg *G); /* delete conflict graph descriptor */ -#endif -#ifdef GLP_UNDOC int glp_clq_cut(glp_prob *P, glp_cfg *G, int ind[], double val[]); /* generate clique cut from conflict graph */ -#endif +#endif /* GLP_UNDOC */ void glp_init_mpscp(glp_mpscp *parm); /* initialize MPS format control parameters */ diff --git a/src/glpk/cglib/cfg.c b/src/glpk/intopt/cfg.c similarity index 100% rename from src/glpk/cglib/cfg.c rename to src/glpk/intopt/cfg.c diff --git a/src/glpk/cglib/cfg.h b/src/glpk/intopt/cfg.h similarity index 100% rename from src/glpk/cglib/cfg.h rename to src/glpk/intopt/cfg.h diff --git a/src/glpk/cglib/cfg1.c b/src/glpk/intopt/cfg1.c similarity index 99% rename from src/glpk/cglib/cfg1.c rename to src/glpk/intopt/cfg1.c index cc2b6df..80a2e83 100644 --- a/src/glpk/cglib/cfg1.c +++ b/src/glpk/intopt/cfg1.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2012-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2012-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "cfg.h" #include "env.h" #include "prob.h" @@ -201,11 +197,7 @@ struct term { int ind; double val; }; /* term a[j] * z[j] used to sort a[j]'s */ -#ifndef _MSC_VER -static int fcmp(const void *e1, const void *e2) -#else -static int __cdecl fcmp(const void *e1, const void *e2) -#endif +static int CDECL fcmp(const void *e1, const void *e2) { /* auxiliary routine called from qsort */ const struct term *t1 = e1, *t2 = e2; if (t1->val > t2->val) diff --git a/src/glpk/cglib/cfg2.c b/src/glpk/intopt/cfg2.c similarity index 100% rename from src/glpk/cglib/cfg2.c rename to src/glpk/intopt/cfg2.c diff --git a/src/glpk/cglib/clqcut.c b/src/glpk/intopt/clqcut.c similarity index 100% rename from src/glpk/cglib/clqcut.c rename to src/glpk/intopt/clqcut.c diff --git a/src/glpk/intopt/covgen.c b/src/glpk/intopt/covgen.c new file mode 100644 index 0000000..427c3aa --- /dev/null +++ b/src/glpk/intopt/covgen.c @@ -0,0 +1,885 @@ +/* covgen.c */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2017-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "env.h" +#include "fvs.h" +#include "ks.h" +#include "prob.h" + +struct glp_cov +{ /* cover cut generator working area */ + int n; + /* number of columns (variables) */ + glp_prob *set; + /* set of globally valid 0-1 knapsack inequalities chosen from + * the root problem; each inequality is either original row or + * its relaxation (surrogate 0-1 knapsack) which is constructed + * by substitution of lower/upper single/variable bounds for + * continuous and general integer (non-binary) variables */ +}; + +struct bnd +{ /* simple or variable bound */ + /* if z = 0, it is a simple bound x >= or <= b; if b = -DBL_MAX + * (b = +DBL_MAX), x has no lower (upper) bound; otherwise, if + * z != 0, it is a variable bound x >= or <= a * z + b */ + int z; + /* number of binary variable or 0 */ + double a, b; + /* bound parameters */ +}; + +struct csa +{ /* common storage area */ + glp_prob *P; + /* original (root) MIP */ + struct bnd *l; /* struct bnd l[1+P->n]; */ + /* lower simple/variable bounds of variables */ + struct bnd *u; /* struct bnd u[1+P->n]; */ + /* upper simple/variable bounds of variables */ + glp_prob *set; + /* see struct glp_cov above */ +}; + +/*********************************************************************** +* init_bounds - initialize bounds of variables with simple bounds +* +* This routine initializes lower and upper bounds of all variables +* with simple bounds specified in the original mip. */ + +static void init_bounds(struct csa *csa) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + int j; + for (j = 1; j <= P->n; j++) + { l[j].z = u[j].z = 0; + l[j].a = u[j].a = 0; + l[j].b = glp_get_col_lb(P, j); + u[j].b = glp_get_col_ub(P, j); + } + return; +} + +/*********************************************************************** +* check_vb - check variable bound +* +* This routine checks if the specified i-th row has the form +* +* a1 * x + a2 * z >= or <= rhs, (1) +* +* where x is a non-fixed continuous or general integer variable, and +* z is a binary variable. If it is, the routine converts the row to +* the following variable lower/upper bound (VLB/VUB) of x: +* +* x >= or <= a * z + b, (2) +* +* where a = - a2 / a1, b = rhs / a1. Note that the inequality type is +* changed to opposite one when a1 < 0. +* +* If the row is identified as a variable bound, the routine returns +* GLP_LO for VLB or GLP_UP for VUB and provides the reference numbers +* of variables x and z and values of a and b. Otherwise, the routine +* returns zero. */ + +static int check_vb(struct csa *csa, int i, int *x, int *z, double *a, + double *b) +{ glp_prob *P = csa->P; + GLPROW *row; + GLPAIJ *a1, *a2; + int type; + double rhs; + xassert(1 <= i && i <= P->m); + row = P->row[i]; + /* check row type */ + switch (row->type) + { case GLP_LO: + case GLP_UP: + break; + default: + return 0; + } + /* take first term of the row */ + a1 = row->ptr; + if (a1 == NULL) + return 0; + /* take second term of the row */ + a2 = a1->r_next; + if (a2 == NULL) + return 0; + /* there should be exactly two terms in the row */ + if (a2->r_next != NULL) + return 0; + /* if first term is a binary variable, swap the terms */ + if (glp_get_col_kind(P, a1->col->j) == GLP_BV) + { GLPAIJ *a; + a = a1, a1 = a2, a2 = a; + } + /* now first term should be a non-fixed continuous or general + * integer variable */ + if (a1->col->type == GLP_FX) + return 0; + if (glp_get_col_kind(P, a1->col->j) == GLP_BV) + return 0; + /* and second term should be a binary variable */ + if (glp_get_col_kind(P, a2->col->j) != GLP_BV) + return 0; + /* VLB/VUB row has been identified */ + switch (row->type) + { case GLP_LO: + type = a1->val > 0 ? GLP_LO : GLP_UP; + rhs = row->lb; + break; + case GLP_UP: + type = a1->val > 0 ? GLP_UP : GLP_LO; + rhs = row->ub; + break; + default: + xassert(type != type); + } + *x = a1->col->j; + *z = a2->col->j; + *a = - a2->val / a1->val; + *b = rhs / a1->val; + return type; +} + +/*********************************************************************** +* set_vb - set variable bound +* +* This routine sets lower or upper variable bound specified as +* +* x >= a * z + b (type = GLP_LO) +* +* x <= a * z + b (type = GLP_UP) */ + +static void set_vb(struct csa *csa, int type, int x, int z, double a, + double b) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + xassert(glp_get_col_type(P, x) != GLP_FX); + xassert(glp_get_col_kind(P, x) != GLP_BV); + xassert(glp_get_col_kind(P, z) == GLP_BV); + xassert(a != 0); + switch (type) + { case GLP_LO: + /* FIXME: check existing simple lower bound? */ + l[x].z = z, l[x].a = a, l[x].b = b; + break; + case GLP_UP: + /* FIXME: check existing simple upper bound? */ + u[x].z = z, u[x].a = a, u[x].b = b; + break; + default: + xassert(type != type); + } + return; +} + +/*********************************************************************** +* obtain_vbs - obtain and set variable bounds +* +* This routine walks thru all rows of the original mip, identifies +* rows specifying variable lower/upper bounds, and sets these bounds +* for corresponding (non-binary) variables. */ + +static void obtain_vbs(struct csa *csa) +{ glp_prob *P = csa->P; + int i, x, z, type, save; + double a, b; + for (i = 1; i <= P->m; i++) + { switch (P->row[i]->type) + { case GLP_FR: + break; + case GLP_LO: + case GLP_UP: + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + break; + case GLP_DB: + case GLP_FX: + /* double-side inequality l <= ... <= u and equality + * ... = l = u are considered as two single inequalities + * ... >= l and ... <= u */ + save = P->row[i]->type; + P->row[i]->type = GLP_LO; + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + P->row[i]->type = GLP_UP; + type = check_vb(csa, i, &x, &z, &a, &b); + if (type) + set_vb(csa, type, x, z, a, b); + P->row[i]->type = save; + break; + default: + xassert(P != P); + } + } + return; +} + +/*********************************************************************** +* add_term - add term to sparse vector +* +* This routine computes the following linear combination: +* +* v := v + a * e[j], +* +* where v is a sparse vector in full storage format, a is a non-zero +* scalar, e[j] is j-th column of unity matrix. */ + +static void add_term(FVS *v, int j, double a) +{ xassert(1 <= j && j <= v->n); + xassert(a != 0); + if (v->vec[j] == 0) + { /* create j-th component */ + v->nnz++; + xassert(v->nnz <= v->n); + v->ind[v->nnz] = j; + } + /* perform addition */ + v->vec[j] += a; + if (fabs(v->vec[j]) < 1e-9 * (1 + fabs(a))) + { /* remove j-th component */ + v->vec[j] = DBL_MIN; + } + return; +} + +/*********************************************************************** +* build_ks - build "0-1 knapsack" inequality +* +* Given an inequality of "not greater" type: +* +* sum{j in 1..n} a[j]*x[j] <= b, (1) +* +* this routine attempts to transform it to equivalent or relaxed "0-1 +* knapsack" inequality that contains only binary variables. +* +* If x[j] is a binary variable, the term a[j]*x[j] is not changed. +* Otherwise, if x[j] is a continuous or integer non-binary variable, +* it is replaced by its lower (if a[j] > 0) or upper (if a[j] < 0) +* single or variable bound. In the latter case, if x[j] is a non-fixed +* variable, this results in a relaxation of original inequality known +* as "surrogate knapsack". Thus, if the specified inequality is valid +* for the original mip, the resulting inequality is also valid. +* +* Note that in both source and resulting inequalities coefficients +* a[j] can have any sign. +* +* On entry to the routine the source inequality is specified by the +* parameters n, ind (contains original numbers of x[j]), a, and b. The +* parameter v is a working sparse vector whose components are assumed +* to be zero. +* +* On exit the routine stores the resulting "0-1 knapsack" inequality +* in the parameters ind, a, and b, and returns n which is the number +* of terms in the resulting inequality. Zero content of the vector v +* is restored before exit. +* +* If the resulting inequality cannot be constructed due to missing +* lower/upper bounds of some variable, the routine returns a negative +* value. */ + +static int build_ks(struct csa *csa, int n, int ind[], double a[], + double *b, FVS *v) +{ glp_prob *P = csa->P; + struct bnd *l = csa->l, *u = csa->u; + int j, k; + /* check that v = 0 */ +#ifdef GLP_DEBUG + fvs_check_vec(v); +#endif + xassert(v->nnz == 0); + /* walk thru terms of original inequality */ + for (j = 1; j <= n; j++) + { /* process term a[j]*x[j] */ + k = ind[j]; /* original number of x[j] in mip */ + if (glp_get_col_kind(P, k) == GLP_BV) + { /* x[j] is a binary variable */ + /* include its term into resulting inequality */ + add_term(v, k, a[j]); + } + else if (a[j] > 0) + { /* substitute x[j] by its lower bound */ + if (l[k].b == -DBL_MAX) + { /* x[j] has no lower bound */ + n = -1; + goto skip; + } + else if (l[k].z == 0) + { /* x[j] has simple lower bound */ + *b -= a[j] * l[k].b; + } + else + { /* x[j] has variable lower bound (a * z + b) */ + add_term(v, l[k].z, a[j] * l[k].a); + *b -= a[j] * l[k].b; + } + } + else /* a[j] < 0 */ + { /* substitute x[j] by its upper bound */ + if (u[k].b == +DBL_MAX) + { /* x[j] has no upper bound */ + n = -1; + goto skip; + } + else if (u[k].z == 0) + { /* x[j] has simple upper bound */ + *b -= a[j] * u[k].b; + } + else + { /* x[j] has variable upper bound (a * z + b) */ + add_term(v, u[k].z, a[j] * u[k].a); + *b -= a[j] * u[k].b; + } + } + } + /* replace tiny coefficients by exact zeros (see add_term) */ + fvs_adjust_vec(v, 2 * DBL_MIN); + /* copy terms of resulting inequality */ + xassert(v->nnz <= n); + n = v->nnz; + for (j = 1; j <= n; j++) + { ind[j] = v->ind[j]; + a[j] = v->vec[ind[j]]; + } +skip: /* restore zero content of v */ + fvs_clear_vec(v); + return n; +} + +/*********************************************************************** +* can_be_active - check if inequality can be active +* +* This routine checks if the specified "0-1 knapsack" inequality +* +* sum{j in 1..n} a[j]*x[j] <= b +* +* can be active. If so, the routine returns true, otherwise false. */ + +static int can_be_active(int n, const double a[], double b) +{ int j; + double s; + s = 0; + for (j = 1; j <= n; j++) + { if (a[j] > 0) + s += a[j]; + } + return s > b + .001 * (1 + fabs(b)); +} + +/*********************************************************************** +* is_sos_ineq - check if inequality is packing (SOS) constraint +* +* This routine checks if the specified "0-1 knapsack" inequality +* +* sum{j in 1..n} a[j]*x[j] <= b (1) +* +* is equivalent to packing inequality (Padberg calls such inequalities +* special ordered set or SOS constraints) +* +* sum{j in J'} x[j] - sum{j in J"} x[j] <= 1 - |J"|. (2) +* +* If so, the routine returns true, otherwise false. +* +* Note that if X is a set of feasible binary points satisfying to (2), +* its convex hull conv(X) equals to the set of feasible points of LP +* relaxation of (2), which is a n-dimensional simplex, so inequalities +* (2) are useless for generating cover cuts (due to unimodularity). +* +* ALGORITHM +* +* First, we make all a[j] positive by complementing x[j] = 1 - x'[j] +* in (1). This is performed implicitly (i.e. actually the array a is +* not changed), but b is replaced by b - sum{j : a[j] < 0}. +* +* Then we find two smallest coefficients a[p] = min{j in 1..n} a[j] +* and a[q] = min{j in 1..n : j != p} a[j]. It is obvious that if +* a[p] + a[q] > b, then a[i] + a[j] > b for all i != j, from which it +* follows that x[i] + x[j] <= 1 for all i != j. But the latter means +* that the original inequality (with all a[j] > 0) is equivalent to +* packing inequality +* +* sum{j in 1..n} x[j] <= 1. (3) +* +* Returning to original (uncomplemented) variables x'[j] = 1 - x[j] +* we have that the original inequality is equivalent to (2), where +* J' = {j : a[j] > 0} and J" = {j : a[j] < 0}. */ + +static int is_sos_ineq(int n, const double a[], double b) +{ int j, p, q; + xassert(n >= 2); + /* compute b := b - sum{j : a[j] < 0} */ + for (j = 1; j <= n; j++) + { if (a[j] < 0) + b -= a[j]; + } + /* find a[p] = min{j in 1..n} a[j] */ + p = 1; + for (j = 2; j <= n; j++) + { if (fabs(a[p]) > fabs(a[j])) + p = j; + } + /* find a[q] = min{j in 1..n : j != p} a[j] */ + q = 0; + for (j = 1; j <= n; j++) + { if (j != p) + { if (q == 0 || fabs(a[q]) > fabs(a[j])) + q = j; + } + } + xassert(q != 0); + /* check condition a[p] + a[q] > b */ + return fabs(a[p]) + fabs(a[q]) > b + .001 * (1 + fabs(b)); +} + +/*********************************************************************** +* process_ineq - basic inequality processing +* +* This routine performs basic processing of an inequality of "not +* greater" type +* +* sum{j in 1..n} a[j]*x[j] <= b +* +* specified by the parameters, n, ind, a, and b. +* +* If the inequality can be transformed to "0-1 knapsack" ineqiality +* suitable for generating cover cuts, the routine adds it to the set +* of "0-1 knapsack" inequalities. +* +* Note that the arrays ind and a are not saved on exit. */ + +static void process_ineq(struct csa *csa, int n, int ind[], double a[], + double b, FVS *v) +{ int i; + /* attempt to transform the specified inequality to equivalent or + * relaxed "0-1 knapsack" inequality */ + n = build_ks(csa, n, ind, a, &b, v); + if (n <= 1) + { /* uninteresting inequality (in principle, such inequalities + * should be removed by the preprocessor) */ + goto done; + } + if (!can_be_active(n, a, b)) + { /* inequality is redundant (i.e. cannot be active) */ + goto done; + } + if (is_sos_ineq(n, a, b)) + { /* packing (SOS) inequality is useless for generating cover + * cuts; currently such inequalities are just ignored */ + goto done; + } + /* add resulting "0-1 knapsack" inequality to the set */ + i = glp_add_rows(csa->set, 1); + glp_set_mat_row(csa->set, i, n, ind, a); + glp_set_row_bnds(csa->set, i, GLP_UP, b, b); +done: return; +} + +/**********************************************************************/ + +glp_cov *glp_cov_init(glp_prob *P) +{ /* create and initialize cover cut generator */ + glp_cov *cov; + struct csa csa; + int i, k, len, *ind; + double rhs, *val; + FVS fvs; + csa.P = P; + csa.l = talloc(1+P->n, struct bnd); + csa.u = talloc(1+P->n, struct bnd); + csa.set = glp_create_prob(); + glp_add_cols(csa.set, P->n); + /* initialize bounds of variables with simple bounds */ + init_bounds(&csa); + /* obtain and set variable bounds */ + obtain_vbs(&csa); + /* allocate working arrays */ + ind = talloc(1+P->n, int); + val = talloc(1+P->n, double); + fvs_alloc_vec(&fvs, P->n); + /* process all rows of the root mip */ + for (i = 1; i <= P->m; i++) + { switch (P->row[i]->type) + { case GLP_FR: + break; + case GLP_LO: + /* obtain row of ">=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->lb; + /* transforms it to row of "<=" type */ + for (k = 1; k <= len; k++) + val[k] = - val[k]; + rhs = - rhs; + /* process the row */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + case GLP_UP: + /* obtain row of "<=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->ub; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + case GLP_DB: + case GLP_FX: + /* double-sided inequalitiy and equality constraints are + * processed as two separate inequalities */ + /* obtain row as if it were of ">=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->lb; + /* transforms it to row of "<=" type */ + for (k = 1; k <= len; k++) + val[k] = - val[k]; + rhs = - rhs; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + /* obtain the same row as if it were of "<=" type */ + len = glp_get_mat_row(P, i, ind, val); + rhs = P->row[i]->ub; + /* and process it */ + process_ineq(&csa, len, ind, val, rhs, &fvs); + break; + default: + xassert(P != P); + } + } + /* free working arrays */ + tfree(ind); + tfree(val); + fvs_check_vec(&fvs); + fvs_free_vec(&fvs); + /* the set of "0-1 knapsack" inequalities has been built */ + if (csa.set->m == 0) + { /* the set is empty */ + xprintf("No 0-1 knapsack inequalities detected\n"); + cov = NULL; + glp_delete_prob(csa.set); + } + else + { /* create the cover cut generator working area */ + xprintf("Number of 0-1 knapsack inequalities = %d\n", + csa.set->m); + cov = talloc(1, glp_cov); + cov->n = P->n; + cov->set = csa.set; +#if 0 + glp_write_lp(cov->set, 0, "set.lp"); +#endif + } + tfree(csa.l); + tfree(csa.u); + return cov; +} + +/*********************************************************************** +* solve_ks - solve 0-1 knapsack problem +* +* This routine finds (sub)optimal solution to 0-1 knapsack problem: +* +* maximize z = sum{j in 1..n} c[j]x[j] (1) +* +* s.t. sum{j in 1..n} a[j]x[j] <= b (2) +* +* x[j] in {0, 1} for all j in 1..n (3) +* +* It is assumed that the instance is non-normalized, i.e. parameters +* a, b, and c may have any sign. +* +* On exit the routine stores the (sub)optimal point found in locations +* x[1], ..., x[n] and returns the optimal objective value. However, if +* the instance is infeasible, the routine returns INT_MIN. */ + +static int solve_ks(int n, const int a[], int b, const int c[], + char x[]) +{ int z; + /* surprisingly, even for some small instances (n = 50-100) + * MT1 routine takes too much time, so it is used only for tiny + * instances */ + if (n <= 16) +#if 0 + z = ks_enum(n, a, b, c, x); +#else + z = ks_mt1(n, a, b, c, x); +#endif + else + z = ks_greedy(n, a, b, c, x); + return z; +} + +/*********************************************************************** +* simple_cover - find simple cover cut +* +* Given a 0-1 knapsack inequality (which may be globally as well as +* locally valid) +* +* sum{j in 1..n} a[j]x[j] <= b, (1) +* +* where all x[j] are binary variables and all a[j] are positive, and +* a fractional point x~{j in 1..n}, which is feasible to LP relaxation +* of (1), this routine attempts to find a simple cover inequality +* +* sum{j in C} (1 - x[j]) >= 1, (2) +* +* which is valid for (1) and violated at x~. +* +* Actually, the routine finds a cover C, i.e. a subset of {1, ..., n} +* such that +* +* sum{j in C} a[j] > b, (3) +* +* and which minimizes the left-hand side of (2) at x~ +* +* zeta = sum{j in C} (1 - x~[j]). (4) +* +* On exit the routine stores the characteritic vector z{j in 1..n} +* of the cover found (i.e. z[j] = 1 means j in C, and z[j] = 0 means +* j not in C), and returns corresponding minimal value of zeta (4). +* However, if no cover is found, the routine returns DBL_MAX. +* +* ALGORITHM +* +* The separation problem (3)-(4) is converted to 0-1 knapsack problem +* as follows. +* +* First, note that the constraint (3) is equivalent to +* +* sum{j in 1..n} a[j]z[j] >= b + eps, (5) +* +* where eps > 0 is a sufficiently small number (in case of integral +* a and b we may take eps = 1). Multiplying both sides of (5) by (-1) +* gives +* +* sum{j in 1..n} (-a[j])z[j] <= - b - eps. (6) +* +* To make all coefficients in (6) positive, z[j] is complemented by +* substitution z[j] = 1 - z'[j] that finally gives +* +* sum{j in 1..n} a[j]z'[j] <= sum{j in 1..n} a[j] - b - eps. (7) +* +* Minimization of zeta (4) is equivalent to maximization of +* +* -zeta = sum{j in 1..n} (x~[j] - 1)z[j]. (8) +* +* Substitution z[j] = 1 - z'[j] gives +* +* -zeta = sum{j in 1..n} (1 - x~[j])z'[j] - zeta0, (9) +* +* where zeta0 = sum{j in 1..n} (1 - x~[j]) is a constant term. +* +* Thus, the 0-1 knapsack problem to be solved is the following: +* +* maximize +* +* -zeta = sum{j in 1..n} (1 - x~[j])z'[j] - zeta0 (10) +* +* subject to +* +* sum{j in 1..n} a[j]z'[j] <= sum{j in 1..n} a[j] - b - eps (11) +* +* z'[j] in {0,1} for all j = 1,...,n (12) +* +* (The constant term zeta0 doesn't affect the solution, so it can be +* dropped.) */ + +static double simple_cover(int n, const double a[], double b, const + double x[], char z[]) +{ int j, *aa, bb, *cc; + double max_aj, min_aj, s, eps; + xassert(n >= 3); + /* allocate working arrays */ + aa = talloc(1+n, int); + cc = talloc(1+n, int); + /* compute max{j in 1..n} a[j] and min{j in 1..n} a[j] */ + max_aj = 0, min_aj = DBL_MAX; + for (j = 1; j <= n; j++) + { xassert(a[j] > 0); + if (max_aj < a[j]) + max_aj = a[j]; + if (min_aj > a[j]) + min_aj = a[j]; + } + /* scale and round constraint parameters to make them integral; + * note that we make the resulting inequality stronger than (11), + * so a[j]'s are rounded up while rhs is rounded down */ + s = 0; + for (j = 1; j <= n; j++) + { s += a[j]; + aa[j] = ceil(a[j] / max_aj * 1000); + } + bb = floor((s - b) / max_aj * 1000) - 1; + /* scale and round obj. coefficients to make them integral; + * again we make the objective function stronger than (10), so + * the coefficients are rounded down */ + for (j = 1; j <= n; j++) + { xassert(0 <= x[j] && x[j] <= 1); + cc[j] = floor((1 - x[j]) * 1000); + } + /* solve separation problem */ + if (solve_ks(n, aa, bb, cc, z) == INT_MIN) + { /* no cover exists */ + s = DBL_MAX; + goto skip; + } + /* determine z[j] = 1 - z'[j] */ + for (j = 1; j <= n; j++) + { xassert(z[j] == 0 || z[j] == 1); + z[j] ^= 1; + } + /* check condition (11) for original (non-scaled) parameters */ + s = 0; + for (j = 1; j <= n; j++) + { if (z[j]) + s += a[j]; + } + eps = 0.01 * (min_aj >= 1 ? min_aj : 1); + if (!(s >= b + eps)) + { /* no cover found within a precision req'd */ + s = DBL_MAX; + goto skip; + } + /* compute corresponding zeta (4) for cover found */ + s = 0; + for (j = 1; j <= n; j++) + { if (z[j]) + s += 1 - x[j]; + } +skip: /* free working arrays */ + tfree(aa); + tfree(cc); + return s; +} + +/**********************************************************************/ + +void glp_cov_gen1(glp_prob *P, glp_cov *cov, glp_prob *pool) +{ /* generate locally valid simple cover cuts */ + int i, k, len, new_len, *ind; + double *val, rhs, *x, zeta; + char *z; + xassert(P->n == cov->n && P->n == cov->set->n); + xassert(glp_get_status(P) == GLP_OPT); + /* allocate working arrays */ + ind = talloc(1+P->n, int); + val = talloc(1+P->n, double); + x = talloc(1+P->n, double); + z = talloc(1+P->n, char); + /* walk thru 0-1 knapsack inequalities */ + for (i = 1; i <= cov->set->m; i++) + { /* retrieve 0-1 knapsack inequality */ + len = glp_get_mat_row(cov->set, i, ind, val); + rhs = glp_get_row_ub(cov->set, i); + xassert(rhs != +DBL_MAX); + /* FIXME: skip, if slack is too large? */ + /* substitute and eliminate binary variables which have been + * fixed in the current subproblem (this makes the inequality + * only locally valid) */ + new_len = 0; + for (k = 1; k <= len; k++) + { if (glp_get_col_type(P, ind[k]) == GLP_FX) + rhs -= val[k] * glp_get_col_prim(P, ind[k]); + else + { new_len++; + ind[new_len] = ind[k]; + val[new_len] = val[k]; + } + } + len = new_len; + /* we need at least 3 binary variables in the inequality */ + if (len <= 2) + continue; + /* obtain values of binary variables from optimal solution to + * LP relaxation of current subproblem */ + for (k = 1; k <= len; k++) + { xassert(glp_get_col_kind(P, ind[k]) == GLP_BV); + x[k] = glp_get_col_prim(P, ind[k]); + if (x[k] < 0.00001) + x[k] = 0; + else if (x[k] > 0.99999) + x[k] = 1; + /* if val[k] < 0, perform substitution x[k] = 1 - x'[k] to + * make all coefficients positive */ + if (val[k] < 0) + { ind[k] = - ind[k]; /* x[k] is complemented */ + val[k] = - val[k]; + rhs += val[k]; + x[k] = 1 - x[k]; + } + } + /* find locally valid simple cover cut */ + zeta = simple_cover(len, val, rhs, x, z); + if (zeta > 0.95) + { /* no violation or insufficient violation; see (2) */ + continue; + } + /* construct cover inequality (2) for the cover found, which + * for original binary variables x[k] is equivalent to: + * sum{k in C'} x[k] + sum{k in C"} x'[k] <= |C| - 1 + * or + * sum{k in C'} x[k] + sum{k in C"} (1 - x[k]) <= |C| - 1 + * or + * sum{k in C'} x[k] - sum{k in C"} x[k] <= |C'| - 1 + * since |C| - |C"| = |C'| */ + new_len = 0; + rhs = -1; + for (k = 1; k <= len; k++) + { if (z[k]) + { new_len++; + if (ind[k] > 0) + { ind[new_len] = +ind[k]; + val[new_len] = +1; + rhs++; + } + else /* ind[k] < 0 */ + { ind[new_len] = -ind[k]; + val[new_len] = -1; + } + } + } + len = new_len; + /* add the cover inequality to the local cut pool */ + k = glp_add_rows(pool, 1); + glp_set_mat_row(pool, k, len, ind, val); + glp_set_row_bnds(pool, k, GLP_UP, rhs, rhs); + } + /* free working arrays */ + tfree(ind); + tfree(val); + tfree(x); + tfree(z); + return; +} + +/**********************************************************************/ + +void glp_cov_free(glp_cov *cov) +{ /* delete cover cut generator workspace */ + xassert(cov != NULL); + glp_delete_prob(cov->set); + tfree(cov); + return; +} + +/* eof */ diff --git a/src/glpk/glpios10.c b/src/glpk/intopt/fpump.c similarity index 96% rename from src/glpk/glpios10.c rename to src/glpk/intopt/fpump.c index 29f0b03..0bdd6d4 100644 --- a/src/glpk/glpios10.c +++ b/src/glpk/intopt/fpump.c @@ -1,10 +1,9 @@ -/* glpios10.c (feasibility pump heuristic) */ +/* fpump.c (feasibility pump heuristic) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013, 2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -22,12 +21,8 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "env.h" -#include "glpios.h" +#include "ios.h" #include "rng.h" /*********************************************************************** @@ -60,11 +55,7 @@ struct VAR /* sorting key */ }; -#ifndef _MSC_VER -static int fcmp(const void *x, const void *y) -#else -static int __cdecl fcmp(const void *x, const void *y) -#endif +static int CDECL fcmp(const void *x, const void *y) { /* comparison routine */ const struct VAR *vx = x, *vy = y; if (vx->d > vy->d) diff --git a/src/glpk/cglib/gmicut.c b/src/glpk/intopt/gmicut.c similarity index 100% rename from src/glpk/cglib/gmicut.c rename to src/glpk/intopt/gmicut.c diff --git a/src/glpk/cglib/gmigen.c b/src/glpk/intopt/gmigen.c similarity index 95% rename from src/glpk/cglib/gmigen.c rename to src/glpk/intopt/gmigen.c index 0051068..627682c 100644 --- a/src/glpk/cglib/gmigen.c +++ b/src/glpk/intopt/gmigen.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2002-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2002-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "env.h" #include "prob.h" @@ -64,11 +60,7 @@ struct var { int j; double f; }; -#ifndef _MSC_VER -static int fcmp(const void *p1, const void *p2) -#else -static int __cdecl fcmp(const void *p1, const void *p2) -#endif +static int CDECL fcmp(const void *p1, const void *p2) { const struct var *v1 = p1, *v2 = p2; if (v1->f > v2->f) return -1; if (v1->f < v2->f) return +1; diff --git a/src/glpk/cglib/mirgen.c b/src/glpk/intopt/mirgen.c similarity index 99% rename from src/glpk/cglib/mirgen.c rename to src/glpk/intopt/mirgen.c index 80b5286..45a0a55 100644 --- a/src/glpk/cglib/mirgen.c +++ b/src/glpk/intopt/mirgen.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2007-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2007-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #if 1 /* 29/II-2016 by Chris */ /*---------------------------------------------------------------------- Subject: Mir cut generation performance improvement @@ -802,11 +798,7 @@ static int cmir_ineq(const int n, const double a[], const double b, struct vset { int j; double v; }; -#ifndef _MSC_VER -static int cmir_cmp(const void *p1, const void *p2) -#else -static int __cdecl cmir_cmp(const void *p1, const void *p2) -#endif +static int CDECL cmir_cmp(const void *p1, const void *p2) { const struct vset *v1 = p1, *v2 = p2; if (v1->v < v2->v) return -1; if (v1->v > v2->v) return +1; diff --git a/src/glpk/cglib/spv.c b/src/glpk/intopt/spv.c similarity index 100% rename from src/glpk/cglib/spv.c rename to src/glpk/intopt/spv.c diff --git a/src/glpk/cglib/spv.h b/src/glpk/intopt/spv.h similarity index 100% rename from src/glpk/cglib/spv.h rename to src/glpk/intopt/spv.h diff --git a/src/glpk/bflib/fvs.c b/src/glpk/misc/fvs.c similarity index 100% rename from src/glpk/bflib/fvs.c rename to src/glpk/misc/fvs.c diff --git a/src/glpk/bflib/fvs.h b/src/glpk/misc/fvs.h similarity index 97% rename from src/glpk/bflib/fvs.h rename to src/glpk/misc/fvs.h index 5da0a19..abfed8c 100644 --- a/src/glpk/bflib/fvs.h +++ b/src/glpk/misc/fvs.h @@ -36,7 +36,7 @@ struct FVS /* ind[0] is not used; * ind[k] = j, 1 <= k <= nnz, means that vec[j] != 0 * non-zero indices in the array ind are stored in arbitrary - * order; if vec[j] = 0, its index j should NOT be presented in + * order; if vec[j] = 0, its index j SHOULD NOT be presented in * the array ind */ double *vec; /* double vec[1+n]; */ /* vec[0] is not used; diff --git a/src/glpk/misc/ks.c b/src/glpk/misc/ks.c new file mode 100644 index 0000000..0720cc9 --- /dev/null +++ b/src/glpk/misc/ks.c @@ -0,0 +1,466 @@ +/* ks.c (0-1 knapsack problem) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2017-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "env.h" +#include "ks.h" +#include "mt1.h" + +/*********************************************************************** +* 0-1 knapsack problem has the following formulation: +* +* maximize z = sum{j in 1..n} c[j]x[j] (1) +* +* s.t. sum{j in 1..n} a[j]x[j] <= b (2) +* +* x[j] in {0, 1} for all j in 1..n (3) +* +* In general case it is assumed that the instance is non-normalized, +* i.e. parameters a, b, and c may have any sign. +***********************************************************************/ + +/*********************************************************************** +* ks_enum - solve 0-1 knapsack problem by complete enumeration +* +* This routine finds optimal solution to 0-1 knapsack problem (1)-(3) +* by complete enumeration. It is intended mainly for testing purposes. +* +* The instance to be solved is specified by parameters n, a, b, and c. +* Note that these parameters can have any sign, i.e. normalization is +* not needed. +* +* On exit the routine stores the optimal point found in locations +* x[1], ..., x[n] and returns the optimal objective value. However, if +* the instance is infeasible, the routine returns INT_MIN. +* +* Since the complete enumeration is inefficient, this routine can be +* used only for small instances (n <= 20-30). */ + +#define N_MAX 40 + +int ks_enum(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]) +{ int j, s, z, z_best; + char x_best[1+N_MAX]; + xassert(0 <= n && n <= N_MAX); + /* initialization */ + memset(&x[1], 0, n * sizeof(char)); + z_best = INT_MIN; +loop: /* compute constraint and objective at current x */ + s = z = 0; + for (j = 1; j <= n; j++) + { if (x[j]) + s += a[j], z += c[j]; + } + /* check constraint violation */ + if (s > b) + goto next; + /* check objective function */ + if (z_best < z) + { /* better solution has been found */ + memcpy(&x_best[1], &x[1], n * sizeof(char)); + z_best = z; + } +next: /* generate next x */ + for (j = 1; j <= n; j++) + { if (!x[j]) + { x[j] = 1; + goto loop; + } + x[j] = 0; + } + /* report best (optimal) solution */ + memcpy(&x[1], &x_best[1], n * sizeof(char)); + return z_best; +} + +/*********************************************************************** +* reduce - prepare reduced instance of 0-1 knapsack +* +* Given original instance of 0-1 knapsack (1)-(3) specified by the +* parameters n, a, b, and c this routine transforms it to equivalent +* reduced instance in the same format. The reduced instance is +* normalized, i.e. the following additional conditions are met: +* +* n >= 2 (4) +* +* 1 <= a[j] <= b for all j in 1..n (5) +* +* sum{j in 1..n} a[j] >= b+1 (6) +* +* c[j] >= 1 for all j in 1..n (7) +* +* The routine creates the structure ks and stores there parameters n, +* a, b, and c of the reduced instance as well as template of solution +* to original instance. +* +* Normally the routine returns a pointer to the structure ks created. +* However, if the original instance is infeasible, the routine returns +* a null pointer. */ + +struct ks +{ int orig_n; + /* original problem dimension */ + int n; + /* reduced problem dimension */ + int *a; /* int a[1+orig_n]; */ + /* a{j in 1..n} are constraint coefficients (2) */ + int b; + /* b is constraint right-hand side (2) */ + int *c; /* int c[1+orig_n]; */ + /* c{j in 1..n} are objective coefficients (1) */ + int c0; + /* c0 is objective constant term */ + char *x; /* char x[1+orig_n]; */ + /* x{j in 1..orig_n} is solution template to original instance: + * x[j] = 0 x[j] is fixed at 0 + * x[j] = 1 x[j] is fixed at 1 + * x[j] = 0x10 x[j] = x[j'] + * x[j] = 0x11 x[j] = 1 - x[j'] + * where x[j'] is corresponding solution to reduced instance */ +}; + +static void free_ks(struct ks *ks); + +static struct ks *reduce(const int n, const int a[/*1+n*/], int b, + const int c[/*1+n*/]) +{ struct ks *ks; + int j, s; + xassert(n >= 0); + /* initially reduced instance is the same as original one */ + ks = talloc(1, struct ks); + ks->orig_n = n; + ks->n = 0; + ks->a = talloc(1+n, int); + memcpy(&ks->a[1], &a[1], n * sizeof(int)); + ks->b = b; + ks->c = talloc(1+n, int); + memcpy(&ks->c[1], &c[1], n * sizeof(int)); + ks->c0 = 0; + ks->x = talloc(1+n, char); + /* make all a[j] non-negative */ + for (j = 1; j <= n; j++) + { if (a[j] >= 0) + { /* keep original x[j] */ + ks->x[j] = 0x10; + } + else /* a[j] < 0 */ + { /* substitute x[j] = 1 - x'[j] */ + ks->x[j] = 0x11; + /* ... + a[j]x[j] + ... <= b + * ... + a[j](1 - x'[j]) + ... <= b + * ... - a[j]x'[j] + ... <= b - a[j] */ + ks->a[j] = - ks->a[j]; + ks->b += ks->a[j]; + /* z = ... + c[j]x[j] + ... + c0 = + * = ... + c[j](1 - x'[j]) + ... + c0 = + * = ... - c[j]x'[j] + ... + (c0 + c[j]) */ + ks->c0 += ks->c[j]; + ks->c[j] = - ks->c[j]; + } + } + /* now a[j] >= 0 for all j in 1..n */ + if (ks->b < 0) + { /* instance is infeasible */ + free_ks(ks); + return NULL; + } + /* build reduced instance */ + for (j = 1; j <= n; j++) + { if (ks->a[j] == 0) + { if (ks->c[j] <= 0) + { /* fix x[j] at 0 */ + ks->x[j] ^= 0x10; + } + else + { /* fix x[j] at 1 */ + ks->x[j] ^= 0x11; + ks->c0 += ks->c[j]; + } + } + else if (ks->a[j] > ks->b || ks->c[j] <= 0) + { /* fix x[j] at 0 */ + ks->x[j] ^= 0x10; + } + else + { /* include x[j] in reduced instance */ + ks->n++; + ks->a[ks->n] = ks->a[j]; + ks->c[ks->n] = ks->c[j]; + } + } + /* now conditions (5) and (7) are met */ + /* check condition (6) */ + s = 0; + for (j = 1; j <= ks->n; j++) + { xassert(1 <= ks->a[j] && ks->a[j] <= ks->b); + xassert(ks->c[j] >= 1); + s += ks->a[j]; + } + if (s <= ks->b) + { /* sum{j in 1..n} a[j] <= b */ + /* fix all remaining x[j] at 1 to obtain trivial solution */ + for (j = 1; j <= n; j++) + { if (ks->x[j] & 0x10) + ks->x[j] ^= 0x11; + } + for (j = 1; j <= ks->n; j++) + ks->c0 += ks->c[j]; + /* reduced instance is empty */ + ks->n = 0; + } + /* here n = 0 or n >= 2 due to condition (6) */ + xassert(ks->n == 0 || ks->n >= 2); + return ks; +} + +/*********************************************************************** +* restore - restore solution to original 0-1 knapsack instance +* +* Given optimal solution x{j in 1..ks->n} to the reduced 0-1 knapsack +* instance (previously prepared by the routine reduce) this routine +* constructs optimal solution to the original instance and stores it +* in the array ks->x{j in 1..ks->orig_n}. +* +* On exit the routine returns optimal objective value for the original +* instance. +* +* NOTE: This operation should be performed only once. */ + +static int restore(struct ks *ks, char x[]) +{ int j, k, z; + z = ks->c0; + for (j = 1, k = 0; j <= ks->orig_n; j++) + { if (ks->x[j] & 0x10) + { k++; + xassert(k <= ks->n); + xassert(x[k] == 0 || x[k] == 1); + if (ks->x[j] & 1) + ks->x[j] = 1 - x[k]; + else + ks->x[j] = x[k]; + if (x[k]) + z += ks->c[k]; + } + } + xassert(k == ks->n); + return z; +} + +/*********************************************************************** +* free_ks - deallocate structure ks +* +* This routine frees memory previously allocated to the structure ks +* and all its components. */ + +static void free_ks(struct ks *ks) +{ xassert(ks != NULL); + tfree(ks->a); + tfree(ks->c); + tfree(ks->x); + tfree(ks); +} + +/*********************************************************************** +* ks_mt1 - solve 0-1 knapsack problem with Martello & Toth algorithm +* +* This routine finds optimal solution to 0-1 knapsack problem (1)-(3) +* with Martello & Toth algorithm MT1. +* +* The instance to be solved is specified by parameters n, a, b, and c. +* Note that these parameters can have any sign, i.e. normalization is +* not needed. +* +* On exit the routine stores the optimal point found in locations +* x[1], ..., x[n] and returns the optimal objective value. However, if +* the instance is infeasible, the routine returns INT_MIN. +* +* REFERENCES +* +* S.Martello, P.Toth. Knapsack Problems: Algorithms and Computer Imp- +* lementations. John Wiley & Sons, 1990. */ + +struct mt +{ int j; + float r; /* r[j] = c[j] / a[j] */ +}; + +static int CDECL fcmp(const void *p1, const void *p2) +{ if (((struct mt *)p1)->r > ((struct mt *)p2)->r) + return -1; + else if (((struct mt *)p1)->r < ((struct mt *)p2)->r) + return +1; + else + return 0; +} + +static int mt1a(int n, const int a[], int b, const int c[], char x[]) +{ /* interface routine to MT1 */ + struct mt *mt; + int j, z, *p, *w, *x1, *xx, *min, *psign, *wsign, *zsign; + xassert(n >= 2); + /* allocate working arrays */ + mt = talloc(1+n, struct mt); + p = talloc(1+n+1, int); + w = talloc(1+n+1, int); + x1 = talloc(1+n+1, int); + xx = talloc(1+n+1, int); + min = talloc(1+n+1, int); + psign = talloc(1+n+1, int); + wsign = talloc(1+n+1, int); + zsign = talloc(1+n+1, int); + /* reorder items to provide c[j] / a[j] >= a[j+1] / a[j+1] */ + for (j = 1; j <= n; j++) + { mt[j].j = j; + mt[j].r = (float)c[j] / (float)a[j]; + } + qsort(&mt[1], n, sizeof(struct mt), fcmp); + /* load instance parameters */ + for (j = 1; j <= n; j++) + { p[j] = c[mt[j].j]; + w[j] = a[mt[j].j]; + } + /* find optimal solution */ + z = mt1(n, p, w, b, x1, 1, xx, min, psign, wsign, zsign); + xassert(z >= 0); + /* store optimal point found */ + for (j = 1; j <= n; j++) + { xassert(x1[j] == 0 || x1[j] == 1); + x[mt[j].j] = x1[j]; + } + /* free working arrays */ + tfree(mt); + tfree(p); + tfree(w); + tfree(x1); + tfree(xx); + tfree(min); + tfree(psign); + tfree(wsign); + tfree(zsign); + return z; +} + +int ks_mt1(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]) +{ struct ks *ks; + int j, s1, s2, z; + xassert(n >= 0); + /* prepare reduced instance */ + ks = reduce(n, a, b, c); + if (ks == NULL) + { /* original instance is infeasible */ + return INT_MIN; + } + /* find optimal solution to reduced instance */ + if (ks->n > 0) + mt1a(ks->n, ks->a, ks->b, ks->c, x); + /* restore solution to original instance */ + z = restore(ks, x); + memcpy(&x[1], &ks->x[1], n * sizeof(char)); + free_ks(ks); + /* check solution found */ + s1 = s2 = 0; + for (j = 1; j <= n; j++) + { xassert(x[j] == 0 || x[j] == 1); + if (x[j]) + s1 += a[j], s2 += c[j]; + } + xassert(s1 <= b); + xassert(s2 == z); + return z; +} + +/*********************************************************************** +* ks_greedy - solve 0-1 knapsack problem with greedy heuristic +* +* This routine finds (sub)optimal solution to 0-1 knapsack problem +* (1)-(3) with greedy heuristic. +* +* The instance to be solved is specified by parameters n, a, b, and c. +* Note that these parameters can have any sign, i.e. normalization is +* not needed. +* +* On exit the routine stores the optimal point found in locations +* x[1], ..., x[n] and returns the optimal objective value. However, if +* the instance is infeasible, the routine returns INT_MIN. */ + +static int greedy(int n, const int a[], int b, const int c[], char x[]) +{ /* core routine for normalized 0-1 knapsack instance */ + struct mt *mt; + int j, s, z; + xassert(n >= 2); + /* reorder items to provide c[j] / a[j] >= a[j+1] / a[j+1] */ + mt = talloc(1+n, struct mt); + for (j = 1; j <= n; j++) + { mt[j].j = j; + mt[j].r = (float)c[j] / (float)a[j]; + } + qsort(&mt[1], n, sizeof(struct mt), fcmp); + /* take items starting from most valuable ones until the knapsack + * is full */ + s = z = 0; + for (j = 1; j <= n; j++) + { if (s + a[mt[j].j] > b) + break; + x[mt[j].j] = 1; + s += a[mt[j].j]; + z += c[mt[j].j]; + } + /* don't take remaining items */ + for (j = j; j <= n; j++) + x[mt[j].j] = 0; + tfree(mt); + return z; +} + +int ks_greedy(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]) +{ struct ks *ks; + int j, s1, s2, z; + xassert(n >= 0); + /* prepare reduced instance */ + ks = reduce(n, a, b, c); + if (ks == NULL) + { /* original instance is infeasible */ + return INT_MIN; + } + /* find suboptimal solution to reduced instance */ + if (ks->n > 0) + greedy(ks->n, ks->a, ks->b, ks->c, x); + /* restore solution to original instance */ + z = restore(ks, x); + memcpy(&x[1], &ks->x[1], n * sizeof(char)); + free_ks(ks); + /* check solution found */ + s1 = s2 = 0; + for (j = 1; j <= n; j++) + { xassert(x[j] == 0 || x[j] == 1); + if (x[j]) + s1 += a[j], s2 += c[j]; + } + xassert(s1 <= b); + xassert(s2 == z); + return z; +} + +/* eof */ diff --git a/src/glpk/misc/ks.h b/src/glpk/misc/ks.h new file mode 100644 index 0000000..d607dc4 --- /dev/null +++ b/src/glpk/misc/ks.h @@ -0,0 +1,44 @@ +/* ks.h (0-1 knapsack problem) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2017-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#ifndef KS_H +#define KS_H + +#define ks_enum _glp_ks_enum +int ks_enum(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]); +/* solve 0-1 knapsack problem by complete enumeration */ + +#define ks_mt1 _glp_ks_mt1 +int ks_mt1(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]); +/* solve 0-1 knapsack problem with Martello & Toth algorithm */ + +#define ks_greedy _glp_ks_greedy +int ks_greedy(int n, const int a[/*1+n*/], int b, const int c[/*1+n*/], + char x[/*1+n*/]); +/* solve 0-1 knapsack problem with greedy heuristic */ + +#endif + +/* eof */ diff --git a/src/glpk/misc/mt1.c b/src/glpk/misc/mt1.c new file mode 100644 index 0000000..63a0f80 --- /dev/null +++ b/src/glpk/misc/mt1.c @@ -0,0 +1,1110 @@ +/* mt1.c (0-1 knapsack problem; Martello & Toth algorithm) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* THIS CODE IS THE RESULT OF TRANSLATION OF THE FORTRAN SUBROUTINES +* MT1 FROM THE BOOK: +* +* SILVANO MARTELLO, PAOLO TOTH. KNAPSACK PROBLEMS: ALGORITHMS AND +* COMPUTER IMPLEMENTATIONS. JOHN WILEY & SONS, 1990. +* +* THE TRANSLATION HAS BEEN DONE WITH THE PERMISSION OF THE AUTHORS OF +* THE ORIGINAL FORTRAN SUBROUTINES: SILVANO MARTELLO AND PAOLO TOTH. +* +* The translation was made by Andrew Makhorin . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#line 1 "" +/* -- translated by f2c (version 20100827). + You must link the resulting object file with libf2c: + on Microsoft Windows system, link with libf2c.lib; + on Linux or Unix systems, link with .../path/to/libf2c.a -lm + or, if you install libf2c.a in a standard place, with -lf2c -lm + -- in that order, at the end of the command line, as in + cc *.o -lf2c -lm + Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., + + http://www.netlib.org/f2c/libf2c.zip +*/ + +#if 0 /* by mao */ +#include "f2c.h" +#else +#include "env.h" +#include "mt1.h" + +typedef int integer; +typedef float real; +#endif + +#line 1 "" +/*< SUBROUTINE MT1(N,P,W,C,Z,X,JDIM,JCK,XX,MIN,PSIGN,WSIGN,ZSIGN) >*/ +#if 1 /* by mao */ +static int chmt1_(int *, int *, int *, int *, int *, int *); + +static +#endif +/* Subroutine */ int mt1_(integer *n, integer *p, integer *w, integer *c__, + integer *z__, integer *x, integer *jdim, integer *jck, integer *xx, + integer *min__, integer *psign, integer *wsign, integer *zsign) +{ + /* System generated locals */ + integer i__1; + + /* Local variables */ + static real a, b; + static integer j, r__, t, j1, n1, ch, ii, jj, kk, in, ll, ip, nn, iu, ii1, + chs, lim, lim1, diff, lold, mink; + extern /* Subroutine */ int chmt1_(integer *, integer *, integer *, + integer *, integer *, integer *); + static integer profit; + + +/* THIS SUBROUTINE SOLVES THE 0-1 SINGLE KNAPSACK PROBLEM */ + +/* MAXIMIZE Z = P(1)*X(1) + ... + P(N)*X(N) */ + +/* SUBJECT TO: W(1)*X(1) + ... + W(N)*X(N) .LE. C , */ +/* X(J) = 0 OR 1 FOR J=1,...,N. */ + +/* THE PROGRAM IS INCLUDED IN THE VOLUME */ +/* S. MARTELLO, P. TOTH, "KNAPSACK PROBLEMS: ALGORITHMS */ +/* AND COMPUTER IMPLEMENTATIONS", JOHN WILEY, 1990 */ +/* AND IMPLEMENTS THE BRANCH-AND-BOUND ALGORITHM DESCRIBED IN */ +/* SECTION 2.5.2 . */ +/* THE PROGRAM DERIVES FROM AN EARLIER CODE PRESENTED IN */ +/* S. MARTELLO, P. TOTH, "ALGORITHM FOR THE SOLUTION OF THE 0-1 SINGLE */ +/* KNAPSACK PROBLEM", COMPUTING, 1978. */ + +/* THE INPUT PROBLEM MUST SATISFY THE CONDITIONS */ + +/* 1) 2 .LE. N .LE. JDIM - 1 ; */ +/* 2) P(J), W(J), C POSITIVE INTEGERS; */ +/* 3) MAX (W(J)) .LE. C ; */ +/* 4) W(1) + ... + W(N) .GT. C ; */ +/* 5) P(J)/W(J) .GE. P(J+1)/W(J+1) FOR J=1,...,N-1. */ + +/* MT1 CALLS 1 PROCEDURE: CHMT1. */ + +/* THE PROGRAM IS COMPLETELY SELF-CONTAINED AND COMMUNICATION TO IT IS */ +/* ACHIEVED SOLELY THROUGH THE PARAMETER LIST OF MT1. */ +/* NO MACHINE-DEPENDENT CONSTANT IS USED. */ +/* THE PROGRAM IS WRITTEN IN 1967 AMERICAN NATIONAL STANDARD FORTRAN */ +/* AND IS ACCEPTED BY THE PFORT VERIFIER (PFORT IS THE PORTABLE */ +/* SUBSET OF ANSI DEFINED BY THE ASSOCIATION FOR COMPUTING MACHINERY). */ +/* THE PROGRAM HAS BEEN TESTED ON A DIGITAL VAX 11/780 AND AN H.P. */ +/* 9000/840. */ + +/* MT1 NEEDS 8 ARRAYS ( P , W , X , XX , MIN , PSIGN , WSIGN */ +/* AND ZSIGN ) OF LENGTH AT LEAST N + 1 . */ + +/* MEANING OF THE INPUT PARAMETERS: */ +/* N = NUMBER OF ITEMS; */ +/* P(J) = PROFIT OF ITEM J (J=1,...,N); */ +/* W(J) = WEIGHT OF ITEM J (J=1,...,N); */ +/* C = CAPACITY OF THE KNAPSACK; */ +/* JDIM = DIMENSION OF THE 8 ARRAYS; */ +/* JCK = 1 IF CHECK ON THE INPUT DATA IS DESIRED, */ +/* = 0 OTHERWISE. */ + +/* MEANING OF THE OUTPUT PARAMETERS: */ +/* Z = VALUE OF THE OPTIMAL SOLUTION IF Z .GT. 0 , */ +/* = ERROR IN THE INPUT DATA (WHEN JCK=1) IF Z .LT. 0 : CONDI- */ +/* TION - Z IS VIOLATED; */ +/* X(J) = 1 IF ITEM J IS IN THE OPTIMAL SOLUTION, */ +/* = 0 OTHERWISE. */ + +/* ARRAYS XX, MIN, PSIGN, WSIGN AND ZSIGN ARE DUMMY. */ + +/* ALL THE PARAMETERS ARE INTEGER. ON RETURN OF MT1 ALL THE INPUT */ +/* PARAMETERS ARE UNCHANGED. */ + +/*< INTEGER P(JDIM),W(JDIM),X(JDIM),C,Z >*/ +/*< INTEGER XX(JDIM),MIN(JDIM),PSIGN(JDIM),WSIGN(JDIM),ZSIGN(JDIM) >*/ +/*< INTEGER CH,CHS,DIFF,PROFIT,R,T >*/ +/*< Z = 0 >*/ +#line 65 "" + /* Parameter adjustments */ +#line 65 "" + --zsign; +#line 65 "" + --wsign; +#line 65 "" + --psign; +#line 65 "" + --min__; +#line 65 "" + --xx; +#line 65 "" + --x; +#line 65 "" + --w; +#line 65 "" + --p; +#line 65 "" + +#line 65 "" + /* Function Body */ +#line 65 "" + *z__ = 0; +/*< IF ( JCK .EQ. 1 ) CALL CHMT1(N,P,W,C,Z,JDIM) >*/ +#line 66 "" + if (*jck == 1) { +#line 66 "" + chmt1_(n, &p[1], &w[1], c__, z__, jdim); +#line 66 "" + } +/*< IF ( Z .LT. 0 ) RETURN >*/ +#line 67 "" + if (*z__ < 0) { +#line 67 "" + return 0; +#line 67 "" + } +/* INITIALIZE. */ +/*< CH = C >*/ +#line 69 "" + ch = *c__; +/*< IP = 0 >*/ +#line 70 "" + ip = 0; +/*< CHS = CH >*/ +#line 71 "" + chs = ch; +/*< DO 10 LL=1,N >*/ +#line 72 "" + i__1 = *n; +#line 72 "" + for (ll = 1; ll <= i__1; ++ll) { +/*< IF ( W(LL) .GT. CHS ) GO TO 20 >*/ +#line 73 "" + if (w[ll] > chs) { +#line 73 "" + goto L20; +#line 73 "" + } +/*< IP = IP + P(LL) >*/ +#line 74 "" + ip += p[ll]; +/*< CHS = CHS - W(LL) >*/ +#line 75 "" + chs -= w[ll]; +/*< 10 CONTINUE >*/ +#line 76 "" +/* L10: */ +#line 76 "" + } +/*< 20 LL = LL - 1 >*/ +#line 77 "" +L20: +#line 77 "" + --ll; +/*< IF ( CHS .EQ. 0 ) GO TO 50 >*/ +#line 78 "" + if (chs == 0) { +#line 78 "" + goto L50; +#line 78 "" + } +/*< P(N+1) = 0 >*/ +#line 79 "" + p[*n + 1] = 0; +/*< W(N+1) = CH + 1 >*/ +#line 80 "" + w[*n + 1] = ch + 1; +/*< LIM = IP + CHS*P(LL+2)/W(LL+2) >*/ +#line 81 "" + lim = ip + chs * p[ll + 2] / w[ll + 2]; +/*< A = W(LL+1) - CHS >*/ +#line 82 "" + a = (real) (w[ll + 1] - chs); +/*< B = IP + P(LL+1) >*/ +#line 83 "" + b = (real) (ip + p[ll + 1]); +/*< LIM1 = B - A*FLOAT(P(LL))/FLOAT(W(LL)) >*/ +#line 84 "" + lim1 = b - a * (real) p[ll] / (real) w[ll]; +/*< IF ( LIM1 .GT. LIM ) LIM = LIM1 >*/ +#line 85 "" + if (lim1 > lim) { +#line 85 "" + lim = lim1; +#line 85 "" + } +/*< MINK = CH + 1 >*/ +#line 86 "" + mink = ch + 1; +/*< MIN(N) = MINK >*/ +#line 87 "" + min__[*n] = mink; +/*< DO 30 J=2,N >*/ +#line 88 "" + i__1 = *n; +#line 88 "" + for (j = 2; j <= i__1; ++j) { +/*< KK = N + 2 - J >*/ +#line 89 "" + kk = *n + 2 - j; +/*< IF ( W(KK) .LT. MINK ) MINK = W(KK) >*/ +#line 90 "" + if (w[kk] < mink) { +#line 90 "" + mink = w[kk]; +#line 90 "" + } +/*< MIN(KK-1) = MINK >*/ +#line 91 "" + min__[kk - 1] = mink; +/*< 30 CONTINUE >*/ +#line 92 "" +/* L30: */ +#line 92 "" + } +/*< DO 40 J=1,N >*/ +#line 93 "" + i__1 = *n; +#line 93 "" + for (j = 1; j <= i__1; ++j) { +/*< XX(J) = 0 >*/ +#line 94 "" + xx[j] = 0; +/*< 40 CONTINUE >*/ +#line 95 "" +/* L40: */ +#line 95 "" + } +/*< Z = 0 >*/ +#line 96 "" + *z__ = 0; +/*< PROFIT = 0 >*/ +#line 97 "" + profit = 0; +/*< LOLD = N >*/ +#line 98 "" + lold = *n; +/*< II = 1 >*/ +#line 99 "" + ii = 1; +/*< GO TO 170 >*/ +#line 100 "" + goto L170; +/*< 50 Z = IP >*/ +#line 101 "" +L50: +#line 101 "" + *z__ = ip; +/*< DO 60 J=1,LL >*/ +#line 102 "" + i__1 = ll; +#line 102 "" + for (j = 1; j <= i__1; ++j) { +/*< X(J) = 1 >*/ +#line 103 "" + x[j] = 1; +/*< 60 CONTINUE >*/ +#line 104 "" +/* L60: */ +#line 104 "" + } +/*< NN = LL + 1 >*/ +#line 105 "" + nn = ll + 1; +/*< DO 70 J=NN,N >*/ +#line 106 "" + i__1 = *n; +#line 106 "" + for (j = nn; j <= i__1; ++j) { +/*< X(J) = 0 >*/ +#line 107 "" + x[j] = 0; +/*< 70 CONTINUE >*/ +#line 108 "" +/* L70: */ +#line 108 "" + } +/*< RETURN >*/ +#line 109 "" + return 0; +/* TRY TO INSERT THE II-TH ITEM INTO THE CURRENT SOLUTION. */ +/*< 80 IF ( W(II) .LE. CH ) GO TO 90 >*/ +#line 111 "" +L80: +#line 111 "" + if (w[ii] <= ch) { +#line 111 "" + goto L90; +#line 111 "" + } +/*< II1 = II + 1 >*/ +#line 112 "" + ii1 = ii + 1; +/*< IF ( Z .GE. CH*P(II1)/W(II1) + PROFIT ) GO TO 280 >*/ +#line 113 "" + if (*z__ >= ch * p[ii1] / w[ii1] + profit) { +#line 113 "" + goto L280; +#line 113 "" + } +/*< II = II1 >*/ +#line 114 "" + ii = ii1; +/*< GO TO 80 >*/ +#line 115 "" + goto L80; +/* BUILD A NEW CURRENT SOLUTION. */ +/*< 90 IP = PSIGN(II) >*/ +#line 117 "" +L90: +#line 117 "" + ip = psign[ii]; +/*< CHS = CH - WSIGN(II) >*/ +#line 118 "" + chs = ch - wsign[ii]; +/*< IN = ZSIGN(II) >*/ +#line 119 "" + in = zsign[ii]; +/*< DO 100 LL=IN,N >*/ +#line 120 "" + i__1 = *n; +#line 120 "" + for (ll = in; ll <= i__1; ++ll) { +/*< IF ( W(LL) .GT. CHS ) GO TO 160 >*/ +#line 121 "" + if (w[ll] > chs) { +#line 121 "" + goto L160; +#line 121 "" + } +/*< IP = IP + P(LL) >*/ +#line 122 "" + ip += p[ll]; +/*< CHS = CHS - W(LL) >*/ +#line 123 "" + chs -= w[ll]; +/*< 100 CONTINUE >*/ +#line 124 "" +/* L100: */ +#line 124 "" + } +/*< LL = N >*/ +#line 125 "" + ll = *n; +/*< 110 IF ( Z .GE. IP + PROFIT ) GO TO 280 >*/ +#line 126 "" +L110: +#line 126 "" + if (*z__ >= ip + profit) { +#line 126 "" + goto L280; +#line 126 "" + } +/*< Z = IP + PROFIT >*/ +#line 127 "" + *z__ = ip + profit; +/*< NN = II - 1 >*/ +#line 128 "" + nn = ii - 1; +/*< DO 120 J=1,NN >*/ +#line 129 "" + i__1 = nn; +#line 129 "" + for (j = 1; j <= i__1; ++j) { +/*< X(J) = XX(J) >*/ +#line 130 "" + x[j] = xx[j]; +/*< 120 CONTINUE >*/ +#line 131 "" +/* L120: */ +#line 131 "" + } +/*< DO 130 J=II,LL >*/ +#line 132 "" + i__1 = ll; +#line 132 "" + for (j = ii; j <= i__1; ++j) { +/*< X(J) = 1 >*/ +#line 133 "" + x[j] = 1; +/*< 130 CONTINUE >*/ +#line 134 "" +/* L130: */ +#line 134 "" + } +/*< IF ( LL .EQ. N ) GO TO 150 >*/ +#line 135 "" + if (ll == *n) { +#line 135 "" + goto L150; +#line 135 "" + } +/*< NN = LL + 1 >*/ +#line 136 "" + nn = ll + 1; +/*< DO 140 J=NN,N >*/ +#line 137 "" + i__1 = *n; +#line 137 "" + for (j = nn; j <= i__1; ++j) { +/*< X(J) = 0 >*/ +#line 138 "" + x[j] = 0; +/*< 140 CONTINUE >*/ +#line 139 "" +/* L140: */ +#line 139 "" + } +/*< 150 IF ( Z .NE. LIM ) GO TO 280 >*/ +#line 140 "" +L150: +#line 140 "" + if (*z__ != lim) { +#line 140 "" + goto L280; +#line 140 "" + } +/*< RETURN >*/ +#line 141 "" + return 0; +/*< 160 IU = CHS*P(LL)/W(LL) >*/ +#line 142 "" +L160: +#line 142 "" + iu = chs * p[ll] / w[ll]; +/*< LL = LL - 1 >*/ +#line 143 "" + --ll; +/*< IF ( IU .EQ. 0 ) GO TO 110 >*/ +#line 144 "" + if (iu == 0) { +#line 144 "" + goto L110; +#line 144 "" + } +/*< IF ( Z .GE. PROFIT + IP + IU ) GO TO 280 >*/ +#line 145 "" + if (*z__ >= profit + ip + iu) { +#line 145 "" + goto L280; +#line 145 "" + } +/* SAVE THE CURRENT SOLUTION. */ +/*< 170 WSIGN(II) = CH - CHS >*/ +#line 147 "" +L170: +#line 147 "" + wsign[ii] = ch - chs; +/*< PSIGN(II) = IP >*/ +#line 148 "" + psign[ii] = ip; +/*< ZSIGN(II) = LL + 1 >*/ +#line 149 "" + zsign[ii] = ll + 1; +/*< XX(II) = 1 >*/ +#line 150 "" + xx[ii] = 1; +/*< NN = LL - 1 >*/ +#line 151 "" + nn = ll - 1; +/*< IF ( NN .LT. II) GO TO 190 >*/ +#line 152 "" + if (nn < ii) { +#line 152 "" + goto L190; +#line 152 "" + } +/*< DO 180 J=II,NN >*/ +#line 153 "" + i__1 = nn; +#line 153 "" + for (j = ii; j <= i__1; ++j) { +/*< WSIGN(J+1) = WSIGN(J) - W(J) >*/ +#line 154 "" + wsign[j + 1] = wsign[j] - w[j]; +/*< PSIGN(J+1) = PSIGN(J) - P(J) >*/ +#line 155 "" + psign[j + 1] = psign[j] - p[j]; +/*< ZSIGN(J+1) = LL + 1 >*/ +#line 156 "" + zsign[j + 1] = ll + 1; +/*< XX(J+1) = 1 >*/ +#line 157 "" + xx[j + 1] = 1; +/*< 180 CONTINUE >*/ +#line 158 "" +/* L180: */ +#line 158 "" + } +/*< 190 J1 = LL + 1 >*/ +#line 159 "" +L190: +#line 159 "" + j1 = ll + 1; +/*< DO 200 J=J1,LOLD >*/ +#line 160 "" + i__1 = lold; +#line 160 "" + for (j = j1; j <= i__1; ++j) { +/*< WSIGN(J) = 0 >*/ +#line 161 "" + wsign[j] = 0; +/*< PSIGN(J) = 0 >*/ +#line 162 "" + psign[j] = 0; +/*< ZSIGN(J) = J >*/ +#line 163 "" + zsign[j] = j; +/*< 200 CONTINUE >*/ +#line 164 "" +/* L200: */ +#line 164 "" + } +/*< LOLD = LL >*/ +#line 165 "" + lold = ll; +/*< CH = CHS >*/ +#line 166 "" + ch = chs; +/*< PROFIT = PROFIT + IP >*/ +#line 167 "" + profit += ip; +/*< IF ( LL - (N - 2) ) 240, 220, 210 >*/ +#line 168 "" + if ((i__1 = ll - (*n - 2)) < 0) { +#line 168 "" + goto L240; +#line 168 "" + } else if (i__1 == 0) { +#line 168 "" + goto L220; +#line 168 "" + } else { +#line 168 "" + goto L210; +#line 168 "" + } +/*< 210 II = N >*/ +#line 169 "" +L210: +#line 169 "" + ii = *n; +/*< GO TO 250 >*/ +#line 170 "" + goto L250; +/*< 220 IF ( CH .LT. W(N) ) GO TO 230 >*/ +#line 171 "" +L220: +#line 171 "" + if (ch < w[*n]) { +#line 171 "" + goto L230; +#line 171 "" + } +/*< CH = CH - W(N) >*/ +#line 172 "" + ch -= w[*n]; +/*< PROFIT = PROFIT + P(N) >*/ +#line 173 "" + profit += p[*n]; +/*< XX(N) = 1 >*/ +#line 174 "" + xx[*n] = 1; +/*< 230 II = N - 1 >*/ +#line 175 "" +L230: +#line 175 "" + ii = *n - 1; +/*< GO TO 250 >*/ +#line 176 "" + goto L250; +/*< 240 II = LL + 2 >*/ +#line 177 "" +L240: +#line 177 "" + ii = ll + 2; +/*< IF ( CH .GE. MIN(II-1) ) GO TO 80 >*/ +#line 178 "" + if (ch >= min__[ii - 1]) { +#line 178 "" + goto L80; +#line 178 "" + } +/* SAVE THE CURRENT OPTIMAL SOLUTION. */ +/*< 250 IF ( Z .GE. PROFIT ) GO TO 270 >*/ +#line 180 "" +L250: +#line 180 "" + if (*z__ >= profit) { +#line 180 "" + goto L270; +#line 180 "" + } +/*< Z = PROFIT >*/ +#line 181 "" + *z__ = profit; +/*< DO 260 J=1,N >*/ +#line 182 "" + i__1 = *n; +#line 182 "" + for (j = 1; j <= i__1; ++j) { +/*< X(J) = XX(J) >*/ +#line 183 "" + x[j] = xx[j]; +/*< 260 CONTINUE >*/ +#line 184 "" +/* L260: */ +#line 184 "" + } +/*< IF ( Z .EQ. LIM ) RETURN >*/ +#line 185 "" + if (*z__ == lim) { +#line 185 "" + return 0; +#line 185 "" + } +/*< 270 IF ( XX(N) .EQ. 0 ) GO TO 280 >*/ +#line 186 "" +L270: +#line 186 "" + if (xx[*n] == 0) { +#line 186 "" + goto L280; +#line 186 "" + } +/*< XX(N) = 0 >*/ +#line 187 "" + xx[*n] = 0; +/*< CH = CH + W(N) >*/ +#line 188 "" + ch += w[*n]; +/*< PROFIT = PROFIT - P(N) >*/ +#line 189 "" + profit -= p[*n]; +/* BACKTRACK. */ +/*< 280 NN = II - 1 >*/ +#line 191 "" +L280: +#line 191 "" + nn = ii - 1; +/*< IF ( NN .EQ. 0 ) RETURN >*/ +#line 192 "" + if (nn == 0) { +#line 192 "" + return 0; +#line 192 "" + } +/*< DO 290 J=1,NN >*/ +#line 193 "" + i__1 = nn; +#line 193 "" + for (j = 1; j <= i__1; ++j) { +/*< KK = II - J >*/ +#line 194 "" + kk = ii - j; +/*< IF ( XX(KK) .EQ. 1 ) GO TO 300 >*/ +#line 195 "" + if (xx[kk] == 1) { +#line 195 "" + goto L300; +#line 195 "" + } +/*< 290 CONTINUE >*/ +#line 196 "" +/* L290: */ +#line 196 "" + } +/*< RETURN >*/ +#line 197 "" + return 0; +/*< 300 R = CH >*/ +#line 198 "" +L300: +#line 198 "" + r__ = ch; +/*< CH = CH + W(KK) >*/ +#line 199 "" + ch += w[kk]; +/*< PROFIT = PROFIT - P(KK) >*/ +#line 200 "" + profit -= p[kk]; +/*< XX(KK) = 0 >*/ +#line 201 "" + xx[kk] = 0; +/*< IF ( R .LT. MIN(KK) ) GO TO 310 >*/ +#line 202 "" + if (r__ < min__[kk]) { +#line 202 "" + goto L310; +#line 202 "" + } +/*< II = KK + 1 >*/ +#line 203 "" + ii = kk + 1; +/*< GO TO 80 >*/ +#line 204 "" + goto L80; +/*< 310 NN = KK + 1 >*/ +#line 205 "" +L310: +#line 205 "" + nn = kk + 1; +/*< II = KK >*/ +#line 206 "" + ii = kk; +/* TRY TO SUBSTITUTE THE NN-TH ITEM FOR THE KK-TH. */ +/*< 320 IF ( Z .GE. PROFIT + CH*P(NN)/W(NN) ) GO TO 280 >*/ +#line 208 "" +L320: +#line 208 "" + if (*z__ >= profit + ch * p[nn] / w[nn]) { +#line 208 "" + goto L280; +#line 208 "" + } +/*< DIFF = W(NN) - W(KK) >*/ +#line 209 "" + diff = w[nn] - w[kk]; +/*< IF ( DIFF ) 370, 330, 340 >*/ +#line 210 "" + if (diff < 0) { +#line 210 "" + goto L370; +#line 210 "" + } else if (diff == 0) { +#line 210 "" + goto L330; +#line 210 "" + } else { +#line 210 "" + goto L340; +#line 210 "" + } +/*< 330 NN = NN + 1 >*/ +#line 211 "" +L330: +#line 211 "" + ++nn; +/*< GO TO 320 >*/ +#line 212 "" + goto L320; +/*< 340 IF ( DIFF .GT. R ) GO TO 330 >*/ +#line 213 "" +L340: +#line 213 "" + if (diff > r__) { +#line 213 "" + goto L330; +#line 213 "" + } +/*< IF ( Z .GE. PROFIT + P(NN) ) GO TO 330 >*/ +#line 214 "" + if (*z__ >= profit + p[nn]) { +#line 214 "" + goto L330; +#line 214 "" + } +/*< Z = PROFIT + P(NN) >*/ +#line 215 "" + *z__ = profit + p[nn]; +/*< DO 350 J=1,KK >*/ +#line 216 "" + i__1 = kk; +#line 216 "" + for (j = 1; j <= i__1; ++j) { +/*< X(J) = XX(J) >*/ +#line 217 "" + x[j] = xx[j]; +/*< 350 CONTINUE >*/ +#line 218 "" +/* L350: */ +#line 218 "" + } +/*< JJ = KK + 1 >*/ +#line 219 "" + jj = kk + 1; +/*< DO 360 J=JJ,N >*/ +#line 220 "" + i__1 = *n; +#line 220 "" + for (j = jj; j <= i__1; ++j) { +/*< X(J) = 0 >*/ +#line 221 "" + x[j] = 0; +/*< 360 CONTINUE >*/ +#line 222 "" +/* L360: */ +#line 222 "" + } +/*< X(NN) = 1 >*/ +#line 223 "" + x[nn] = 1; +/*< IF ( Z .EQ. LIM ) RETURN >*/ +#line 224 "" + if (*z__ == lim) { +#line 224 "" + return 0; +#line 224 "" + } +/*< R = R - DIFF >*/ +#line 225 "" + r__ -= diff; +/*< KK = NN >*/ +#line 226 "" + kk = nn; +/*< NN = NN + 1 >*/ +#line 227 "" + ++nn; +/*< GO TO 320 >*/ +#line 228 "" + goto L320; +/*< 370 T = R - DIFF >*/ +#line 229 "" +L370: +#line 229 "" + t = r__ - diff; +/*< IF ( T .LT. MIN(NN) ) GO TO 330 >*/ +#line 230 "" + if (t < min__[nn]) { +#line 230 "" + goto L330; +#line 230 "" + } +/*< IF ( Z .GE. PROFIT + P(NN) + T*P(NN+1)/W(NN+1)) GO TO 280 >*/ +#line 231 "" + if (*z__ >= profit + p[nn] + t * p[nn + 1] / w[nn + 1]) { +#line 231 "" + goto L280; +#line 231 "" + } +/*< CH = CH - W(NN) >*/ +#line 232 "" + ch -= w[nn]; +/*< PROFIT = PROFIT + P(NN) >*/ +#line 233 "" + profit += p[nn]; +/*< XX(NN) = 1 >*/ +#line 234 "" + xx[nn] = 1; +/*< II = NN + 1 >*/ +#line 235 "" + ii = nn + 1; +/*< WSIGN(NN) = W(NN) >*/ +#line 236 "" + wsign[nn] = w[nn]; +/*< PSIGN(NN) = P(NN) >*/ +#line 237 "" + psign[nn] = p[nn]; +/*< ZSIGN(NN) = II >*/ +#line 238 "" + zsign[nn] = ii; +/*< N1 = NN + 1 >*/ +#line 239 "" + n1 = nn + 1; +/*< DO 380 J=N1,LOLD >*/ +#line 240 "" + i__1 = lold; +#line 240 "" + for (j = n1; j <= i__1; ++j) { +/*< WSIGN(J) = 0 >*/ +#line 241 "" + wsign[j] = 0; +/*< PSIGN(J) = 0 >*/ +#line 242 "" + psign[j] = 0; +/*< ZSIGN(J) = J >*/ +#line 243 "" + zsign[j] = j; +/*< 380 CONTINUE >*/ +#line 244 "" +/* L380: */ +#line 244 "" + } +/*< LOLD = NN >*/ +#line 245 "" + lold = nn; +/*< GO TO 80 >*/ +#line 246 "" + goto L80; +/*< END >*/ +} /* mt1_ */ + +/*< SUBROUTINE CHMT1(N,P,W,C,Z,JDIM) >*/ +#if 1 /* by mao */ +static +#endif +/* Subroutine */ int chmt1_(integer *n, integer *p, integer *w, integer *c__, + integer *z__, integer *jdim) +{ + /* System generated locals */ + integer i__1; + + /* Local variables */ + static integer j; + static real r__, rr; + static integer jsw; + + +/* CHECK THE INPUT DATA. */ + +/*< INTEGER P(JDIM),W(JDIM),C,Z >*/ +/*< IF ( N .GE. 2 .AND. N .LE. JDIM - 1 ) GO TO 10 >*/ +#line 253 "" + /* Parameter adjustments */ +#line 253 "" + --w; +#line 253 "" + --p; +#line 253 "" + +#line 253 "" + /* Function Body */ +#line 253 "" + if (*n >= 2 && *n <= *jdim - 1) { +#line 253 "" + goto L10; +#line 253 "" + } +/*< Z = - 1 >*/ +#line 254 "" + *z__ = -1; +/*< RETURN >*/ +#line 255 "" + return 0; +/*< 10 IF ( C .GT. 0 ) GO TO 30 >*/ +#line 256 "" +L10: +#line 256 "" + if (*c__ > 0) { +#line 256 "" + goto L30; +#line 256 "" + } +/*< 20 Z = - 2 >*/ +#line 257 "" +L20: +#line 257 "" + *z__ = -2; +/*< RETURN >*/ +#line 258 "" + return 0; +/*< 30 JSW = 0 >*/ +#line 259 "" +L30: +#line 259 "" + jsw = 0; +/*< RR = FLOAT(P(1))/FLOAT(W(1)) >*/ +#line 260 "" + rr = (real) p[1] / (real) w[1]; +/*< DO 50 J=1,N >*/ +#line 261 "" + i__1 = *n; +#line 261 "" + for (j = 1; j <= i__1; ++j) { +/*< R = RR >*/ +#line 262 "" + r__ = rr; +/*< IF ( P(J) .LE. 0 ) GO TO 20 >*/ +#line 263 "" + if (p[j] <= 0) { +#line 263 "" + goto L20; +#line 263 "" + } +/*< IF ( W(J) .LE. 0 ) GO TO 20 >*/ +#line 264 "" + if (w[j] <= 0) { +#line 264 "" + goto L20; +#line 264 "" + } +/*< JSW = JSW + W(J) >*/ +#line 265 "" + jsw += w[j]; +/*< IF ( W(J) .LE. C ) GO TO 40 >*/ +#line 266 "" + if (w[j] <= *c__) { +#line 266 "" + goto L40; +#line 266 "" + } +/*< Z = - 3 >*/ +#line 267 "" + *z__ = -3; +/*< RETURN >*/ +#line 268 "" + return 0; +/*< 40 RR = FLOAT(P(J))/FLOAT(W(J)) >*/ +#line 269 "" +L40: +#line 269 "" + rr = (real) p[j] / (real) w[j]; +/*< IF ( RR .LE. R ) GO TO 50 >*/ +#line 270 "" + if (rr <= r__) { +#line 270 "" + goto L50; +#line 270 "" + } +/*< Z = - 5 >*/ +#line 271 "" + *z__ = -5; +/*< RETURN >*/ +#line 272 "" + return 0; +/*< 50 CONTINUE >*/ +#line 273 "" +L50: +#line 273 "" + ; +#line 273 "" + } +/*< IF ( JSW .GT. C ) RETURN >*/ +#line 274 "" + if (jsw > *c__) { +#line 274 "" + return 0; +#line 274 "" + } +/*< Z = - 4 >*/ +#line 275 "" + *z__ = -4; +/*< RETURN >*/ +#line 276 "" + return 0; +/*< END >*/ +} /* chmt1_ */ + +#if 1 /* by mao */ +int mt1(int n, int p[], int w[], int c, int x[], int jck, int xx[], + int min[], int psign[], int wsign[], int zsign[]) +{ /* solve 0-1 knapsack problem */ + int z, jdim = n+1, j, s1, s2; + mt1_(&n, &p[1], &w[1], &c, &z, &x[1], &jdim, &jck, &xx[1], + &min[1], &psign[1], &wsign[1], &zsign[1]); + /* check solution found */ + s1 = s2 = 0; + for (j = 1; j <= n; j++) + { xassert(x[j] == 0 || x[j] == 1); + if (x[j]) + s1 += p[j], s2 += w[j]; + } + xassert(s1 == z); + xassert(s2 <= c); + return z; +} +#endif + +/* eof */ diff --git a/src/glpk/misc/mt1.h b/src/glpk/misc/mt1.h new file mode 100644 index 0000000..cceebba --- /dev/null +++ b/src/glpk/misc/mt1.h @@ -0,0 +1,34 @@ +/* mt1.h (0-1 knapsack problem; Martello & Toth algorithm) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2017-2018 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#ifndef MT1_H +#define MT1_H + +#define mt1 _glp_mt1 +int mt1(int n, int p[], int w[], int c, int x[], int jck, int xx[], + int min[], int psign[], int wsign[], int zsign[]); +/* solve 0-1 single knapsack problem */ + +#endif + +/* eof */ diff --git a/src/glpk/misc/wclique1.c b/src/glpk/misc/wclique1.c index a6fa960..a3d8954 100644 --- a/src/glpk/misc/wclique1.c +++ b/src/glpk/misc/wclique1.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2012-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2012-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "env.h" #include "wclique1.h" @@ -68,11 +64,7 @@ struct vertex { int i; double cw; }; -#ifndef _MSC_VER -static int fcmp(const void *xx, const void *yy) -#else -static int __cdecl fcmp(const void *xx, const void *yy) -#endif +static int CDECL fcmp(const void *xx, const void *yy) { const struct vertex *x = xx, *y = yy; if (x->cw > y->cw) return -1; if (x->cw < y->cw) return +1; diff --git a/src/glpk/glpnpp.h b/src/glpk/npp/npp.h similarity index 98% rename from src/glpk/glpnpp.h rename to src/glpk/npp/npp.h index 6aaebe8..428cb23 100644 --- a/src/glpk/glpnpp.h +++ b/src/glpk/npp/npp.h @@ -1,10 +1,9 @@ -/* glpnpp.h (LP/MIP preprocessor) */ +/* npp.h (LP/MIP preprocessor) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -22,19 +21,27 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifndef GLPNPP_H -#define GLPNPP_H +#ifndef NPP_H +#define NPP_H #include "prob.h" +#if 0 /* 20/XI-2017 */ typedef struct NPP NPP; +#else +typedef struct glp_prep NPP; +#endif typedef struct NPPROW NPPROW; typedef struct NPPCOL NPPCOL; typedef struct NPPAIJ NPPAIJ; typedef struct NPPTSE NPPTSE; typedef struct NPPLFE NPPLFE; +#if 0 /* 20/XI-2017 */ struct NPP +#else +struct glp_prep +#endif { /* LP/MIP preprocessor workspace */ /*--------------------------------------------------------------*/ /* original problem segment */ diff --git a/src/glpk/glpnpp01.c b/src/glpk/npp/npp1.c similarity index 99% rename from src/glpk/glpnpp01.c rename to src/glpk/npp/npp1.c index 9eb59f2..51758ba 100644 --- a/src/glpk/glpnpp01.c +++ b/src/glpk/npp/npp1.c @@ -1,10 +1,9 @@ -/* glpnpp01.c */ +/* npp1.c */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpnpp.h" +#include "npp.h" NPP *npp_create_wksp(void) { /* create LP/MIP preprocessor workspace */ diff --git a/src/glpk/glpnpp02.c b/src/glpk/npp/npp2.c similarity index 99% rename from src/glpk/glpnpp02.c rename to src/glpk/npp/npp2.c index ec4a455..4efcf1d 100644 --- a/src/glpk/glpnpp02.c +++ b/src/glpk/npp/npp2.c @@ -1,10 +1,9 @@ -/* glpnpp02.c */ +/* npp2.c */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpnpp.h" +#include "npp.h" /*********************************************************************** * NAME diff --git a/src/glpk/glpnpp03.c b/src/glpk/npp/npp3.c similarity index 99% rename from src/glpk/glpnpp03.c rename to src/glpk/npp/npp3.c index 0c869ee..883af12 100644 --- a/src/glpk/glpnpp03.c +++ b/src/glpk/npp/npp3.c @@ -1,10 +1,9 @@ -/* glpnpp03.c */ +/* npp3.c */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpnpp.h" +#include "npp.h" /*********************************************************************** * NAME diff --git a/src/glpk/glpnpp04.c b/src/glpk/npp/npp4.c similarity index 99% rename from src/glpk/glpnpp04.c rename to src/glpk/npp/npp4.c index 925b1ba..d7dd0e8 100644 --- a/src/glpk/glpnpp04.c +++ b/src/glpk/npp/npp4.c @@ -1,10 +1,9 @@ -/* glpnpp04.c */ +/* npp4.c */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpnpp.h" +#include "npp.h" /*********************************************************************** * NAME diff --git a/src/glpk/glpnpp05.c b/src/glpk/npp/npp5.c similarity index 99% rename from src/glpk/glpnpp05.c rename to src/glpk/npp/npp5.c index 8c88185..2fad496 100644 --- a/src/glpk/glpnpp05.c +++ b/src/glpk/npp/npp5.c @@ -1,10 +1,9 @@ -/* glpnpp05.c */ +/* npp5.c */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, -* 2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2009-2017 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -23,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpnpp.h" +#include "npp.h" /*********************************************************************** * NAME diff --git a/src/glpk/proxy/proxy1.c b/src/glpk/proxy/proxy1.c index f657864..5f9850d 100644 --- a/src/glpk/proxy/proxy1.c +++ b/src/glpk/proxy/proxy1.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2013 Andrew Makhorin, Department for Applied +* Copyright (C) 2013, 2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -22,7 +22,7 @@ ***********************************************************************/ #include "env.h" -#include "glpios.h" +#include "ios.h" #include "proxy.h" void ios_proxy_heur(glp_tree *T) diff --git a/src/glpk/simplex/spxchuzr.c b/src/glpk/simplex/spxchuzr.c index 471928b..8bef77b 100644 --- a/src/glpk/simplex/spxchuzr.c +++ b/src/glpk/simplex/spxchuzr.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2015-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2015-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -536,11 +536,7 @@ int spx_ls_eval_bp(SPXLP *lp, const double beta[/*1+m*/], * On exit the routine also replaces the parameter slope with a new * value that corresponds to the new last break-point bp[num1]. */ -#ifndef _MSC_VER -static int fcmp(const void *v1, const void *v2) -#else -static int __cdecl fcmp(const void *v1, const void *v2) -#endif +static int CDECL fcmp(const void *v1, const void *v2) { const SPXBP *p1 = v1, *p2 = v2; if (p1->teta < p2->teta) return -1; diff --git a/src/glpk/simplex/spychuzc.c b/src/glpk/simplex/spychuzc.c index a209b0c..b922129 100644 --- a/src/glpk/simplex/spychuzc.c +++ b/src/glpk/simplex/spychuzc.c @@ -3,7 +3,7 @@ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * -* Copyright (C) 2015-2017 Andrew Makhorin, Department for Applied +* Copyright (C) 2015-2018 Andrew Makhorin, Department for Applied * Informatics, Moscow Aviation Institute, Moscow, Russia. All rights * reserved. E-mail: . * @@ -21,10 +21,6 @@ * along with GLPK. If not, see . ***********************************************************************/ -#ifdef HAVE_CONFIG_H -#include -#endif - #include "env.h" #include "spychuzc.h" @@ -279,11 +275,7 @@ done: return q; * array bp only the break-points that correspond to positive increment * of the dual objective. */ -#ifndef _MSC_VER -static int fcmp(const void *v1, const void *v2) -#else -static int __cdecl fcmp(const void *v1, const void *v2) -#endif +static int CDECL fcmp(const void *v1, const void *v2) { const SPYBP *p1 = v1, *p2 = v2; if (p1->teta < p2->teta) return -1; @@ -505,11 +497,7 @@ int spy_ls_eval_bp(SPXLP *lp, const double d[/*1+n-m*/], * On exit the routine also replaces the parameter slope with a new * value that corresponds to the new last break-point bp[num1]. */ -#ifndef _MSC_VER -static int fcmp(const void *v1, const void *v2) -#else -static int __cdecl fcmp(const void *v1, const void *v2) -#endif +static int CDECL fcmp(const void *v1, const void *v2) { const SPYBP *p1 = v1, *p2 = v2; if (p1->teta < p2->teta) return -1; diff --git a/src/glpk/simplex/spydual.c b/src/glpk/simplex/spydual.c index 04d918f..89d98db 100644 --- a/src/glpk/simplex/spydual.c +++ b/src/glpk/simplex/spydual.c @@ -914,13 +914,10 @@ static void play_coef(struct csa *csa, int all) char *flag = lp->flag; double *orig_c = csa->orig_c; double *d = csa->d; -#if 0 /* 29/III-2016 */ - double *trow = csa->trow; /* was used to update d = (d[j]) */ -#else const double *trow = csa->trow.vec; -#endif + /* this vector was used to update d = (d[j]) */ int j, k; - double temp; + static const double eps = 1e-9; /* reduced costs d = (d[j]) should be valid */ xassert(csa->d_st); /* walk thru the list of non-basic variables xN = (xN[j]) */ @@ -933,45 +930,33 @@ static void play_coef(struct csa *csa, int all) /* d[j] may have any sign */ } else if (l[k] == -DBL_MAX && u[k] == +DBL_MAX) - { /* xN[j] is free variable */ + { /* xN[j] is free (unbounded) variable */ /* strong feasibility means d[j] = 0 */ c[k] -= d[j], d[j] = 0.0; + /* in this case dual degeneracy is not critical, since + * if xN[j] enters the basis, it never leaves it */ } else if (!flag[j]) { /* xN[j] has its lower bound active */ - /* strong feasibility means d[j] >= 0 */ - if (d[j] < 0.0) -#if 1 /* 12/VII-2017 */ - c[k] -= d[j], d[j] = 0.0; -#else - c[k] -= d[j] - 1e-9, d[j] = +1e-9; -#endif - else if (c[k] > orig_c[k]) - { /* remove/reduce perturbation of cN[j] */ - temp = c[k] - orig_c[k]; /* > 0 */ - if (temp < d[j]) - c[k] = orig_c[k], d[j] -= temp; - else - c[k] -= d[j], d[j] = 0.0; - } + xassert(l[k] != -DBL_MAX); + /* first, we remove current perturbation to provide + * c[k] = orig_c[k] */ + d[j] -= c[k] - orig_c[k], c[k] = orig_c[k]; + /* strong feasibility means d[j] >= 0, but we provide + * d[j] >= +eps to prevent dual degeneracy */ + if (d[j] < +eps) + c[k] -= d[j] - eps, d[j] = +eps; } else { /* xN[j] has its upper bound active */ - /* strong feasibility means d[j] <= 0 */ - if (d[j] > 0.0) -#if 1 /* 12/VII-2017 */ - c[k] -= d[j], d[j] = 0.0; -#else - c[k] -= d[j] + 1e-9, d[j] = -1e-9; -#endif - else if (c[k] < orig_c[k]) - { /* remove/reduce perturbation of cN[j] */ - temp = c[k] - orig_c[k]; /* < 0 */ - if (temp > d[j]) - c[k] = orig_c[k], d[j] -= temp; - else - c[k] -= d[j], d[j] = 0.0; - } + xassert(u[k] != +DBL_MAX); + /* similarly, we remove current perturbation to provide + * c[k] = orig_c[k] */ + d[j] -= c[k] - orig_c[k], c[k] = orig_c[k]; + /* strong feasibility means d[j] <= 0, but we provide + * d[j] <= -eps to prevent dual degeneracy */ + if (d[j] > -eps) + c[k] -= d[j] + eps, d[j] = -eps; } } }