THESE/codes/graphe/GCmex1.9/maxflow.cpp

   1 /* maxflow.cpp */\r
   2 /*\r
   3     Copyright 2001 Vladimir Kolmogorov (vnk@cs.cornell.edu), Yuri Boykov (yuri@csd.uwo.ca).\r
   4 \r
   5     This program is free software; you can redistribute it and/or modify\r
   6     it under the terms of the GNU General Public License as published by\r
   7     the Free Software Foundation; either version 2 of the License, or\r
   8     (at your option) any later version.\r
   9 \r
  10     This program is distributed in the hope that it will be useful,\r
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of\r
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
  13     GNU General Public License for more details.\r
  14 \r
  15     You should have received a copy of the GNU General Public License\r
  16     along with this program; if not, write to the Free Software\r
  17     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA\r
  18 */\r
  19 \r
  20 \r
  21 #include <stdio.h>\r
  22 #include "graph.h"\r
  23 \r
  24 /*\r
  25         special constants for node->parent\r
  26 */\r
  27 #define TERMINAL ( (arc_forward *) 1 )          /* to terminal */\r
  28 #define ORPHAN   ( (arc_forward *) 2 )          /* orphan */\r
  29 \r
  30 #define INFINITE_D 1000000000           /* infinite distance to the terminal */\r
  31 \r
  32 /***********************************************************************/\r
  33 \r
  34 /*\r
  35         Functions for processing active list.\r
  36         i->next points to the next node in the list\r
  37         (or to i, if i is the last node in the list).\r
  38         If i->next is NULL iff i is not in the list.\r
  39 \r
  40         There are two queues. Active nodes are added\r
  41         to the end of the second queue and read from\r
  42         the front of the first queue. If the first queue\r
  43         is empty, it is replaced by the second queue\r
  44         (and the second queue becomes empty).\r
  45 */\r
  46 \r
  47 inline void Graph::set_active(node *i)\r
  48 {\r
  49         if (!i->next)\r
  50         {\r
  51                 /* it's not in the list yet */\r
  52                 if (queue_last[1]) queue_last[1] -> next = i;\r
  53                 else               queue_first[1]        = i;\r
  54                 queue_last[1] = i;\r
  55                 i -> next = i;\r
  56         }\r
  57 }\r
  58 \r
  59 /*\r
  60         Returns the next active node.\r
  61         If it is connected to the sink, it stays in the list,\r
  62         otherwise it is removed from the list\r
  63 */\r
  64 inline Graph::node * Graph::next_active()\r
  65 {\r
  66         node *i;\r
  67 \r
  68         while ( 1 )\r
  69         {\r
  70                 if (!(i=queue_first[0]))\r
  71                 {\r
  72                         queue_first[0] = i = queue_first[1];\r
  73                         queue_last[0]  = queue_last[1];\r
  74                         queue_first[1] = NULL;\r
  75                         queue_last[1]  = NULL;\r
  76                         if (!i) return NULL;\r
  77                 }\r
  78 \r
  79                 /* remove it from the active list */\r
  80                 if (i->next == i) queue_first[0] = queue_last[0] = NULL;\r
  81                 else              queue_first[0] = i -> next;\r
  82                 i -> next = NULL;\r
  83 \r
  84                 /* a node in the list is active iff it has a parent */\r
  85                 if (i->parent) return i;\r
  86         }\r
  87 }\r
  88 \r
  89 /***********************************************************************/\r
  90 \r
  91 void Graph::maxflow_init()\r
  92 {\r
  93         node *i;\r
  94         node_block *nb;\r
  95 \r
  96         queue_first[0] = queue_last[0] = NULL;\r
  97         queue_first[1] = queue_last[1] = NULL;\r
  98         orphan_first = NULL;\r
  99 \r
 100         for (nb=node_block_first; nb; nb=nb->next)\r
 101         for (i=&nb->nodes[0]; i<nb->current; i++)\r
 102         {\r
 103                 i -> next = NULL;\r
 104                 i -> TS = 0;\r
 105                 if (i->tr_cap > 0)\r
 106                 {\r
 107                         /* i is connected to the source */\r
 108                         i -> is_sink = 0;\r
 109                         i -> parent = TERMINAL;\r
 110                         set_active(i);\r
 111                         i -> TS = 0;\r
 112                         i -> DIST = 1;\r
 113                 }\r
 114                 else if (i->tr_cap < 0)\r
 115                 {\r
 116                         /* i is connected to the sink */\r
 117                         i -> is_sink = 1;\r
 118                         i -> parent = TERMINAL;\r
 119                         set_active(i);\r
 120                         i -> TS = 0;\r
 121                         i -> DIST = 1;\r
 122                 }\r
 123                 else\r
 124                 {\r
 125                         i -> parent = NULL;\r
 126                 }\r
 127         }\r
 128         TIME = 0;\r
 129 }\r
 130 \r
 131 /***********************************************************************/\r
 132 \r
 133 void Graph::augment(node *s_start, node *t_start, captype *cap_middle, captype *rev_cap_middle)\r
 134 {\r
 135         node *i;\r
 136         arc_forward *a;\r
 137         captype bottleneck;\r
 138         nodeptr *np;\r
 139 \r
 140 \r
 141         /* 1. Finding bottleneck capacity */\r
 142         /* 1a - the source tree */\r
 143         bottleneck = *cap_middle;\r
 144         for (i=s_start; ; )\r
 145         {\r
 146                 a = i -> parent;\r
 147                 if (a == TERMINAL) break;\r
 148                 if (IS_ODD(a))\r
 149                 {\r
 150                         a = MAKE_EVEN(a);\r
 151                         if (bottleneck > a->r_cap) bottleneck = a -> r_cap;\r
 152                         i = NEIGHBOR_NODE_REV(i, a -> shift);\r
 153                 }\r
 154                 else\r
 155                 {\r
 156                         if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;\r
 157                         i = NEIGHBOR_NODE(i, a -> shift);\r
 158                 }\r
 159         }\r
 160         if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;\r
 161         /* 1b - the sink tree */\r
 162         for (i=t_start; ; )\r
 163         {\r
 164                 a = i -> parent;\r
 165                 if (a == TERMINAL) break;\r
 166                 if (IS_ODD(a))\r
 167                 {\r
 168                         a = MAKE_EVEN(a);\r
 169                         if (bottleneck > a->r_rev_cap) bottleneck = a -> r_rev_cap;\r
 170                         i = NEIGHBOR_NODE_REV(i, a -> shift);\r
 171                 }\r
 172                 else\r
 173                 {\r
 174                         if (bottleneck > a->r_cap) bottleneck = a -> r_cap;\r
 175                         i = NEIGHBOR_NODE(i, a -> shift);\r
 176                 }\r
 177         }\r
 178         if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;\r
 179 \r
 180 \r
 181         /* 2. Augmenting */\r
 182         /* 2a - the source tree */\r
 183         *rev_cap_middle += bottleneck;\r
 184         *cap_middle -= bottleneck;\r
 185         for (i=s_start; ; )\r
 186         {\r
 187                 a = i -> parent;\r
 188                 if (a == TERMINAL) break;\r
 189                 if (IS_ODD(a))\r
 190                 {\r
 191                         a = MAKE_EVEN(a);\r
 192                         a -> r_rev_cap += bottleneck;\r
 193                         a -> r_cap -= bottleneck;\r
 194                         if (!a->r_cap)\r
 195                         {\r
 196                                 /* add i to the adoption list */\r
 197                                 i -> parent = ORPHAN;\r
 198                                 np = nodeptr_block -> New();\r
 199                                 np -> ptr = i;\r
 200                                 np -> next = orphan_first;\r
 201                                 orphan_first = np;\r
 202                         }\r
 203                         i = NEIGHBOR_NODE_REV(i, a -> shift);\r
 204                 }\r
 205                 else\r
 206                 {\r
 207                         a -> r_cap += bottleneck;\r
 208                         a -> r_rev_cap -= bottleneck;\r
 209                         if (!a->r_rev_cap)\r
 210                         {\r
 211                                 /* add i to the adoption list */\r
 212                                 i -> parent = ORPHAN;\r
 213                                 np = nodeptr_block -> New();\r
 214                                 np -> ptr = i;\r
 215                                 np -> next = orphan_first;\r
 216                                 orphan_first = np;\r
 217                         }\r
 218                         i = NEIGHBOR_NODE(i, a -> shift);\r
 219                 }\r
 220         }\r
 221         i -> tr_cap -= bottleneck;\r
 222         if (!i->tr_cap)\r
 223         {\r
 224                 /* add i to the adoption list */\r
 225                 i -> parent = ORPHAN;\r
 226                 np = nodeptr_block -> New();\r
 227                 np -> ptr = i;\r
 228                 np -> next = orphan_first;\r
 229                 orphan_first = np;\r
 230         }\r
 231         /* 2b - the sink tree */\r
 232         for (i=t_start; ; )\r
 233         {\r
 234                 a = i -> parent;\r
 235                 if (a == TERMINAL) break;\r
 236                 if (IS_ODD(a))\r
 237                 {\r
 238                         a = MAKE_EVEN(a);\r
 239                         a -> r_cap += bottleneck;\r
 240                         a -> r_rev_cap -= bottleneck;\r
 241                         if (!a->r_rev_cap)\r
 242                         {\r
 243                                 /* add i to the adoption list */\r
 244                                 i -> parent = ORPHAN;\r
 245                                 np = nodeptr_block -> New();\r
 246                                 np -> ptr = i;\r
 247                                 np -> next = orphan_first;\r
 248                                 orphan_first = np;\r
 249                         }\r
 250                         i = NEIGHBOR_NODE_REV(i, a -> shift);\r
 251                 }\r
 252                 else\r
 253                 {\r
 254                         a -> r_rev_cap += bottleneck;\r
 255                         a -> r_cap -= bottleneck;\r
 256                         if (!a->r_cap)\r
 257                         {\r
 258                                 /* add i to the adoption list */\r
 259                                 i -> parent = ORPHAN;\r
 260                                 np = nodeptr_block -> New();\r
 261                                 np -> ptr = i;\r
 262                                 np -> next = orphan_first;\r
 263                                 orphan_first = np;\r
 264                         }\r
 265                         i = NEIGHBOR_NODE(i, a -> shift);\r
 266                 }\r
 267         }\r
 268         i -> tr_cap += bottleneck;\r
 269         if (!i->tr_cap)\r
 270         {\r
 271                 /* add i to the adoption list */\r
 272                 i -> parent = ORPHAN;\r
 273                 np = nodeptr_block -> New();\r
 274                 np -> ptr = i;\r
 275                 np -> next = orphan_first;\r
 276                 orphan_first = np;\r
 277         }\r
 278 \r
 279 \r
 280         flow += bottleneck;\r
 281 }\r
 282 \r
 283 /***********************************************************************/\r
 284 \r
 285 void Graph::process_source_orphan(node *i)\r
 286 {\r
 287         node *j;\r
 288         arc_forward *a0_for, *a0_for_first, *a0_for_last;\r
 289         arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;\r
 290         arc_forward *a0_min = NULL, *a;\r
 291         nodeptr *np;\r
 292         int d, d_min = INFINITE_D;\r
 293 \r
 294         /* trying to find a new parent */\r
 295         a0_for_first = i -> first_out;\r
 296         if (IS_ODD(a0_for_first))\r
 297         {\r
 298                 a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);\r
 299                 a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);\r
 300         }\r
 301         else a0_for_last = (i + 1) -> first_out;\r
 302         a0_rev_first = i -> first_in;\r
 303         if (IS_ODD(a0_rev_first))\r
 304         {\r
 305                 a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);\r
 306                 a0_rev_last  = (arc_reverse *) ((a0_rev_first ++) -> sister);\r
 307         }\r
 308         else a0_rev_last = (i + 1) -> first_in;\r
 309 \r
 310 \r
 311         for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)\r
 312         if (a0_for->r_rev_cap)\r
 313         {\r
 314                 j = NEIGHBOR_NODE(i, a0_for -> shift);\r
 315                 if (!j->is_sink && (a=j->parent))\r
 316                 {\r
 317                         /* checking the origin of j */\r
 318                         d = 0;\r
 319                         while ( 1 )\r
 320                         {\r
 321                                 if (j->TS == TIME)\r
 322                                 {\r
 323                                         d += j -> DIST;\r
 324                                         break;\r
 325                                 }\r
 326                                 a = j -> parent;\r
 327                                 d ++;\r
 328                                 if (a==TERMINAL)\r
 329                                 {\r
 330                                         j -> TS = TIME;\r
 331                                         j -> DIST = 1;\r
 332                                         break;\r
 333                                 }\r
 334                                 if (a==ORPHAN) { d = INFINITE_D; break; }\r
 335                                 if (IS_ODD(a))\r
 336                                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 337                                 else\r
 338                                         j = NEIGHBOR_NODE(j, a -> shift);\r
 339                         }\r
 340                         if (d<INFINITE_D) /* j originates from the source - done */\r
 341                         {\r
 342                                 if (d<d_min)\r
 343                                 {\r
 344                                         a0_min = a0_for;\r
 345                                         d_min = d;\r
 346                                 }\r
 347                                 /* set marks along the path */\r
 348                                 for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; )\r
 349                                 {\r
 350                                         j -> TS = TIME;\r
 351                                         j -> DIST = d --;\r
 352                                         a = j->parent;\r
 353                                         if (IS_ODD(a))\r
 354                                                 j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 355                                         else\r
 356                                                 j = NEIGHBOR_NODE(j, a -> shift);\r
 357                                 }\r
 358                         }\r
 359                 }\r
 360         }\r
 361         for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)\r
 362         {\r
 363                 a0_for = a0_rev -> sister;\r
 364                 if (a0_for->r_cap)\r
 365                 {\r
 366                         j = NEIGHBOR_NODE_REV(i, a0_for -> shift);\r
 367                         if (!j->is_sink && (a=j->parent))\r
 368                         {\r
 369                                 /* checking the origin of j */\r
 370                                 d = 0;\r
 371                                 while ( 1 )\r
 372                                 {\r
 373                                         if (j->TS == TIME)\r
 374                                         {\r
 375                                                 d += j -> DIST;\r
 376                                                 break;\r
 377                                         }\r
 378                                         a = j -> parent;\r
 379                                         d ++;\r
 380                                         if (a==TERMINAL)\r
 381                                         {\r
 382                                                 j -> TS = TIME;\r
 383                                                 j -> DIST = 1;\r
 384                                                 break;\r
 385                                         }\r
 386                                         if (a==ORPHAN) { d = INFINITE_D; break; }\r
 387                                         if (IS_ODD(a))\r
 388                                                 j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 389                                         else\r
 390                                                 j = NEIGHBOR_NODE(j, a -> shift);\r
 391                                 }\r
 392                                 if (d<INFINITE_D) /* j originates from the source - done */\r
 393                                 {\r
 394                                         if (d<d_min)\r
 395                                         {\r
 396                                                 a0_min = MAKE_ODD(a0_for);\r
 397                                                 d_min = d;\r
 398                                         }\r
 399                                         /* set marks along the path */\r
 400                                         for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; )\r
 401                                         {\r
 402                                                 j -> TS = TIME;\r
 403                                                 j -> DIST = d --;\r
 404                                                 a = j->parent;\r
 405                                                 if (IS_ODD(a))\r
 406                                                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 407                                                 else\r
 408                                                         j = NEIGHBOR_NODE(j, a -> shift);\r
 409                                         }\r
 410                                 }\r
 411                         }\r
 412                 }\r
 413         }\r
 414 \r
 415         if (i->parent = a0_min)\r
 416         {\r
 417                 i -> TS = TIME;\r
 418                 i -> DIST = d_min + 1;\r
 419         }\r
 420         else\r
 421         {\r
 422                 /* no parent is found */\r
 423                 i -> TS = 0;\r
 424 \r
 425                 /* process neighbors */\r
 426                 for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)\r
 427                 {\r
 428                         j = NEIGHBOR_NODE(i, a0_for -> shift);\r
 429                         if (!j->is_sink && (a=j->parent))\r
 430                         {\r
 431                                 if (a0_for->r_rev_cap) set_active(j);\r
 432                                 if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)\r
 433                                 {\r
 434                                         /* add j to the adoption list */\r
 435                                         j -> parent = ORPHAN;\r
 436                                         np = nodeptr_block -> New();\r
 437                                         np -> ptr = j;\r
 438                                         if (orphan_last) orphan_last -> next = np;\r
 439                                         else             orphan_first        = np;\r
 440                                         orphan_last = np;\r
 441                                         np -> next = NULL;\r
 442                                 }\r
 443                         }\r
 444                 }\r
 445                 for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)\r
 446                 {\r
 447                         a0_for = a0_rev -> sister;\r
 448                         j = NEIGHBOR_NODE_REV(i, a0_for -> shift);\r
 449                         if (!j->is_sink && (a=j->parent))\r
 450                         {\r
 451                                 if (a0_for->r_cap) set_active(j);\r
 452                                 if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)\r
 453                                 {\r
 454                                         /* add j to the adoption list */\r
 455                                         j -> parent = ORPHAN;\r
 456                                         np = nodeptr_block -> New();\r
 457                                         np -> ptr = j;\r
 458                                         if (orphan_last) orphan_last -> next = np;\r
 459                                         else             orphan_first        = np;\r
 460                                         orphan_last = np;\r
 461                                         np -> next = NULL;\r
 462                                 }\r
 463                         }\r
 464                 }\r
 465         }\r
 466 }\r
 467 \r
 468 void Graph::process_sink_orphan(node *i)\r
 469 {\r
 470         node *j;\r
 471         arc_forward *a0_for, *a0_for_first, *a0_for_last;\r
 472         arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;\r
 473         arc_forward *a0_min = NULL, *a;\r
 474         nodeptr *np;\r
 475         int d, d_min = INFINITE_D;\r
 476 \r
 477         /* trying to find a new parent */\r
 478         a0_for_first = i -> first_out;\r
 479         if (IS_ODD(a0_for_first))\r
 480         {\r
 481                 a0_for_first = (arc_forward *) (((char *)a0_for_first) + 1);\r
 482                 a0_for_last = (arc_forward *) ((a0_for_first ++) -> shift);\r
 483         }\r
 484         else a0_for_last = (i + 1) -> first_out;\r
 485         a0_rev_first = i -> first_in;\r
 486         if (IS_ODD(a0_rev_first))\r
 487         {\r
 488                 a0_rev_first = (arc_reverse *) (((char *)a0_rev_first) + 1);\r
 489                 a0_rev_last  = (arc_reverse *) ((a0_rev_first ++) -> sister);\r
 490         }\r
 491         else a0_rev_last = (i + 1) -> first_in;\r
 492 \r
 493 \r
 494         for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)\r
 495         if (a0_for->r_cap)\r
 496         {\r
 497                 j = NEIGHBOR_NODE(i, a0_for -> shift);\r
 498                 if (j->is_sink && (a=j->parent))\r
 499                 {\r
 500                         /* checking the origin of j */\r
 501                         d = 0;\r
 502                         while ( 1 )\r
 503                         {\r
 504                                 if (j->TS == TIME)\r
 505                                 {\r
 506                                         d += j -> DIST;\r
 507                                         break;\r
 508                                 }\r
 509                                 a = j -> parent;\r
 510                                 d ++;\r
 511                                 if (a==TERMINAL)\r
 512                                 {\r
 513                                         j -> TS = TIME;\r
 514                                         j -> DIST = 1;\r
 515                                         break;\r
 516                                 }\r
 517                                 if (a==ORPHAN) { d = INFINITE_D; break; }\r
 518                                 if (IS_ODD(a))\r
 519                                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 520                                 else\r
 521                                         j = NEIGHBOR_NODE(j, a -> shift);\r
 522                         }\r
 523                         if (d<INFINITE_D) /* j originates from the sink - done */\r
 524                         {\r
 525                                 if (d<d_min)\r
 526                                 {\r
 527                                         a0_min = a0_for;\r
 528                                         d_min = d;\r
 529                                 }\r
 530                                 /* set marks along the path */\r
 531                                 for (j=NEIGHBOR_NODE(i, a0_for->shift); j->TS!=TIME; )\r
 532                                 {\r
 533                                         j -> TS = TIME;\r
 534                                         j -> DIST = d --;\r
 535                                         a = j->parent;\r
 536                                         if (IS_ODD(a))\r
 537                                                 j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 538                                         else\r
 539                                                 j = NEIGHBOR_NODE(j, a -> shift);\r
 540                                 }\r
 541                         }\r
 542                 }\r
 543         }\r
 544         for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)\r
 545         {\r
 546                 a0_for = a0_rev -> sister;\r
 547                 if (a0_for->r_rev_cap)\r
 548                 {\r
 549                         j = NEIGHBOR_NODE_REV(i, a0_for -> shift);\r
 550                         if (j->is_sink && (a=j->parent))\r
 551                         {\r
 552                                 /* checking the origin of j */\r
 553                                 d = 0;\r
 554                                 while ( 1 )\r
 555                                 {\r
 556                                         if (j->TS == TIME)\r
 557                                         {\r
 558                                                 d += j -> DIST;\r
 559                                                 break;\r
 560                                         }\r
 561                                         a = j -> parent;\r
 562                                         d ++;\r
 563                                         if (a==TERMINAL)\r
 564                                         {\r
 565                                                 j -> TS = TIME;\r
 566                                                 j -> DIST = 1;\r
 567                                                 break;\r
 568                                         }\r
 569                                         if (a==ORPHAN) { d = INFINITE_D; break; }\r
 570                                         if (IS_ODD(a))\r
 571                                                 j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 572                                         else\r
 573                                                 j = NEIGHBOR_NODE(j, a -> shift);\r
 574                                 }\r
 575                                 if (d<INFINITE_D) /* j originates from the sink - done */\r
 576                                 {\r
 577                                         if (d<d_min)\r
 578                                         {\r
 579                                                 a0_min = MAKE_ODD(a0_for);\r
 580                                                 d_min = d;\r
 581                                         }\r
 582                                         /* set marks along the path */\r
 583                                         for (j=NEIGHBOR_NODE_REV(i,a0_for->shift); j->TS!=TIME; )\r
 584                                         {\r
 585                                                 j -> TS = TIME;\r
 586                                                 j -> DIST = d --;\r
 587                                                 a = j->parent;\r
 588                                                 if (IS_ODD(a))\r
 589                                                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);\r
 590                                                 else\r
 591                                                         j = NEIGHBOR_NODE(j, a -> shift);\r
 592                                         }\r
 593                                 }\r
 594                         }\r
 595                 }\r
 596         }\r
 597 \r
 598         if (i->parent = a0_min)\r
 599         {\r
 600                 i -> TS = TIME;\r
 601                 i -> DIST = d_min + 1;\r
 602         }\r
 603         else\r
 604         {\r
 605                 /* no parent is found */\r
 606                 i -> TS = 0;\r
 607 \r
 608                 /* process neighbors */\r
 609                 for (a0_for=a0_for_first; a0_for<a0_for_last; a0_for++)\r
 610                 {\r
 611                         j = NEIGHBOR_NODE(i, a0_for -> shift);\r
 612                         if (j->is_sink && (a=j->parent))\r
 613                         {\r
 614                                 if (a0_for->r_cap) set_active(j);\r
 615                                 if (a!=TERMINAL && a!=ORPHAN && IS_ODD(a) && NEIGHBOR_NODE_REV(j, MAKE_EVEN(a)->shift)==i)\r
 616                                 {\r
 617                                         /* add j to the adoption list */\r
 618                                         j -> parent = ORPHAN;\r
 619                                         np = nodeptr_block -> New();\r
 620                                         np -> ptr = j;\r
 621                                         if (orphan_last) orphan_last -> next = np;\r
 622                                         else             orphan_first        = np;\r
 623                                         orphan_last = np;\r
 624                                         np -> next = NULL;\r
 625                                 }\r
 626                         }\r
 627                 }\r
 628                 for (a0_rev=a0_rev_first; a0_rev<a0_rev_last; a0_rev++)\r
 629                 {\r
 630                         a0_for = a0_rev -> sister;\r
 631                         j = NEIGHBOR_NODE_REV(i, a0_for -> shift);\r
 632                         if (j->is_sink && (a=j->parent))\r
 633                         {\r
 634                                 if (a0_for->r_rev_cap) set_active(j);\r
 635                                 if (a!=TERMINAL && a!=ORPHAN && !IS_ODD(a) && NEIGHBOR_NODE(j, a->shift)==i)\r
 636                                 {\r
 637                                         /* add j to the adoption list */\r
 638                                         j -> parent = ORPHAN;\r
 639                                         np = nodeptr_block -> New();\r
 640                                         np -> ptr = j;\r
 641                                         if (orphan_last) orphan_last -> next = np;\r
 642                                         else             orphan_first        = np;\r
 643                                         orphan_last = np;\r
 644                                         np -> next = NULL;\r
 645                                 }\r
 646                         }\r
 647                 }\r
 648         }\r
 649 }\r
 650 \r
 651 /***********************************************************************/\r
 652 \r
 653 Graph::flowtype Graph::maxflow()\r
 654 {\r
 655         node *i, *j, *current_node = NULL, *s_start, *t_start;\r
 656         captype *cap_middle, *rev_cap_middle;\r
 657         arc_forward *a_for, *a_for_first, *a_for_last;\r
 658         arc_reverse *a_rev, *a_rev_first, *a_rev_last;\r
 659         nodeptr *np, *np_next;\r
 660 \r
 661         prepare_graph();\r
 662         maxflow_init();\r
 663         nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);\r
 664 \r
 665         while ( 1 )\r
 666         {\r
 667                 if (i=current_node)\r
 668                 {\r
 669                         i -> next = NULL; /* remove active flag */\r
 670                         if (!i->parent) i = NULL;\r
 671                 }\r
 672                 if (!i)\r
 673                 {\r
 674                         if (!(i = next_active())) break;\r
 675                 }\r
 676 \r
 677                 /* growth */\r
 678                 s_start = NULL;\r
 679 \r
 680                 a_for_first = i -> first_out;\r
 681                 if (IS_ODD(a_for_first))\r
 682                 {\r
 683                         a_for_first = (arc_forward *) (((char *)a_for_first) + 1);\r
 684                         a_for_last = (arc_forward *) ((a_for_first ++) -> shift);\r
 685                 }\r
 686                 else a_for_last = (i + 1) -> first_out;\r
 687                 a_rev_first = i -> first_in;\r
 688                 if (IS_ODD(a_rev_first))\r
 689                 {\r
 690                         a_rev_first = (arc_reverse *) (((char *)a_rev_first) + 1);\r
 691                         a_rev_last = (arc_reverse *) ((a_rev_first ++) -> sister);\r
 692                 }\r
 693                 else a_rev_last = (i + 1) -> first_in;\r
 694 \r
 695                 if (!i->is_sink)\r
 696                 {\r
 697                         /* grow source tree */\r
 698                         for (a_for=a_for_first; a_for<a_for_last; a_for++)\r
 699                         if (a_for->r_cap)\r
 700                         {\r
 701                                 j = NEIGHBOR_NODE(i, a_for -> shift);\r
 702                                 if (!j->parent)\r
 703                                 {\r
 704                                         j -> is_sink = 0;\r
 705                                         j -> parent = MAKE_ODD(a_for);\r
 706                                         j -> TS = i -> TS;\r
 707                                         j -> DIST = i -> DIST + 1;\r
 708                                         set_active(j);\r
 709                                 }\r
 710                                 else if (j->is_sink)\r
 711                                 {\r
 712                                         s_start = i;\r
 713                                         t_start = j;\r
 714                                         cap_middle     = & ( a_for -> r_cap );\r
 715                                         rev_cap_middle = & ( a_for -> r_rev_cap );\r
 716                                         break;\r
 717                                 }\r
 718                                 else if (j->TS <= i->TS &&\r
 719                                          j->DIST > i->DIST)\r
 720                                 {\r
 721                                         /* heuristic - trying to make the distance from j to the source shorter */\r
 722                                         j -> parent = MAKE_ODD(a_for);\r
 723                                         j -> TS = i -> TS;\r
 724                                         j -> DIST = i -> DIST + 1;\r
 725                                 }\r
 726                         }\r
 727                         if (!s_start)\r
 728                         for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)\r
 729                         {\r
 730                                 a_for = a_rev -> sister;\r
 731                                 if (a_for->r_rev_cap)\r
 732                                 {\r
 733                                         j = NEIGHBOR_NODE_REV(i, a_for -> shift);\r
 734                                         if (!j->parent)\r
 735                                         {\r
 736                                                 j -> is_sink = 0;\r
 737                                                 j -> parent = a_for;\r
 738                                                 j -> TS = i -> TS;\r
 739                                                 j -> DIST = i -> DIST + 1;\r
 740                                                 set_active(j);\r
 741                                         }\r
 742                                         else if (j->is_sink)\r
 743                                         {\r
 744                                                 s_start = i;\r
 745                                                 t_start = j;\r
 746                                                 cap_middle     = & ( a_for -> r_rev_cap );\r
 747                                                 rev_cap_middle = & ( a_for -> r_cap );\r
 748                                                 break;\r
 749                                         }\r
 750                                         else if (j->TS <= i->TS &&\r
 751                                                          j->DIST > i->DIST)\r
 752                                         {\r
 753                                                 /* heuristic - trying to make the distance from j to the source shorter */\r
 754                                                 j -> parent = a_for;\r
 755                                                 j -> TS = i -> TS;\r
 756                                                 j -> DIST = i -> DIST + 1;\r
 757                                         }\r
 758                                 }\r
 759                         }\r
 760                 }\r
 761                 else\r
 762                 {\r
 763                         /* grow sink tree */\r
 764                         for (a_for=a_for_first; a_for<a_for_last; a_for++)\r
 765                         if (a_for->r_rev_cap)\r
 766                         {\r
 767                                 j = NEIGHBOR_NODE(i, a_for -> shift);\r
 768                                 if (!j->parent)\r
 769                                 {\r
 770                                         j -> is_sink = 1;\r
 771                                         j -> parent = MAKE_ODD(a_for);\r
 772                                         j -> TS = i -> TS;\r
 773                                         j -> DIST = i -> DIST + 1;\r
 774                                         set_active(j);\r
 775                                 }\r
 776                                 else if (!j->is_sink)\r
 777                                 {\r
 778                                         s_start = j;\r
 779                                         t_start = i;\r
 780                                         cap_middle     = & ( a_for -> r_rev_cap );\r
 781                                         rev_cap_middle = & ( a_for -> r_cap );\r
 782                                         break;\r
 783                                 }\r
 784                                 else if (j->TS <= i->TS &&\r
 785                                          j->DIST > i->DIST)\r
 786                                 {\r
 787                                         /* heuristic - trying to make the distance from j to the sink shorter */\r
 788                                         j -> parent = MAKE_ODD(a_for);\r
 789                                         j -> TS = i -> TS;\r
 790                                         j -> DIST = i -> DIST + 1;\r
 791                                 }\r
 792                         }\r
 793                         for (a_rev=a_rev_first; a_rev<a_rev_last; a_rev++)\r
 794                         {\r
 795                                 a_for = a_rev -> sister;\r
 796                                 if (a_for->r_cap)\r
 797                                 {\r
 798                                         j = NEIGHBOR_NODE_REV(i, a_for -> shift);\r
 799                                         if (!j->parent)\r
 800                                         {\r
 801                                                 j -> is_sink = 1;\r
 802                                                 j -> parent = a_for;\r
 803                                                 j -> TS = i -> TS;\r
 804                                                 j -> DIST = i -> DIST + 1;\r
 805                                                 set_active(j);\r
 806                                         }\r
 807                                         else if (!j->is_sink)\r
 808                                         {\r
 809                                                 s_start = j;\r
 810                                                 t_start = i;\r
 811                                                 cap_middle     = & ( a_for -> r_cap );\r
 812                                                 rev_cap_middle = & ( a_for -> r_rev_cap );\r
 813                                                 break;\r
 814                                         }\r
 815                                         else if (j->TS <= i->TS &&\r
 816                                                          j->DIST > i->DIST)\r
 817                                         {\r
 818                                                 /* heuristic - trying to make the distance from j to the sink shorter */\r
 819                                                 j -> parent = a_for;\r
 820                                                 j -> TS = i -> TS;\r
 821                                                 j -> DIST = i -> DIST + 1;\r
 822                                         }\r
 823                                 }\r
 824                         }\r
 825                 }\r
 826 \r
 827                 TIME ++;\r
 828 \r
 829                 if (s_start)\r
 830                 {\r
 831                         i -> next = i; /* set active flag */\r
 832                         current_node = i;\r
 833 \r
 834                         /* augmentation */\r
 835                         augment(s_start, t_start, cap_middle, rev_cap_middle);\r
 836                         /* augmentation end */\r
 837 \r
 838                         /* adoption */\r
 839                         while (np=orphan_first)\r
 840                         {\r
 841                                 np_next = np -> next;\r
 842                                 np -> next = NULL;\r
 843 \r
 844                                 while (np=orphan_first)\r
 845                                 {\r
 846                                         orphan_first = np -> next;\r
 847                                         i = np -> ptr;\r
 848                                         nodeptr_block -> Delete(np);\r
 849                                         if (!orphan_first) orphan_last = NULL;\r
 850                                         if (i->is_sink) process_sink_orphan(i);\r
 851                                         else            process_source_orphan(i);\r
 852                                 }\r
 853 \r
 854                                 orphan_first = np_next;\r
 855                         }\r
 856                         /* adoption end */\r
 857                 }\r
 858                 else current_node = NULL;\r
 859         }\r
 860 \r
 861         delete nodeptr_block;\r
 862 \r
 863         return flow;\r
 864 }\r
 865 \r
 866 /***********************************************************************/\r
 867 \r
 868 Graph::termtype Graph::what_segment(node_id i)\r
 869 {\r
 870         if (((node*)i)->parent && !((node*)i)->is_sink) return SOURCE;\r
 871         return SINK;\r
 872 }\r
 873 \r