glpios01.c 52 KB

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