lemon-project-template-glpk
comparison deps/glpk/src/glpios01.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 21:43:29 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:7530bd28198f |
---|---|
1 /* glpios01.c */ | |
2 | |
3 /*********************************************************************** | |
4 * This code is part of GLPK (GNU Linear Programming Kit). | |
5 * | |
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, | |
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, | |
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. | |
9 * E-mail: <mao@gnu.org>. | |
10 * | |
11 * GLPK is free software: you can redistribute it and/or modify it | |
12 * under the terms of the GNU General Public License as published by | |
13 * the Free Software Foundation, either version 3 of the License, or | |
14 * (at your option) any later version. | |
15 * | |
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT | |
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY | |
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public | |
19 * License for more details. | |
20 * | |
21 * You should have received a copy of the GNU General Public License | |
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. | |
23 ***********************************************************************/ | |
24 | |
25 #include "glpios.h" | |
26 | |
27 /*********************************************************************** | |
28 * NAME | |
29 * | |
30 * ios_create_tree - create branch-and-bound tree | |
31 * | |
32 * SYNOPSIS | |
33 * | |
34 * #include "glpios.h" | |
35 * glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm); | |
36 * | |
37 * DESCRIPTION | |
38 * | |
39 * The routine ios_create_tree creates the branch-and-bound tree. | |
40 * | |
41 * Being created the tree consists of the only root subproblem whose | |
42 * reference number is 1. Note that initially the root subproblem is in | |
43 * frozen state and therefore needs to be revived. | |
44 * | |
45 * RETURNS | |
46 * | |
47 * The routine returns a pointer to the tree created. */ | |
48 | |
49 static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent); | |
50 | |
51 glp_tree *ios_create_tree(glp_prob *mip, const glp_iocp *parm) | |
52 { int m = mip->m; | |
53 int n = mip->n; | |
54 glp_tree *tree; | |
55 int i, j; | |
56 xassert(mip->tree == NULL); | |
57 mip->tree = tree = xmalloc(sizeof(glp_tree)); | |
58 tree->pool = dmp_create_pool(); | |
59 tree->n = n; | |
60 /* save original problem components */ | |
61 tree->orig_m = m; | |
62 tree->orig_type = xcalloc(1+m+n, sizeof(char)); | |
63 tree->orig_lb = xcalloc(1+m+n, sizeof(double)); | |
64 tree->orig_ub = xcalloc(1+m+n, sizeof(double)); | |
65 tree->orig_stat = xcalloc(1+m+n, sizeof(char)); | |
66 tree->orig_prim = xcalloc(1+m+n, sizeof(double)); | |
67 tree->orig_dual = xcalloc(1+m+n, sizeof(double)); | |
68 for (i = 1; i <= m; i++) | |
69 { GLPROW *row = mip->row[i]; | |
70 tree->orig_type[i] = (char)row->type; | |
71 tree->orig_lb[i] = row->lb; | |
72 tree->orig_ub[i] = row->ub; | |
73 tree->orig_stat[i] = (char)row->stat; | |
74 tree->orig_prim[i] = row->prim; | |
75 tree->orig_dual[i] = row->dual; | |
76 } | |
77 for (j = 1; j <= n; j++) | |
78 { GLPCOL *col = mip->col[j]; | |
79 tree->orig_type[m+j] = (char)col->type; | |
80 tree->orig_lb[m+j] = col->lb; | |
81 tree->orig_ub[m+j] = col->ub; | |
82 tree->orig_stat[m+j] = (char)col->stat; | |
83 tree->orig_prim[m+j] = col->prim; | |
84 tree->orig_dual[m+j] = col->dual; | |
85 } | |
86 tree->orig_obj = mip->obj_val; | |
87 /* initialize the branch-and-bound tree */ | |
88 tree->nslots = 0; | |
89 tree->avail = 0; | |
90 tree->slot = NULL; | |
91 tree->head = tree->tail = NULL; | |
92 tree->a_cnt = tree->n_cnt = tree->t_cnt = 0; | |
93 /* the root subproblem is not solved yet, so its final components | |
94 are unknown so far */ | |
95 tree->root_m = 0; | |
96 tree->root_type = NULL; | |
97 tree->root_lb = tree->root_ub = NULL; | |
98 tree->root_stat = NULL; | |
99 /* the current subproblem does not exist yet */ | |
100 tree->curr = NULL; | |
101 tree->mip = mip; | |
102 /*tree->solved = 0;*/ | |
103 tree->non_int = xcalloc(1+n, sizeof(char)); | |
104 memset(&tree->non_int[1], 0, n); | |
105 /* arrays to save parent subproblem components will be allocated | |
106 later */ | |
107 tree->pred_m = tree->pred_max = 0; | |
108 tree->pred_type = NULL; | |
109 tree->pred_lb = tree->pred_ub = NULL; | |
110 tree->pred_stat = NULL; | |
111 /* cut generator */ | |
112 tree->local = ios_create_pool(tree); | |
113 /*tree->first_attempt = 1;*/ | |
114 /*tree->max_added_cuts = 0;*/ | |
115 /*tree->min_eff = 0.0;*/ | |
116 /*tree->miss = 0;*/ | |
117 /*tree->just_selected = 0;*/ | |
118 tree->mir_gen = NULL; | |
119 tree->clq_gen = NULL; | |
120 /*tree->round = 0;*/ | |
121 #if 0 | |
122 /* create the conflict graph */ | |
123 tree->n_ref = xcalloc(1+n, sizeof(int)); | |
124 memset(&tree->n_ref[1], 0, n * sizeof(int)); | |
125 tree->c_ref = xcalloc(1+n, sizeof(int)); | |
126 memset(&tree->c_ref[1], 0, n * sizeof(int)); | |
127 tree->g = scg_create_graph(0); | |
128 tree->j_ref = xcalloc(1+tree->g->n_max, sizeof(int)); | |
129 #endif | |
130 /* pseudocost branching */ | |
131 tree->pcost = NULL; | |
132 tree->iwrk = xcalloc(1+n, sizeof(int)); | |
133 tree->dwrk = xcalloc(1+n, sizeof(double)); | |
134 /* initialize control parameters */ | |
135 tree->parm = parm; | |
136 tree->tm_beg = xtime(); | |
137 tree->tm_lag = xlset(0); | |
138 tree->sol_cnt = 0; | |
139 /* initialize advanced solver interface */ | |
140 tree->reason = 0; | |
141 tree->reopt = 0; | |
142 tree->reinv = 0; | |
143 tree->br_var = 0; | |
144 tree->br_sel = 0; | |
145 tree->child = 0; | |
146 tree->next_p = 0; | |
147 /*tree->btrack = NULL;*/ | |
148 tree->stop = 0; | |
149 /* create the root subproblem, which initially is identical to | |
150 the original MIP */ | |
151 new_node(tree, NULL); | |
152 return tree; | |
153 } | |
154 | |
155 /*********************************************************************** | |
156 * NAME | |
157 * | |
158 * ios_revive_node - revive specified subproblem | |
159 * | |
160 * SYNOPSIS | |
161 * | |
162 * #include "glpios.h" | |
163 * void ios_revive_node(glp_tree *tree, int p); | |
164 * | |
165 * DESCRIPTION | |
166 * | |
167 * The routine ios_revive_node revives the specified subproblem, whose | |
168 * reference number is p, and thereby makes it the current subproblem. | |
169 * Note that the specified subproblem must be active. Besides, if the | |
170 * current subproblem already exists, it must be frozen before reviving | |
171 * another subproblem. */ | |
172 | |
173 void ios_revive_node(glp_tree *tree, int p) | |
174 { glp_prob *mip = tree->mip; | |
175 IOSNPD *node, *root; | |
176 /* obtain pointer to the specified subproblem */ | |
177 xassert(1 <= p && p <= tree->nslots); | |
178 node = tree->slot[p].node; | |
179 xassert(node != NULL); | |
180 /* the specified subproblem must be active */ | |
181 xassert(node->count == 0); | |
182 /* the current subproblem must not exist */ | |
183 xassert(tree->curr == NULL); | |
184 /* the specified subproblem becomes current */ | |
185 tree->curr = node; | |
186 /*tree->solved = 0;*/ | |
187 /* obtain pointer to the root subproblem */ | |
188 root = tree->slot[1].node; | |
189 xassert(root != NULL); | |
190 /* at this point problem object components correspond to the root | |
191 subproblem, so if the root subproblem should be revived, there | |
192 is nothing more to do */ | |
193 if (node == root) goto done; | |
194 xassert(mip->m == tree->root_m); | |
195 /* build path from the root to the current node */ | |
196 node->temp = NULL; | |
197 for (node = node; node != NULL; node = node->up) | |
198 { if (node->up == NULL) | |
199 xassert(node == root); | |
200 else | |
201 node->up->temp = node; | |
202 } | |
203 /* go down from the root to the current node and make necessary | |
204 changes to restore components of the current subproblem */ | |
205 for (node = root; node != NULL; node = node->temp) | |
206 { int m = mip->m; | |
207 int n = mip->n; | |
208 /* if the current node is reached, the problem object at this | |
209 point corresponds to its parent, so save attributes of rows | |
210 and columns for the parent subproblem */ | |
211 if (node->temp == NULL) | |
212 { int i, j; | |
213 tree->pred_m = m; | |
214 /* allocate/reallocate arrays, if necessary */ | |
215 if (tree->pred_max < m + n) | |
216 { int new_size = m + n + 100; | |
217 if (tree->pred_type != NULL) xfree(tree->pred_type); | |
218 if (tree->pred_lb != NULL) xfree(tree->pred_lb); | |
219 if (tree->pred_ub != NULL) xfree(tree->pred_ub); | |
220 if (tree->pred_stat != NULL) xfree(tree->pred_stat); | |
221 tree->pred_max = new_size; | |
222 tree->pred_type = xcalloc(1+new_size, sizeof(char)); | |
223 tree->pred_lb = xcalloc(1+new_size, sizeof(double)); | |
224 tree->pred_ub = xcalloc(1+new_size, sizeof(double)); | |
225 tree->pred_stat = xcalloc(1+new_size, sizeof(char)); | |
226 } | |
227 /* save row attributes */ | |
228 for (i = 1; i <= m; i++) | |
229 { GLPROW *row = mip->row[i]; | |
230 tree->pred_type[i] = (char)row->type; | |
231 tree->pred_lb[i] = row->lb; | |
232 tree->pred_ub[i] = row->ub; | |
233 tree->pred_stat[i] = (char)row->stat; | |
234 } | |
235 /* save column attributes */ | |
236 for (j = 1; j <= n; j++) | |
237 { GLPCOL *col = mip->col[j]; | |
238 tree->pred_type[mip->m+j] = (char)col->type; | |
239 tree->pred_lb[mip->m+j] = col->lb; | |
240 tree->pred_ub[mip->m+j] = col->ub; | |
241 tree->pred_stat[mip->m+j] = (char)col->stat; | |
242 } | |
243 } | |
244 /* change bounds of rows and columns */ | |
245 { IOSBND *b; | |
246 for (b = node->b_ptr; b != NULL; b = b->next) | |
247 { if (b->k <= m) | |
248 glp_set_row_bnds(mip, b->k, b->type, b->lb, b->ub); | |
249 else | |
250 glp_set_col_bnds(mip, b->k-m, b->type, b->lb, b->ub); | |
251 } | |
252 } | |
253 /* change statuses of rows and columns */ | |
254 { IOSTAT *s; | |
255 for (s = node->s_ptr; s != NULL; s = s->next) | |
256 { if (s->k <= m) | |
257 glp_set_row_stat(mip, s->k, s->stat); | |
258 else | |
259 glp_set_col_stat(mip, s->k-m, s->stat); | |
260 } | |
261 } | |
262 /* add new rows */ | |
263 if (node->r_ptr != NULL) | |
264 { IOSROW *r; | |
265 IOSAIJ *a; | |
266 int i, len, *ind; | |
267 double *val; | |
268 ind = xcalloc(1+n, sizeof(int)); | |
269 val = xcalloc(1+n, sizeof(double)); | |
270 for (r = node->r_ptr; r != NULL; r = r->next) | |
271 { i = glp_add_rows(mip, 1); | |
272 glp_set_row_name(mip, i, r->name); | |
273 #if 1 /* 20/IX-2008 */ | |
274 xassert(mip->row[i]->level == 0); | |
275 mip->row[i]->level = node->level; | |
276 mip->row[i]->origin = r->origin; | |
277 mip->row[i]->klass = r->klass; | |
278 #endif | |
279 glp_set_row_bnds(mip, i, r->type, r->lb, r->ub); | |
280 len = 0; | |
281 for (a = r->ptr; a != NULL; a = a->next) | |
282 len++, ind[len] = a->j, val[len] = a->val; | |
283 glp_set_mat_row(mip, i, len, ind, val); | |
284 glp_set_rii(mip, i, r->rii); | |
285 glp_set_row_stat(mip, i, r->stat); | |
286 } | |
287 xfree(ind); | |
288 xfree(val); | |
289 } | |
290 #if 0 | |
291 /* add new edges to the conflict graph */ | |
292 /* add new cliques to the conflict graph */ | |
293 /* (not implemented yet) */ | |
294 xassert(node->own_nn == 0); | |
295 xassert(node->own_nc == 0); | |
296 xassert(node->e_ptr == NULL); | |
297 #endif | |
298 } | |
299 /* the specified subproblem has been revived */ | |
300 node = tree->curr; | |
301 /* delete its bound change list */ | |
302 while (node->b_ptr != NULL) | |
303 { IOSBND *b; | |
304 b = node->b_ptr; | |
305 node->b_ptr = b->next; | |
306 dmp_free_atom(tree->pool, b, sizeof(IOSBND)); | |
307 } | |
308 /* delete its status change list */ | |
309 while (node->s_ptr != NULL) | |
310 { IOSTAT *s; | |
311 s = node->s_ptr; | |
312 node->s_ptr = s->next; | |
313 dmp_free_atom(tree->pool, s, sizeof(IOSTAT)); | |
314 } | |
315 #if 1 /* 20/XI-2009 */ | |
316 /* delete its row addition list (additional rows may appear, for | |
317 example, due to branching on GUB constraints */ | |
318 while (node->r_ptr != NULL) | |
319 { IOSROW *r; | |
320 r = node->r_ptr; | |
321 node->r_ptr = r->next; | |
322 xassert(r->name == NULL); | |
323 while (r->ptr != NULL) | |
324 { IOSAIJ *a; | |
325 a = r->ptr; | |
326 r->ptr = a->next; | |
327 dmp_free_atom(tree->pool, a, sizeof(IOSAIJ)); | |
328 } | |
329 dmp_free_atom(tree->pool, r, sizeof(IOSROW)); | |
330 } | |
331 #endif | |
332 done: return; | |
333 } | |
334 | |
335 /*********************************************************************** | |
336 * NAME | |
337 * | |
338 * ios_freeze_node - freeze current subproblem | |
339 * | |
340 * SYNOPSIS | |
341 * | |
342 * #include "glpios.h" | |
343 * void ios_freeze_node(glp_tree *tree); | |
344 * | |
345 * DESCRIPTION | |
346 * | |
347 * The routine ios_freeze_node freezes the current subproblem. */ | |
348 | |
349 void ios_freeze_node(glp_tree *tree) | |
350 { glp_prob *mip = tree->mip; | |
351 int m = mip->m; | |
352 int n = mip->n; | |
353 IOSNPD *node; | |
354 /* obtain pointer to the current subproblem */ | |
355 node = tree->curr; | |
356 xassert(node != NULL); | |
357 if (node->up == NULL) | |
358 { /* freeze the root subproblem */ | |
359 int k; | |
360 xassert(node->p == 1); | |
361 xassert(tree->root_m == 0); | |
362 xassert(tree->root_type == NULL); | |
363 xassert(tree->root_lb == NULL); | |
364 xassert(tree->root_ub == NULL); | |
365 xassert(tree->root_stat == NULL); | |
366 tree->root_m = m; | |
367 tree->root_type = xcalloc(1+m+n, sizeof(char)); | |
368 tree->root_lb = xcalloc(1+m+n, sizeof(double)); | |
369 tree->root_ub = xcalloc(1+m+n, sizeof(double)); | |
370 tree->root_stat = xcalloc(1+m+n, sizeof(char)); | |
371 for (k = 1; k <= m+n; k++) | |
372 { if (k <= m) | |
373 { GLPROW *row = mip->row[k]; | |
374 tree->root_type[k] = (char)row->type; | |
375 tree->root_lb[k] = row->lb; | |
376 tree->root_ub[k] = row->ub; | |
377 tree->root_stat[k] = (char)row->stat; | |
378 } | |
379 else | |
380 { GLPCOL *col = mip->col[k-m]; | |
381 tree->root_type[k] = (char)col->type; | |
382 tree->root_lb[k] = col->lb; | |
383 tree->root_ub[k] = col->ub; | |
384 tree->root_stat[k] = (char)col->stat; | |
385 } | |
386 } | |
387 } | |
388 else | |
389 { /* freeze non-root subproblem */ | |
390 int root_m = tree->root_m; | |
391 int pred_m = tree->pred_m; | |
392 int i, j, k; | |
393 xassert(pred_m <= m); | |
394 /* build change lists for rows and columns which exist in the | |
395 parent subproblem */ | |
396 xassert(node->b_ptr == NULL); | |
397 xassert(node->s_ptr == NULL); | |
398 for (k = 1; k <= pred_m + n; k++) | |
399 { int pred_type, pred_stat, type, stat; | |
400 double pred_lb, pred_ub, lb, ub; | |
401 /* determine attributes in the parent subproblem */ | |
402 pred_type = tree->pred_type[k]; | |
403 pred_lb = tree->pred_lb[k]; | |
404 pred_ub = tree->pred_ub[k]; | |
405 pred_stat = tree->pred_stat[k]; | |
406 /* determine attributes in the current subproblem */ | |
407 if (k <= pred_m) | |
408 { GLPROW *row = mip->row[k]; | |
409 type = row->type; | |
410 lb = row->lb; | |
411 ub = row->ub; | |
412 stat = row->stat; | |
413 } | |
414 else | |
415 { GLPCOL *col = mip->col[k - pred_m]; | |
416 type = col->type; | |
417 lb = col->lb; | |
418 ub = col->ub; | |
419 stat = col->stat; | |
420 } | |
421 /* save type and bounds of a row/column, if changed */ | |
422 if (!(pred_type == type && pred_lb == lb && pred_ub == ub)) | |
423 { IOSBND *b; | |
424 b = dmp_get_atom(tree->pool, sizeof(IOSBND)); | |
425 b->k = k; | |
426 b->type = (unsigned char)type; | |
427 b->lb = lb; | |
428 b->ub = ub; | |
429 b->next = node->b_ptr; | |
430 node->b_ptr = b; | |
431 } | |
432 /* save status of a row/column, if changed */ | |
433 if (pred_stat != stat) | |
434 { IOSTAT *s; | |
435 s = dmp_get_atom(tree->pool, sizeof(IOSTAT)); | |
436 s->k = k; | |
437 s->stat = (unsigned char)stat; | |
438 s->next = node->s_ptr; | |
439 node->s_ptr = s; | |
440 } | |
441 } | |
442 /* save new rows added to the current subproblem */ | |
443 xassert(node->r_ptr == NULL); | |
444 if (pred_m < m) | |
445 { int i, len, *ind; | |
446 double *val; | |
447 ind = xcalloc(1+n, sizeof(int)); | |
448 val = xcalloc(1+n, sizeof(double)); | |
449 for (i = m; i > pred_m; i--) | |
450 { GLPROW *row = mip->row[i]; | |
451 IOSROW *r; | |
452 const char *name; | |
453 r = dmp_get_atom(tree->pool, sizeof(IOSROW)); | |
454 name = glp_get_row_name(mip, i); | |
455 if (name == NULL) | |
456 r->name = NULL; | |
457 else | |
458 { r->name = dmp_get_atom(tree->pool, strlen(name)+1); | |
459 strcpy(r->name, name); | |
460 } | |
461 #if 1 /* 20/IX-2008 */ | |
462 r->origin = row->origin; | |
463 r->klass = row->klass; | |
464 #endif | |
465 r->type = (unsigned char)row->type; | |
466 r->lb = row->lb; | |
467 r->ub = row->ub; | |
468 r->ptr = NULL; | |
469 len = glp_get_mat_row(mip, i, ind, val); | |
470 for (k = 1; k <= len; k++) | |
471 { IOSAIJ *a; | |
472 a = dmp_get_atom(tree->pool, sizeof(IOSAIJ)); | |
473 a->j = ind[k]; | |
474 a->val = val[k]; | |
475 a->next = r->ptr; | |
476 r->ptr = a; | |
477 } | |
478 r->rii = row->rii; | |
479 r->stat = (unsigned char)row->stat; | |
480 r->next = node->r_ptr; | |
481 node->r_ptr = r; | |
482 } | |
483 xfree(ind); | |
484 xfree(val); | |
485 } | |
486 /* remove all rows missing in the root subproblem */ | |
487 if (m != root_m) | |
488 { int nrs, *num; | |
489 nrs = m - root_m; | |
490 xassert(nrs > 0); | |
491 num = xcalloc(1+nrs, sizeof(int)); | |
492 for (i = 1; i <= nrs; i++) num[i] = root_m + i; | |
493 glp_del_rows(mip, nrs, num); | |
494 xfree(num); | |
495 } | |
496 m = mip->m; | |
497 /* and restore attributes of all rows and columns for the root | |
498 subproblem */ | |
499 xassert(m == root_m); | |
500 for (i = 1; i <= m; i++) | |
501 { glp_set_row_bnds(mip, i, tree->root_type[i], | |
502 tree->root_lb[i], tree->root_ub[i]); | |
503 glp_set_row_stat(mip, i, tree->root_stat[i]); | |
504 } | |
505 for (j = 1; j <= n; j++) | |
506 { glp_set_col_bnds(mip, j, tree->root_type[m+j], | |
507 tree->root_lb[m+j], tree->root_ub[m+j]); | |
508 glp_set_col_stat(mip, j, tree->root_stat[m+j]); | |
509 } | |
510 #if 1 | |
511 /* remove all edges and cliques missing in the conflict graph | |
512 for the root subproblem */ | |
513 /* (not implemented yet) */ | |
514 #endif | |
515 } | |
516 /* the current subproblem has been frozen */ | |
517 tree->curr = NULL; | |
518 return; | |
519 } | |
520 | |
521 /*********************************************************************** | |
522 * NAME | |
523 * | |
524 * ios_clone_node - clone specified subproblem | |
525 * | |
526 * SYNOPSIS | |
527 * | |
528 * #include "glpios.h" | |
529 * void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]); | |
530 * | |
531 * DESCRIPTION | |
532 * | |
533 * The routine ios_clone_node clones the specified subproblem, whose | |
534 * reference number is p, creating its nnn exact copies. Note that the | |
535 * specified subproblem must be active and must be in the frozen state | |
536 * (i.e. it must not be the current subproblem). | |
537 * | |
538 * Each clone, an exact copy of the specified subproblem, becomes a new | |
539 * active subproblem added to the end of the active list. After cloning | |
540 * the specified subproblem becomes inactive. | |
541 * | |
542 * The reference numbers of clone subproblems are stored to locations | |
543 * ref[1], ..., ref[nnn]. */ | |
544 | |
545 static int get_slot(glp_tree *tree) | |
546 { int p; | |
547 /* if no free slots are available, increase the room */ | |
548 if (tree->avail == 0) | |
549 { int nslots = tree->nslots; | |
550 IOSLOT *save = tree->slot; | |
551 if (nslots == 0) | |
552 tree->nslots = 20; | |
553 else | |
554 { tree->nslots = nslots + nslots; | |
555 xassert(tree->nslots > nslots); | |
556 } | |
557 tree->slot = xcalloc(1+tree->nslots, sizeof(IOSLOT)); | |
558 if (save != NULL) | |
559 { memcpy(&tree->slot[1], &save[1], nslots * sizeof(IOSLOT)); | |
560 xfree(save); | |
561 } | |
562 /* push more free slots into the stack */ | |
563 for (p = tree->nslots; p > nslots; p--) | |
564 { tree->slot[p].node = NULL; | |
565 tree->slot[p].next = tree->avail; | |
566 tree->avail = p; | |
567 } | |
568 } | |
569 /* pull a free slot from the stack */ | |
570 p = tree->avail; | |
571 tree->avail = tree->slot[p].next; | |
572 xassert(tree->slot[p].node == NULL); | |
573 tree->slot[p].next = 0; | |
574 return p; | |
575 } | |
576 | |
577 static IOSNPD *new_node(glp_tree *tree, IOSNPD *parent) | |
578 { IOSNPD *node; | |
579 int p; | |
580 /* pull a free slot for the new node */ | |
581 p = get_slot(tree); | |
582 /* create descriptor of the new subproblem */ | |
583 node = dmp_get_atom(tree->pool, sizeof(IOSNPD)); | |
584 tree->slot[p].node = node; | |
585 node->p = p; | |
586 node->up = parent; | |
587 node->level = (parent == NULL ? 0 : parent->level + 1); | |
588 node->count = 0; | |
589 node->b_ptr = NULL; | |
590 node->s_ptr = NULL; | |
591 node->r_ptr = NULL; | |
592 node->solved = 0; | |
593 #if 0 | |
594 node->own_nn = node->own_nc = 0; | |
595 node->e_ptr = NULL; | |
596 #endif | |
597 #if 1 /* 04/X-2008 */ | |
598 node->lp_obj = (parent == NULL ? (tree->mip->dir == GLP_MIN ? | |
599 -DBL_MAX : +DBL_MAX) : parent->lp_obj); | |
600 #endif | |
601 node->bound = (parent == NULL ? (tree->mip->dir == GLP_MIN ? | |
602 -DBL_MAX : +DBL_MAX) : parent->bound); | |
603 node->br_var = 0; | |
604 node->br_val = 0.0; | |
605 node->ii_cnt = 0; | |
606 node->ii_sum = 0.0; | |
607 #if 1 /* 30/XI-2009 */ | |
608 node->changed = 0; | |
609 #endif | |
610 if (tree->parm->cb_size == 0) | |
611 node->data = NULL; | |
612 else | |
613 { node->data = dmp_get_atom(tree->pool, tree->parm->cb_size); | |
614 memset(node->data, 0, tree->parm->cb_size); | |
615 } | |
616 node->temp = NULL; | |
617 node->prev = tree->tail; | |
618 node->next = NULL; | |
619 /* add the new subproblem to the end of the active list */ | |
620 if (tree->head == NULL) | |
621 tree->head = node; | |
622 else | |
623 tree->tail->next = node; | |
624 tree->tail = node; | |
625 tree->a_cnt++; | |
626 tree->n_cnt++; | |
627 tree->t_cnt++; | |
628 /* increase the number of child subproblems */ | |
629 if (parent == NULL) | |
630 xassert(p == 1); | |
631 else | |
632 parent->count++; | |
633 return node; | |
634 } | |
635 | |
636 void ios_clone_node(glp_tree *tree, int p, int nnn, int ref[]) | |
637 { IOSNPD *node; | |
638 int k; | |
639 /* obtain pointer to the subproblem to be cloned */ | |
640 xassert(1 <= p && p <= tree->nslots); | |
641 node = tree->slot[p].node; | |
642 xassert(node != NULL); | |
643 /* the specified subproblem must be active */ | |
644 xassert(node->count == 0); | |
645 /* and must be in the frozen state */ | |
646 xassert(tree->curr != node); | |
647 /* remove the specified subproblem from the active list, because | |
648 it becomes inactive */ | |
649 if (node->prev == NULL) | |
650 tree->head = node->next; | |
651 else | |
652 node->prev->next = node->next; | |
653 if (node->next == NULL) | |
654 tree->tail = node->prev; | |
655 else | |
656 node->next->prev = node->prev; | |
657 node->prev = node->next = NULL; | |
658 tree->a_cnt--; | |
659 /* create clone subproblems */ | |
660 xassert(nnn > 0); | |
661 for (k = 1; k <= nnn; k++) | |
662 ref[k] = new_node(tree, node)->p; | |
663 return; | |
664 } | |
665 | |
666 /*********************************************************************** | |
667 * NAME | |
668 * | |
669 * ios_delete_node - delete specified subproblem | |
670 * | |
671 * SYNOPSIS | |
672 * | |
673 * #include "glpios.h" | |
674 * void ios_delete_node(glp_tree *tree, int p); | |
675 * | |
676 * DESCRIPTION | |
677 * | |
678 * The routine ios_delete_node deletes the specified subproblem, whose | |
679 * reference number is p. The subproblem must be active and must be in | |
680 * the frozen state (i.e. it must not be the current subproblem). | |
681 * | |
682 * Note that deletion is performed recursively, i.e. if a subproblem to | |
683 * be deleted is the only child of its parent, the parent subproblem is | |
684 * also deleted, etc. */ | |
685 | |
686 void ios_delete_node(glp_tree *tree, int p) | |
687 { IOSNPD *node, *temp; | |
688 /* obtain pointer to the subproblem to be deleted */ | |
689 xassert(1 <= p && p <= tree->nslots); | |
690 node = tree->slot[p].node; | |
691 xassert(node != NULL); | |
692 /* the specified subproblem must be active */ | |
693 xassert(node->count == 0); | |
694 /* and must be in the frozen state */ | |
695 xassert(tree->curr != node); | |
696 /* remove the specified subproblem from the active list, because | |
697 it is gone from the tree */ | |
698 if (node->prev == NULL) | |
699 tree->head = node->next; | |
700 else | |
701 node->prev->next = node->next; | |
702 if (node->next == NULL) | |
703 tree->tail = node->prev; | |
704 else | |
705 node->next->prev = node->prev; | |
706 node->prev = node->next = NULL; | |
707 tree->a_cnt--; | |
708 loop: /* recursive deletion starts here */ | |
709 /* delete the bound change list */ | |
710 { IOSBND *b; | |
711 while (node->b_ptr != NULL) | |
712 { b = node->b_ptr; | |
713 node->b_ptr = b->next; | |
714 dmp_free_atom(tree->pool, b, sizeof(IOSBND)); | |
715 } | |
716 } | |
717 /* delete the status change list */ | |
718 { IOSTAT *s; | |
719 while (node->s_ptr != NULL) | |
720 { s = node->s_ptr; | |
721 node->s_ptr = s->next; | |
722 dmp_free_atom(tree->pool, s, sizeof(IOSTAT)); | |
723 } | |
724 } | |
725 /* delete the row addition list */ | |
726 while (node->r_ptr != NULL) | |
727 { IOSROW *r; | |
728 r = node->r_ptr; | |
729 if (r->name != NULL) | |
730 dmp_free_atom(tree->pool, r->name, strlen(r->name)+1); | |
731 while (r->ptr != NULL) | |
732 { IOSAIJ *a; | |
733 a = r->ptr; | |
734 r->ptr = a->next; | |
735 dmp_free_atom(tree->pool, a, sizeof(IOSAIJ)); | |
736 } | |
737 node->r_ptr = r->next; | |
738 dmp_free_atom(tree->pool, r, sizeof(IOSROW)); | |
739 } | |
740 #if 0 | |
741 /* delete the edge addition list */ | |
742 /* delete the clique addition list */ | |
743 /* (not implemented yet) */ | |
744 xassert(node->own_nn == 0); | |
745 xassert(node->own_nc == 0); | |
746 xassert(node->e_ptr == NULL); | |
747 #endif | |
748 /* free application-specific data */ | |
749 if (tree->parm->cb_size == 0) | |
750 xassert(node->data == NULL); | |
751 else | |
752 dmp_free_atom(tree->pool, node->data, tree->parm->cb_size); | |
753 /* free the corresponding node slot */ | |
754 p = node->p; | |
755 xassert(tree->slot[p].node == node); | |
756 tree->slot[p].node = NULL; | |
757 tree->slot[p].next = tree->avail; | |
758 tree->avail = p; | |
759 /* save pointer to the parent subproblem */ | |
760 temp = node->up; | |
761 /* delete the subproblem descriptor */ | |
762 dmp_free_atom(tree->pool, node, sizeof(IOSNPD)); | |
763 tree->n_cnt--; | |
764 /* take pointer to the parent subproblem */ | |
765 node = temp; | |
766 if (node != NULL) | |
767 { /* the parent subproblem exists; decrease the number of its | |
768 child subproblems */ | |
769 xassert(node->count > 0); | |
770 node->count--; | |
771 /* if now the parent subproblem has no childs, it also must be | |
772 deleted */ | |
773 if (node->count == 0) goto loop; | |
774 } | |
775 return; | |
776 } | |
777 | |
778 /*********************************************************************** | |
779 * NAME | |
780 * | |
781 * ios_delete_tree - delete branch-and-bound tree | |
782 * | |
783 * SYNOPSIS | |
784 * | |
785 * #include "glpios.h" | |
786 * void ios_delete_tree(glp_tree *tree); | |
787 * | |
788 * DESCRIPTION | |
789 * | |
790 * The routine ios_delete_tree deletes the branch-and-bound tree, which | |
791 * the parameter tree points to, and frees all the memory allocated to | |
792 * this program object. | |
793 * | |
794 * On exit components of the problem object are restored to correspond | |
795 * to the original MIP passed to the routine ios_create_tree. */ | |
796 | |
797 void ios_delete_tree(glp_tree *tree) | |
798 { glp_prob *mip = tree->mip; | |
799 int i, j; | |
800 int m = mip->m; | |
801 int n = mip->n; | |
802 xassert(mip->tree == tree); | |
803 /* remove all additional rows */ | |
804 if (m != tree->orig_m) | |
805 { int nrs, *num; | |
806 nrs = m - tree->orig_m; | |
807 xassert(nrs > 0); | |
808 num = xcalloc(1+nrs, sizeof(int)); | |
809 for (i = 1; i <= nrs; i++) num[i] = tree->orig_m + i; | |
810 glp_del_rows(mip, nrs, num); | |
811 xfree(num); | |
812 } | |
813 m = tree->orig_m; | |
814 /* restore original attributes of rows and columns */ | |
815 xassert(m == tree->orig_m); | |
816 xassert(n == tree->n); | |
817 for (i = 1; i <= m; i++) | |
818 { glp_set_row_bnds(mip, i, tree->orig_type[i], | |
819 tree->orig_lb[i], tree->orig_ub[i]); | |
820 glp_set_row_stat(mip, i, tree->orig_stat[i]); | |
821 mip->row[i]->prim = tree->orig_prim[i]; | |
822 mip->row[i]->dual = tree->orig_dual[i]; | |
823 } | |
824 for (j = 1; j <= n; j++) | |
825 { glp_set_col_bnds(mip, j, tree->orig_type[m+j], | |
826 tree->orig_lb[m+j], tree->orig_ub[m+j]); | |
827 glp_set_col_stat(mip, j, tree->orig_stat[m+j]); | |
828 mip->col[j]->prim = tree->orig_prim[m+j]; | |
829 mip->col[j]->dual = tree->orig_dual[m+j]; | |
830 } | |
831 mip->pbs_stat = mip->dbs_stat = GLP_FEAS; | |
832 mip->obj_val = tree->orig_obj; | |
833 /* delete the branch-and-bound tree */ | |
834 xassert(tree->local != NULL); | |
835 ios_delete_pool(tree, tree->local); | |
836 dmp_delete_pool(tree->pool); | |
837 xfree(tree->orig_type); | |
838 xfree(tree->orig_lb); | |
839 xfree(tree->orig_ub); | |
840 xfree(tree->orig_stat); | |
841 xfree(tree->orig_prim); | |
842 xfree(tree->orig_dual); | |
843 xfree(tree->slot); | |
844 if (tree->root_type != NULL) xfree(tree->root_type); | |
845 if (tree->root_lb != NULL) xfree(tree->root_lb); | |
846 if (tree->root_ub != NULL) xfree(tree->root_ub); | |
847 if (tree->root_stat != NULL) xfree(tree->root_stat); | |
848 xfree(tree->non_int); | |
849 #if 0 | |
850 xfree(tree->n_ref); | |
851 xfree(tree->c_ref); | |
852 xfree(tree->j_ref); | |
853 #endif | |
854 if (tree->pcost != NULL) ios_pcost_free(tree); | |
855 xfree(tree->iwrk); | |
856 xfree(tree->dwrk); | |
857 #if 0 | |
858 scg_delete_graph(tree->g); | |
859 #endif | |
860 if (tree->pred_type != NULL) xfree(tree->pred_type); | |
861 if (tree->pred_lb != NULL) xfree(tree->pred_lb); | |
862 if (tree->pred_ub != NULL) xfree(tree->pred_ub); | |
863 if (tree->pred_stat != NULL) xfree(tree->pred_stat); | |
864 #if 0 | |
865 xassert(tree->cut_gen == NULL); | |
866 #endif | |
867 xassert(tree->mir_gen == NULL); | |
868 xassert(tree->clq_gen == NULL); | |
869 xfree(tree); | |
870 mip->tree = NULL; | |
871 return; | |
872 } | |
873 | |
874 /*********************************************************************** | |
875 * NAME | |
876 * | |
877 * ios_eval_degrad - estimate obj. degrad. for down- and up-branches | |
878 * | |
879 * SYNOPSIS | |
880 * | |
881 * #include "glpios.h" | |
882 * void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up); | |
883 * | |
884 * DESCRIPTION | |
885 * | |
886 * Given optimal basis to LP relaxation of the current subproblem the | |
887 * routine ios_eval_degrad performs the dual ratio test to compute the | |
888 * objective values in the adjacent basis for down- and up-branches, | |
889 * which are stored in locations *dn and *up, assuming that x[j] is a | |
890 * variable chosen to branch upon. */ | |
891 | |
892 void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up) | |
893 { glp_prob *mip = tree->mip; | |
894 int m = mip->m, n = mip->n; | |
895 int len, kase, k, t, stat; | |
896 double alfa, beta, gamma, delta, dz; | |
897 int *ind = tree->iwrk; | |
898 double *val = tree->dwrk; | |
899 /* current basis must be optimal */ | |
900 xassert(glp_get_status(mip) == GLP_OPT); | |
901 /* basis factorization must exist */ | |
902 xassert(glp_bf_exists(mip)); | |
903 /* obtain (fractional) value of x[j] in optimal basic solution | |
904 to LP relaxation of the current subproblem */ | |
905 xassert(1 <= j && j <= n); | |
906 beta = mip->col[j]->prim; | |
907 /* since the value of x[j] is fractional, it is basic; compute | |
908 corresponding row of the simplex table */ | |
909 len = lpx_eval_tab_row(mip, m+j, ind, val); | |
910 /* kase < 0 means down-branch; kase > 0 means up-branch */ | |
911 for (kase = -1; kase <= +1; kase += 2) | |
912 { /* for down-branch we introduce new upper bound floor(beta) | |
913 for x[j]; similarly, for up-branch we introduce new lower | |
914 bound ceil(beta) for x[j]; in the current basis this new | |
915 upper/lower bound is violated, so in the adjacent basis | |
916 x[j] will leave the basis and go to its new upper/lower | |
917 bound; we need to know which non-basic variable x[k] should | |
918 enter the basis to keep dual feasibility */ | |
919 #if 0 /* 23/XI-2009 */ | |
920 k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-7); | |
921 #else | |
922 k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-9); | |
923 #endif | |
924 /* if no variable has been chosen, current basis being primal | |
925 infeasible due to the new upper/lower bound of x[j] is dual | |
926 unbounded, therefore, LP relaxation to corresponding branch | |
927 has no primal feasible solution */ | |
928 if (k == 0) | |
929 { if (mip->dir == GLP_MIN) | |
930 { if (kase < 0) | |
931 *dn = +DBL_MAX; | |
932 else | |
933 *up = +DBL_MAX; | |
934 } | |
935 else if (mip->dir == GLP_MAX) | |
936 { if (kase < 0) | |
937 *dn = -DBL_MAX; | |
938 else | |
939 *up = -DBL_MAX; | |
940 } | |
941 else | |
942 xassert(mip != mip); | |
943 continue; | |
944 } | |
945 xassert(1 <= k && k <= m+n); | |
946 /* row of the simplex table corresponding to specified basic | |
947 variable x[j] is the following: | |
948 x[j] = ... + alfa * x[k] + ... ; | |
949 we need to know influence coefficient, alfa, at non-basic | |
950 variable x[k] chosen with the dual ratio test */ | |
951 for (t = 1; t <= len; t++) | |
952 if (ind[t] == k) break; | |
953 xassert(1 <= t && t <= len); | |
954 alfa = val[t]; | |
955 /* determine status and reduced cost of variable x[k] */ | |
956 if (k <= m) | |
957 { stat = mip->row[k]->stat; | |
958 gamma = mip->row[k]->dual; | |
959 } | |
960 else | |
961 { stat = mip->col[k-m]->stat; | |
962 gamma = mip->col[k-m]->dual; | |
963 } | |
964 /* x[k] cannot be basic or fixed non-basic */ | |
965 xassert(stat == GLP_NL || stat == GLP_NU || stat == GLP_NF); | |
966 /* if the current basis is dual degenerative, some reduced | |
967 costs, which are close to zero, may have wrong sign due to | |
968 round-off errors, so correct the sign of gamma */ | |
969 if (mip->dir == GLP_MIN) | |
970 { if (stat == GLP_NL && gamma < 0.0 || | |
971 stat == GLP_NU && gamma > 0.0 || | |
972 stat == GLP_NF) gamma = 0.0; | |
973 } | |
974 else if (mip->dir == GLP_MAX) | |
975 { if (stat == GLP_NL && gamma > 0.0 || | |
976 stat == GLP_NU && gamma < 0.0 || | |
977 stat == GLP_NF) gamma = 0.0; | |
978 } | |
979 else | |
980 xassert(mip != mip); | |
981 /* determine the change of x[j] in the adjacent basis: | |
982 delta x[j] = new x[j] - old x[j] */ | |
983 delta = (kase < 0 ? floor(beta) : ceil(beta)) - beta; | |
984 /* compute the change of x[k] in the adjacent basis: | |
985 delta x[k] = new x[k] - old x[k] = delta x[j] / alfa */ | |
986 delta /= alfa; | |
987 /* compute the change of the objective in the adjacent basis: | |
988 delta z = new z - old z = gamma * delta x[k] */ | |
989 dz = gamma * delta; | |
990 if (mip->dir == GLP_MIN) | |
991 xassert(dz >= 0.0); | |
992 else if (mip->dir == GLP_MAX) | |
993 xassert(dz <= 0.0); | |
994 else | |
995 xassert(mip != mip); | |
996 /* compute the new objective value in the adjacent basis: | |
997 new z = old z + delta z */ | |
998 if (kase < 0) | |
999 *dn = mip->obj_val + dz; | |
1000 else | |
1001 *up = mip->obj_val + dz; | |
1002 } | |
1003 /*xprintf("obj = %g; dn = %g; up = %g\n", | |
1004 mip->obj_val, *dn, *up);*/ | |
1005 return; | |
1006 } | |
1007 | |
1008 /*********************************************************************** | |
1009 * NAME | |
1010 * | |
1011 * ios_round_bound - improve local bound by rounding | |
1012 * | |
1013 * SYNOPSIS | |
1014 * | |
1015 * #include "glpios.h" | |
1016 * double ios_round_bound(glp_tree *tree, double bound); | |
1017 * | |
1018 * RETURNS | |
1019 * | |
1020 * For the given local bound for any integer feasible solution to the | |
1021 * current subproblem the routine ios_round_bound returns an improved | |
1022 * local bound for the same integer feasible solution. | |
1023 * | |
1024 * BACKGROUND | |
1025 * | |
1026 * Let the current subproblem has the following objective function: | |
1027 * | |
1028 * z = sum c[j] * x[j] + s >= b, (1) | |
1029 * j in J | |
1030 * | |
1031 * where J = {j: c[j] is non-zero and integer, x[j] is integer}, s is | |
1032 * the sum of terms corresponding to fixed variables, b is an initial | |
1033 * local bound (minimization). | |
1034 * | |
1035 * From (1) it follows that: | |
1036 * | |
1037 * d * sum (c[j] / d) * x[j] + s >= b, (2) | |
1038 * j in J | |
1039 * | |
1040 * or, equivalently, | |
1041 * | |
1042 * sum (c[j] / d) * x[j] >= (b - s) / d = h, (3) | |
1043 * j in J | |
1044 * | |
1045 * where d = gcd(c[j]). Since the left-hand side of (3) is integer, | |
1046 * h = (b - s) / d can be rounded up to the nearest integer: | |
1047 * | |
1048 * h' = ceil(h) = (b' - s) / d, (4) | |
1049 * | |
1050 * that gives an rounded, improved local bound: | |
1051 * | |
1052 * b' = d * h' + s. (5) | |
1053 * | |
1054 * In case of maximization '>=' in (1) should be replaced by '<=' that | |
1055 * leads to the following formula: | |
1056 * | |
1057 * h' = floor(h) = (b' - s) / d, (6) | |
1058 * | |
1059 * which should used in the same way as (4). | |
1060 * | |
1061 * NOTE: If b is a valid local bound for a child of the current | |
1062 * subproblem, b' is also valid for that child subproblem. */ | |
1063 | |
1064 double ios_round_bound(glp_tree *tree, double bound) | |
1065 { glp_prob *mip = tree->mip; | |
1066 int n = mip->n; | |
1067 int d, j, nn, *c = tree->iwrk; | |
1068 double s, h; | |
1069 /* determine c[j] and compute s */ | |
1070 nn = 0, s = mip->c0, d = 0; | |
1071 for (j = 1; j <= n; j++) | |
1072 { GLPCOL *col = mip->col[j]; | |
1073 if (col->coef == 0.0) continue; | |
1074 if (col->type == GLP_FX) | |
1075 { /* fixed variable */ | |
1076 s += col->coef * col->prim; | |
1077 } | |
1078 else | |
1079 { /* non-fixed variable */ | |
1080 if (col->kind != GLP_IV) goto skip; | |
1081 if (col->coef != floor(col->coef)) goto skip; | |
1082 if (fabs(col->coef) <= (double)INT_MAX) | |
1083 c[++nn] = (int)fabs(col->coef); | |
1084 else | |
1085 d = 1; | |
1086 } | |
1087 } | |
1088 /* compute d = gcd(c[1],...c[nn]) */ | |
1089 if (d == 0) | |
1090 { if (nn == 0) goto skip; | |
1091 d = gcdn(nn, c); | |
1092 } | |
1093 xassert(d > 0); | |
1094 /* compute new local bound */ | |
1095 if (mip->dir == GLP_MIN) | |
1096 { if (bound != +DBL_MAX) | |
1097 { h = (bound - s) / (double)d; | |
1098 if (h >= floor(h) + 0.001) | |
1099 { /* round up */ | |
1100 h = ceil(h); | |
1101 /*xprintf("d = %d; old = %g; ", d, bound);*/ | |
1102 bound = (double)d * h + s; | |
1103 /*xprintf("new = %g\n", bound);*/ | |
1104 } | |
1105 } | |
1106 } | |
1107 else if (mip->dir == GLP_MAX) | |
1108 { if (bound != -DBL_MAX) | |
1109 { h = (bound - s) / (double)d; | |
1110 if (h <= ceil(h) - 0.001) | |
1111 { /* round down */ | |
1112 h = floor(h); | |
1113 bound = (double)d * h + s; | |
1114 } | |
1115 } | |
1116 } | |
1117 else | |
1118 xassert(mip != mip); | |
1119 skip: return bound; | |
1120 } | |
1121 | |
1122 /*********************************************************************** | |
1123 * NAME | |
1124 * | |
1125 * ios_is_hopeful - check if subproblem is hopeful | |
1126 * | |
1127 * SYNOPSIS | |
1128 * | |
1129 * #include "glpios.h" | |
1130 * int ios_is_hopeful(glp_tree *tree, double bound); | |
1131 * | |
1132 * DESCRIPTION | |
1133 * | |
1134 * Given the local bound of a subproblem the routine ios_is_hopeful | |
1135 * checks if the subproblem can have an integer optimal solution which | |
1136 * is better than the best one currently known. | |
1137 * | |
1138 * RETURNS | |
1139 * | |
1140 * If the subproblem can have a better integer optimal solution, the | |
1141 * routine returns non-zero; otherwise, if the corresponding branch can | |
1142 * be pruned, the routine returns zero. */ | |
1143 | |
1144 int ios_is_hopeful(glp_tree *tree, double bound) | |
1145 { glp_prob *mip = tree->mip; | |
1146 int ret = 1; | |
1147 double eps; | |
1148 if (mip->mip_stat == GLP_FEAS) | |
1149 { eps = tree->parm->tol_obj * (1.0 + fabs(mip->mip_obj)); | |
1150 switch (mip->dir) | |
1151 { case GLP_MIN: | |
1152 if (bound >= mip->mip_obj - eps) ret = 0; | |
1153 break; | |
1154 case GLP_MAX: | |
1155 if (bound <= mip->mip_obj + eps) ret = 0; | |
1156 break; | |
1157 default: | |
1158 xassert(mip != mip); | |
1159 } | |
1160 } | |
1161 else | |
1162 { switch (mip->dir) | |
1163 { case GLP_MIN: | |
1164 if (bound == +DBL_MAX) ret = 0; | |
1165 break; | |
1166 case GLP_MAX: | |
1167 if (bound == -DBL_MAX) ret = 0; | |
1168 break; | |
1169 default: | |
1170 xassert(mip != mip); | |
1171 } | |
1172 } | |
1173 return ret; | |
1174 } | |
1175 | |
1176 /*********************************************************************** | |
1177 * NAME | |
1178 * | |
1179 * ios_best_node - find active node with best local bound | |
1180 * | |
1181 * SYNOPSIS | |
1182 * | |
1183 * #include "glpios.h" | |
1184 * int ios_best_node(glp_tree *tree); | |
1185 * | |
1186 * DESCRIPTION | |
1187 * | |
1188 * The routine ios_best_node finds an active node whose local bound is | |
1189 * best among other active nodes. | |
1190 * | |
1191 * It is understood that the integer optimal solution of the original | |
1192 * mip problem cannot be better than the best bound, so the best bound | |
1193 * is an lower (minimization) or upper (maximization) global bound for | |
1194 * the original problem. | |
1195 * | |
1196 * RETURNS | |
1197 * | |
1198 * The routine ios_best_node returns the subproblem reference number | |
1199 * for the best node. However, if the tree is empty, it returns zero. */ | |
1200 | |
1201 int ios_best_node(glp_tree *tree) | |
1202 { IOSNPD *node, *best = NULL; | |
1203 switch (tree->mip->dir) | |
1204 { case GLP_MIN: | |
1205 /* minimization */ | |
1206 for (node = tree->head; node != NULL; node = node->next) | |
1207 if (best == NULL || best->bound > node->bound) | |
1208 best = node; | |
1209 break; | |
1210 case GLP_MAX: | |
1211 /* maximization */ | |
1212 for (node = tree->head; node != NULL; node = node->next) | |
1213 if (best == NULL || best->bound < node->bound) | |
1214 best = node; | |
1215 break; | |
1216 default: | |
1217 xassert(tree != tree); | |
1218 } | |
1219 return best == NULL ? 0 : best->p; | |
1220 } | |
1221 | |
1222 /*********************************************************************** | |
1223 * NAME | |
1224 * | |
1225 * ios_relative_gap - compute relative mip gap | |
1226 * | |
1227 * SYNOPSIS | |
1228 * | |
1229 * #include "glpios.h" | |
1230 * double ios_relative_gap(glp_tree *tree); | |
1231 * | |
1232 * DESCRIPTION | |
1233 * | |
1234 * The routine ios_relative_gap computes the relative mip gap using the | |
1235 * formula: | |
1236 * | |
1237 * gap = |best_mip - best_bnd| / (|best_mip| + DBL_EPSILON), | |
1238 * | |
1239 * where best_mip is the best integer feasible solution found so far, | |
1240 * best_bnd is the best (global) bound. If no integer feasible solution | |
1241 * has been found yet, rel_gap is set to DBL_MAX. | |
1242 * | |
1243 * RETURNS | |
1244 * | |
1245 * The routine ios_relative_gap returns the relative mip gap. */ | |
1246 | |
1247 double ios_relative_gap(glp_tree *tree) | |
1248 { glp_prob *mip = tree->mip; | |
1249 int p; | |
1250 double best_mip, best_bnd, gap; | |
1251 if (mip->mip_stat == GLP_FEAS) | |
1252 { best_mip = mip->mip_obj; | |
1253 p = ios_best_node(tree); | |
1254 if (p == 0) | |
1255 { /* the tree is empty */ | |
1256 gap = 0.0; | |
1257 } | |
1258 else | |
1259 { best_bnd = tree->slot[p].node->bound; | |
1260 gap = fabs(best_mip - best_bnd) / (fabs(best_mip) + | |
1261 DBL_EPSILON); | |
1262 } | |
1263 } | |
1264 else | |
1265 { /* no integer feasible solution has been found yet */ | |
1266 gap = DBL_MAX; | |
1267 } | |
1268 return gap; | |
1269 } | |
1270 | |
1271 /*********************************************************************** | |
1272 * NAME | |
1273 * | |
1274 * ios_solve_node - solve LP relaxation of current subproblem | |
1275 * | |
1276 * SYNOPSIS | |
1277 * | |
1278 * #include "glpios.h" | |
1279 * int ios_solve_node(glp_tree *tree); | |
1280 * | |
1281 * DESCRIPTION | |
1282 * | |
1283 * The routine ios_solve_node re-optimizes LP relaxation of the current | |
1284 * subproblem using the dual simplex method. | |
1285 * | |
1286 * RETURNS | |
1287 * | |
1288 * The routine returns the code which is reported by glp_simplex. */ | |
1289 | |
1290 int ios_solve_node(glp_tree *tree) | |
1291 { glp_prob *mip = tree->mip; | |
1292 glp_smcp parm; | |
1293 int ret; | |
1294 /* the current subproblem must exist */ | |
1295 xassert(tree->curr != NULL); | |
1296 /* set some control parameters */ | |
1297 glp_init_smcp(&parm); | |
1298 switch (tree->parm->msg_lev) | |
1299 { case GLP_MSG_OFF: | |
1300 parm.msg_lev = GLP_MSG_OFF; break; | |
1301 case GLP_MSG_ERR: | |
1302 parm.msg_lev = GLP_MSG_ERR; break; | |
1303 case GLP_MSG_ON: | |
1304 case GLP_MSG_ALL: | |
1305 parm.msg_lev = GLP_MSG_ON; break; | |
1306 case GLP_MSG_DBG: | |
1307 parm.msg_lev = GLP_MSG_ALL; break; | |
1308 default: | |
1309 xassert(tree != tree); | |
1310 } | |
1311 parm.meth = GLP_DUALP; | |
1312 if (tree->parm->msg_lev < GLP_MSG_DBG) | |
1313 parm.out_dly = tree->parm->out_dly; | |
1314 else | |
1315 parm.out_dly = 0; | |
1316 /* if the incumbent objective value is already known, use it to | |
1317 prematurely terminate the dual simplex search */ | |
1318 if (mip->mip_stat == GLP_FEAS) | |
1319 { switch (tree->mip->dir) | |
1320 { case GLP_MIN: | |
1321 parm.obj_ul = mip->mip_obj; | |
1322 break; | |
1323 case GLP_MAX: | |
1324 parm.obj_ll = mip->mip_obj; | |
1325 break; | |
1326 default: | |
1327 xassert(mip != mip); | |
1328 } | |
1329 } | |
1330 /* try to solve/re-optimize the LP relaxation */ | |
1331 ret = glp_simplex(mip, &parm); | |
1332 tree->curr->solved++; | |
1333 #if 0 | |
1334 xprintf("ret = %d; status = %d; pbs = %d; dbs = %d; some = %d\n", | |
1335 ret, glp_get_status(mip), mip->pbs_stat, mip->dbs_stat, | |
1336 mip->some); | |
1337 lpx_print_sol(mip, "sol"); | |
1338 #endif | |
1339 return ret; | |
1340 } | |
1341 | |
1342 /**********************************************************************/ | |
1343 | |
1344 IOSPOOL *ios_create_pool(glp_tree *tree) | |
1345 { /* create cut pool */ | |
1346 IOSPOOL *pool; | |
1347 #if 0 | |
1348 pool = dmp_get_atom(tree->pool, sizeof(IOSPOOL)); | |
1349 #else | |
1350 xassert(tree == tree); | |
1351 pool = xmalloc(sizeof(IOSPOOL)); | |
1352 #endif | |
1353 pool->size = 0; | |
1354 pool->head = pool->tail = NULL; | |
1355 pool->ord = 0, pool->curr = NULL; | |
1356 return pool; | |
1357 } | |
1358 | |
1359 int ios_add_row(glp_tree *tree, IOSPOOL *pool, | |
1360 const char *name, int klass, int flags, int len, const int ind[], | |
1361 const double val[], int type, double rhs) | |
1362 { /* add row (constraint) to the cut pool */ | |
1363 IOSCUT *cut; | |
1364 IOSAIJ *aij; | |
1365 int k; | |
1366 xassert(pool != NULL); | |
1367 cut = dmp_get_atom(tree->pool, sizeof(IOSCUT)); | |
1368 if (name == NULL || name[0] == '\0') | |
1369 cut->name = NULL; | |
1370 else | |
1371 { for (k = 0; name[k] != '\0'; k++) | |
1372 { if (k == 256) | |
1373 xerror("glp_ios_add_row: cut name too long\n"); | |
1374 if (iscntrl((unsigned char)name[k])) | |
1375 xerror("glp_ios_add_row: cut name contains invalid chara" | |
1376 "cter(s)\n"); | |
1377 } | |
1378 cut->name = dmp_get_atom(tree->pool, strlen(name)+1); | |
1379 strcpy(cut->name, name); | |
1380 } | |
1381 if (!(0 <= klass && klass <= 255)) | |
1382 xerror("glp_ios_add_row: klass = %d; invalid cut class\n", | |
1383 klass); | |
1384 cut->klass = (unsigned char)klass; | |
1385 if (flags != 0) | |
1386 xerror("glp_ios_add_row: flags = %d; invalid cut flags\n", | |
1387 flags); | |
1388 cut->ptr = NULL; | |
1389 if (!(0 <= len && len <= tree->n)) | |
1390 xerror("glp_ios_add_row: len = %d; invalid cut length\n", | |
1391 len); | |
1392 for (k = 1; k <= len; k++) | |
1393 { aij = dmp_get_atom(tree->pool, sizeof(IOSAIJ)); | |
1394 if (!(1 <= ind[k] && ind[k] <= tree->n)) | |
1395 xerror("glp_ios_add_row: ind[%d] = %d; column index out of " | |
1396 "range\n", k, ind[k]); | |
1397 aij->j = ind[k]; | |
1398 aij->val = val[k]; | |
1399 aij->next = cut->ptr; | |
1400 cut->ptr = aij; | |
1401 } | |
1402 if (!(type == GLP_LO || type == GLP_UP || type == GLP_FX)) | |
1403 xerror("glp_ios_add_row: type = %d; invalid cut type\n", | |
1404 type); | |
1405 cut->type = (unsigned char)type; | |
1406 cut->rhs = rhs; | |
1407 cut->prev = pool->tail; | |
1408 cut->next = NULL; | |
1409 if (cut->prev == NULL) | |
1410 pool->head = cut; | |
1411 else | |
1412 cut->prev->next = cut; | |
1413 pool->tail = cut; | |
1414 pool->size++; | |
1415 return pool->size; | |
1416 } | |
1417 | |
1418 IOSCUT *ios_find_row(IOSPOOL *pool, int i) | |
1419 { /* find row (constraint) in the cut pool */ | |
1420 /* (smart linear search) */ | |
1421 xassert(pool != NULL); | |
1422 xassert(1 <= i && i <= pool->size); | |
1423 if (pool->ord == 0) | |
1424 { xassert(pool->curr == NULL); | |
1425 pool->ord = 1; | |
1426 pool->curr = pool->head; | |
1427 } | |
1428 xassert(pool->curr != NULL); | |
1429 if (i < pool->ord) | |
1430 { if (i < pool->ord - i) | |
1431 { pool->ord = 1; | |
1432 pool->curr = pool->head; | |
1433 while (pool->ord != i) | |
1434 { pool->ord++; | |
1435 xassert(pool->curr != NULL); | |
1436 pool->curr = pool->curr->next; | |
1437 } | |
1438 } | |
1439 else | |
1440 { while (pool->ord != i) | |
1441 { pool->ord--; | |
1442 xassert(pool->curr != NULL); | |
1443 pool->curr = pool->curr->prev; | |
1444 } | |
1445 } | |
1446 } | |
1447 else if (i > pool->ord) | |
1448 { if (i - pool->ord < pool->size - i) | |
1449 { while (pool->ord != i) | |
1450 { pool->ord++; | |
1451 xassert(pool->curr != NULL); | |
1452 pool->curr = pool->curr->next; | |
1453 } | |
1454 } | |
1455 else | |
1456 { pool->ord = pool->size; | |
1457 pool->curr = pool->tail; | |
1458 while (pool->ord != i) | |
1459 { pool->ord--; | |
1460 xassert(pool->curr != NULL); | |
1461 pool->curr = pool->curr->prev; | |
1462 } | |
1463 } | |
1464 } | |
1465 xassert(pool->ord == i); | |
1466 xassert(pool->curr != NULL); | |
1467 return pool->curr; | |
1468 } | |
1469 | |
1470 void ios_del_row(glp_tree *tree, IOSPOOL *pool, int i) | |
1471 { /* remove row (constraint) from the cut pool */ | |
1472 IOSCUT *cut; | |
1473 IOSAIJ *aij; | |
1474 xassert(pool != NULL); | |
1475 if (!(1 <= i && i <= pool->size)) | |
1476 xerror("glp_ios_del_row: i = %d; cut number out of range\n", | |
1477 i); | |
1478 cut = ios_find_row(pool, i); | |
1479 xassert(pool->curr == cut); | |
1480 if (cut->next != NULL) | |
1481 pool->curr = cut->next; | |
1482 else if (cut->prev != NULL) | |
1483 pool->ord--, pool->curr = cut->prev; | |
1484 else | |
1485 pool->ord = 0, pool->curr = NULL; | |
1486 if (cut->name != NULL) | |
1487 dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1); | |
1488 if (cut->prev == NULL) | |
1489 { xassert(pool->head == cut); | |
1490 pool->head = cut->next; | |
1491 } | |
1492 else | |
1493 { xassert(cut->prev->next == cut); | |
1494 cut->prev->next = cut->next; | |
1495 } | |
1496 if (cut->next == NULL) | |
1497 { xassert(pool->tail == cut); | |
1498 pool->tail = cut->prev; | |
1499 } | |
1500 else | |
1501 { xassert(cut->next->prev == cut); | |
1502 cut->next->prev = cut->prev; | |
1503 } | |
1504 while (cut->ptr != NULL) | |
1505 { aij = cut->ptr; | |
1506 cut->ptr = aij->next; | |
1507 dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ)); | |
1508 } | |
1509 dmp_free_atom(tree->pool, cut, sizeof(IOSCUT)); | |
1510 pool->size--; | |
1511 return; | |
1512 } | |
1513 | |
1514 void ios_clear_pool(glp_tree *tree, IOSPOOL *pool) | |
1515 { /* remove all rows (constraints) from the cut pool */ | |
1516 xassert(pool != NULL); | |
1517 while (pool->head != NULL) | |
1518 { IOSCUT *cut = pool->head; | |
1519 pool->head = cut->next; | |
1520 if (cut->name != NULL) | |
1521 dmp_free_atom(tree->pool, cut->name, strlen(cut->name)+1); | |
1522 while (cut->ptr != NULL) | |
1523 { IOSAIJ *aij = cut->ptr; | |
1524 cut->ptr = aij->next; | |
1525 dmp_free_atom(tree->pool, aij, sizeof(IOSAIJ)); | |
1526 } | |
1527 dmp_free_atom(tree->pool, cut, sizeof(IOSCUT)); | |
1528 } | |
1529 pool->size = 0; | |
1530 pool->head = pool->tail = NULL; | |
1531 pool->ord = 0, pool->curr = NULL; | |
1532 return; | |
1533 } | |
1534 | |
1535 void ios_delete_pool(glp_tree *tree, IOSPOOL *pool) | |
1536 { /* delete cut pool */ | |
1537 xassert(pool != NULL); | |
1538 ios_clear_pool(tree, pool); | |
1539 xfree(pool); | |
1540 return; | |
1541 } | |
1542 | |
1543 /**********************************************************************/ | |
1544 | |
1545 #if 0 | |
1546 static int refer_to_node(glp_tree *tree, int j) | |
1547 { /* determine node number corresponding to binary variable x[j] or | |
1548 its complement */ | |
1549 glp_prob *mip = tree->mip; | |
1550 int n = mip->n; | |
1551 int *ref; | |
1552 if (j > 0) | |
1553 ref = tree->n_ref; | |
1554 else | |
1555 ref = tree->c_ref, j = - j; | |
1556 xassert(1 <= j && j <= n); | |
1557 if (ref[j] == 0) | |
1558 { /* new node is needed */ | |
1559 SCG *g = tree->g; | |
1560 int n_max = g->n_max; | |
1561 ref[j] = scg_add_nodes(g, 1); | |
1562 if (g->n_max > n_max) | |
1563 { int *save = tree->j_ref; | |
1564 tree->j_ref = xcalloc(1+g->n_max, sizeof(int)); | |
1565 memcpy(&tree->j_ref[1], &save[1], g->n * sizeof(int)); | |
1566 xfree(save); | |
1567 } | |
1568 xassert(ref[j] == g->n); | |
1569 tree->j_ref[ref[j]] = j; | |
1570 xassert(tree->curr != NULL); | |
1571 if (tree->curr->level > 0) tree->curr->own_nn++; | |
1572 } | |
1573 return ref[j]; | |
1574 } | |
1575 #endif | |
1576 | |
1577 #if 0 | |
1578 void ios_add_edge(glp_tree *tree, int j1, int j2) | |
1579 { /* add new edge to the conflict graph */ | |
1580 glp_prob *mip = tree->mip; | |
1581 int n = mip->n; | |
1582 SCGRIB *e; | |
1583 int first, i1, i2; | |
1584 xassert(-n <= j1 && j1 <= +n && j1 != 0); | |
1585 xassert(-n <= j2 && j2 <= +n && j2 != 0); | |
1586 xassert(j1 != j2); | |
1587 /* determine number of the first node, which was added for the | |
1588 current subproblem */ | |
1589 xassert(tree->curr != NULL); | |
1590 first = tree->g->n - tree->curr->own_nn + 1; | |
1591 /* determine node numbers for both endpoints */ | |
1592 i1 = refer_to_node(tree, j1); | |
1593 i2 = refer_to_node(tree, j2); | |
1594 /* add edge (i1,i2) to the conflict graph */ | |
1595 e = scg_add_edge(tree->g, i1, i2); | |
1596 /* if the current subproblem is not the root and both endpoints | |
1597 were created on some previous levels, save the edge */ | |
1598 if (tree->curr->level > 0 && i1 < first && i2 < first) | |
1599 { IOSRIB *rib; | |
1600 rib = dmp_get_atom(tree->pool, sizeof(IOSRIB)); | |
1601 rib->j1 = j1; | |
1602 rib->j2 = j2; | |
1603 rib->e = e; | |
1604 rib->next = tree->curr->e_ptr; | |
1605 tree->curr->e_ptr = rib; | |
1606 } | |
1607 return; | |
1608 } | |
1609 #endif | |
1610 | |
1611 /* eof */ |