lemon-project-template-glpk

annotate deps/glpk/src/glpapi20.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpapi20.c */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpnpp.h"
alpar@9 26
alpar@9 27 int glp_intfeas1(glp_prob *P, int use_bound, int obj_bound)
alpar@9 28 { /* solve integer feasibility problem */
alpar@9 29 NPP *npp = NULL;
alpar@9 30 glp_prob *mip = NULL;
alpar@9 31 int *obj_ind = NULL;
alpar@9 32 double *obj_val = NULL;
alpar@9 33 int obj_row = 0;
alpar@9 34 int i, j, k, obj_len, temp, ret;
alpar@9 35 /* check the problem object */
alpar@9 36 if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@9 37 xerror("glp_intfeas1: P = %p; invalid problem object\n",
alpar@9 38 P);
alpar@9 39 if (P->tree != NULL)
alpar@9 40 xerror("glp_intfeas1: operation not allowed\n");
alpar@9 41 /* integer solution is currently undefined */
alpar@9 42 P->mip_stat = GLP_UNDEF;
alpar@9 43 P->mip_obj = 0.0;
alpar@9 44 /* check columns (variables) */
alpar@9 45 for (j = 1; j <= P->n; j++)
alpar@9 46 { GLPCOL *col = P->col[j];
alpar@9 47 #if 0 /* currently binarization is not yet implemented */
alpar@9 48 if (!(col->kind == GLP_IV || col->type == GLP_FX))
alpar@9 49 { xprintf("glp_intfeas1: column %d: non-integer non-fixed var"
alpar@9 50 "iable not allowed\n", j);
alpar@9 51 #else
alpar@9 52 if (!((col->kind == GLP_IV && col->lb == 0.0 && col->ub == 1.0)
alpar@9 53 || col->type == GLP_FX))
alpar@9 54 { xprintf("glp_intfeas1: column %d: non-binary non-fixed vari"
alpar@9 55 "able not allowed\n", j);
alpar@9 56 #endif
alpar@9 57 ret = GLP_EDATA;
alpar@9 58 goto done;
alpar@9 59 }
alpar@9 60 temp = (int)col->lb;
alpar@9 61 if ((double)temp != col->lb)
alpar@9 62 { if (col->type == GLP_FX)
alpar@9 63 xprintf("glp_intfeas1: column %d: fixed value %g is non-"
alpar@9 64 "integer or out of range\n", j, col->lb);
alpar@9 65 else
alpar@9 66 xprintf("glp_intfeas1: column %d: lower bound %g is non-"
alpar@9 67 "integer or out of range\n", j, col->lb);
alpar@9 68 ret = GLP_EDATA;
alpar@9 69 goto done;
alpar@9 70 }
alpar@9 71 temp = (int)col->ub;
alpar@9 72 if ((double)temp != col->ub)
alpar@9 73 { xprintf("glp_intfeas1: column %d: upper bound %g is non-int"
alpar@9 74 "eger or out of range\n", j, col->ub);
alpar@9 75 ret = GLP_EDATA;
alpar@9 76 goto done;
alpar@9 77 }
alpar@9 78 if (col->type == GLP_DB && col->lb > col->ub)
alpar@9 79 { xprintf("glp_intfeas1: column %d: lower bound %g is greater"
alpar@9 80 " than upper bound %g\n", j, col->lb, col->ub);
alpar@9 81 ret = GLP_EBOUND;
alpar@9 82 goto done;
alpar@9 83 }
alpar@9 84 }
alpar@9 85 /* check rows (constraints) */
alpar@9 86 for (i = 1; i <= P->m; i++)
alpar@9 87 { GLPROW *row = P->row[i];
alpar@9 88 GLPAIJ *aij;
alpar@9 89 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 90 { temp = (int)aij->val;
alpar@9 91 if ((double)temp != aij->val)
alpar@9 92 { xprintf("glp_intfeas1: row = %d, column %d: constraint c"
alpar@9 93 "oefficient %g is non-integer or out of range\n",
alpar@9 94 i, aij->col->j, aij->val);
alpar@9 95 ret = GLP_EDATA;
alpar@9 96 goto done;
alpar@9 97 }
alpar@9 98 }
alpar@9 99 temp = (int)row->lb;
alpar@9 100 if ((double)temp != row->lb)
alpar@9 101 { if (row->type == GLP_FX)
alpar@9 102 xprintf("glp_intfeas1: row = %d: fixed value %g is non-i"
alpar@9 103 "nteger or out of range\n", i, row->lb);
alpar@9 104 else
alpar@9 105 xprintf("glp_intfeas1: row = %d: lower bound %g is non-i"
alpar@9 106 "nteger or out of range\n", i, row->lb);
alpar@9 107 ret = GLP_EDATA;
alpar@9 108 goto done;
alpar@9 109 }
alpar@9 110 temp = (int)row->ub;
alpar@9 111 if ((double)temp != row->ub)
alpar@9 112 { xprintf("glp_intfeas1: row = %d: upper bound %g is non-inte"
alpar@9 113 "ger or out of range\n", i, row->ub);
alpar@9 114 ret = GLP_EDATA;
alpar@9 115 goto done;
alpar@9 116 }
alpar@9 117 if (row->type == GLP_DB && row->lb > row->ub)
alpar@9 118 { xprintf("glp_intfeas1: row %d: lower bound %g is greater th"
alpar@9 119 "an upper bound %g\n", i, row->lb, row->ub);
alpar@9 120 ret = GLP_EBOUND;
alpar@9 121 goto done;
alpar@9 122 }
alpar@9 123 }
alpar@9 124 /* check the objective function */
alpar@9 125 temp = (int)P->c0;
alpar@9 126 if ((double)temp != P->c0)
alpar@9 127 { xprintf("glp_intfeas1: objective constant term %g is non-integ"
alpar@9 128 "er or out of range\n", P->c0);
alpar@9 129 ret = GLP_EDATA;
alpar@9 130 goto done;
alpar@9 131 }
alpar@9 132 for (j = 1; j <= P->n; j++)
alpar@9 133 { temp = (int)P->col[j]->coef;
alpar@9 134 if ((double)temp != P->col[j]->coef)
alpar@9 135 { xprintf("glp_intfeas1: column %d: objective coefficient is "
alpar@9 136 "non-integer or out of range\n", j, P->col[j]->coef);
alpar@9 137 ret = GLP_EDATA;
alpar@9 138 goto done;
alpar@9 139 }
alpar@9 140 }
alpar@9 141 /* save the objective function and set it to zero */
alpar@9 142 obj_ind = xcalloc(1+P->n, sizeof(int));
alpar@9 143 obj_val = xcalloc(1+P->n, sizeof(double));
alpar@9 144 obj_len = 0;
alpar@9 145 obj_ind[0] = 0;
alpar@9 146 obj_val[0] = P->c0;
alpar@9 147 P->c0 = 0.0;
alpar@9 148 for (j = 1; j <= P->n; j++)
alpar@9 149 { if (P->col[j]->coef != 0.0)
alpar@9 150 { obj_len++;
alpar@9 151 obj_ind[obj_len] = j;
alpar@9 152 obj_val[obj_len] = P->col[j]->coef;
alpar@9 153 P->col[j]->coef = 0.0;
alpar@9 154 }
alpar@9 155 }
alpar@9 156 /* add inequality to bound the objective function, if required */
alpar@9 157 if (!use_bound)
alpar@9 158 xprintf("Will search for ANY feasible solution\n");
alpar@9 159 else
alpar@9 160 { xprintf("Will search only for solution not worse than %d\n",
alpar@9 161 obj_bound);
alpar@9 162 obj_row = glp_add_rows(P, 1);
alpar@9 163 glp_set_mat_row(P, obj_row, obj_len, obj_ind, obj_val);
alpar@9 164 if (P->dir == GLP_MIN)
alpar@9 165 glp_set_row_bnds(P, obj_row,
alpar@9 166 GLP_UP, 0.0, (double)obj_bound - obj_val[0]);
alpar@9 167 else if (P->dir == GLP_MAX)
alpar@9 168 glp_set_row_bnds(P, obj_row,
alpar@9 169 GLP_LO, (double)obj_bound - obj_val[0], 0.0);
alpar@9 170 else
alpar@9 171 xassert(P != P);
alpar@9 172 }
alpar@9 173 /* create preprocessor workspace */
alpar@9 174 xprintf("Translating to CNF-SAT...\n");
alpar@9 175 xprintf("Original problem has %d row%s, %d column%s, and %d non-z"
alpar@9 176 "ero%s\n", P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" :
alpar@9 177 "s", P->nnz, P->nnz == 1 ? "" : "s");
alpar@9 178 npp = npp_create_wksp();
alpar@9 179 /* load the original problem into the preprocessor workspace */
alpar@9 180 npp_load_prob(npp, P, GLP_OFF, GLP_MIP, GLP_OFF);
alpar@9 181 /* perform translation to SAT-CNF problem instance */
alpar@9 182 ret = npp_sat_encode_prob(npp);
alpar@9 183 if (ret == 0)
alpar@9 184 ;
alpar@9 185 else if (ret == GLP_ENOPFS)
alpar@9 186 xprintf("PROBLEM HAS NO INTEGER FEASIBLE SOLUTION\n");
alpar@9 187 else if (ret == GLP_ERANGE)
alpar@9 188 xprintf("glp_intfeas1: translation to SAT-CNF failed because o"
alpar@9 189 "f integer overflow\n");
alpar@9 190 else
alpar@9 191 xassert(ret != ret);
alpar@9 192 if (ret != 0)
alpar@9 193 goto done;
alpar@9 194 /* build SAT-CNF problem instance and try to solve it */
alpar@9 195 mip = glp_create_prob();
alpar@9 196 npp_build_prob(npp, mip);
alpar@9 197 ret = glp_minisat1(mip);
alpar@9 198 /* only integer feasible solution can be postprocessed */
alpar@9 199 if (!(mip->mip_stat == GLP_OPT || mip->mip_stat == GLP_FEAS))
alpar@9 200 { P->mip_stat = mip->mip_stat;
alpar@9 201 goto done;
alpar@9 202 }
alpar@9 203 /* postprocess the solution found */
alpar@9 204 npp_postprocess(npp, mip);
alpar@9 205 /* the transformed problem is no longer needed */
alpar@9 206 glp_delete_prob(mip), mip = NULL;
alpar@9 207 /* store solution to the original problem object */
alpar@9 208 npp_unload_sol(npp, P);
alpar@9 209 /* change the solution status to 'integer feasible' */
alpar@9 210 P->mip_stat = GLP_FEAS;
alpar@9 211 /* check integer feasibility */
alpar@9 212 for (i = 1; i <= P->m; i++)
alpar@9 213 { GLPROW *row;
alpar@9 214 GLPAIJ *aij;
alpar@9 215 double sum;
alpar@9 216 row = P->row[i];
alpar@9 217 sum = 0.0;
alpar@9 218 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 219 sum += aij->val * aij->col->mipx;
alpar@9 220 xassert(sum == row->mipx);
alpar@9 221 if (row->type == GLP_LO || row->type == GLP_DB ||
alpar@9 222 row->type == GLP_FX)
alpar@9 223 xassert(sum >= row->lb);
alpar@9 224 if (row->type == GLP_UP || row->type == GLP_DB ||
alpar@9 225 row->type == GLP_FX)
alpar@9 226 xassert(sum <= row->ub);
alpar@9 227 }
alpar@9 228 /* compute value of the original objective function */
alpar@9 229 P->mip_obj = obj_val[0];
alpar@9 230 for (k = 1; k <= obj_len; k++)
alpar@9 231 P->mip_obj += obj_val[k] * P->col[obj_ind[k]]->mipx;
alpar@9 232 xprintf("Objective value = %17.9e\n", P->mip_obj);
alpar@9 233 done: /* delete the transformed problem, if it exists */
alpar@9 234 if (mip != NULL)
alpar@9 235 glp_delete_prob(mip);
alpar@9 236 /* delete the preprocessor workspace, if it exists */
alpar@9 237 if (npp != NULL)
alpar@9 238 npp_delete_wksp(npp);
alpar@9 239 /* remove inequality used to bound the objective function */
alpar@9 240 if (obj_row > 0)
alpar@9 241 { int ind[1+1];
alpar@9 242 ind[1] = obj_row;
alpar@9 243 glp_del_rows(P, 1, ind);
alpar@9 244 }
alpar@9 245 /* restore the original objective function */
alpar@9 246 if (obj_ind != NULL)
alpar@9 247 { P->c0 = obj_val[0];
alpar@9 248 for (k = 1; k <= obj_len; k++)
alpar@9 249 P->col[obj_ind[k]]->coef = obj_val[k];
alpar@9 250 xfree(obj_ind);
alpar@9 251 xfree(obj_val);
alpar@9 252 }
alpar@9 253 return ret;
alpar@9 254 }
alpar@9 255
alpar@9 256 /* eof */