opal/util/bipartite

/* [<][>][^][v][top][bottom][index][help] */
This source file includes following definitions.
check_add64_overflow
edge_constructor
edge_destructor
dump_vec
dump_vec64
dump_flow
get_capacity
set_capacity
free_vertex
opal_bp_graph_create
opal_bp_graph_free
opal_bp_graph_clone
opal_bp_graph_indegree
opal_bp_graph_outdegree
opal_bp_graph_add_edge
opal_bp_graph_add_vertex
opal_bp_graph_order
shrink_flow_matrix
bottleneck_path
opal_bp_graph_bellman_ford
opal_bp_graph_bipartite_to_flow
min_cost_flow_ssp
opal_bp_graph_solve_bipartite_assignment
   1 /*
   2  * Copyright (c) 2014      Cisco Systems, Inc.  All rights reserved.
   3  * Copyright (c) 2017      Amazon.com, Inc. or its affiliates.  All Rights
   4  *                         reserved.
   5  * $COPYRIGHT$
   6  *
   7  * Additional copyrights may follow
   8  *
   9  * $HEADER$
  10  */
  11 
  12 #include "opal_config.h"
  13 
  14 #include <stdlib.h>
  15 #include <stddef.h>
  16 
  17 #include "opal_stdint.h"
  18 #include "opal/constants.h"
  19 #include "opal/class/opal_list.h"
  20 #include "opal/class/opal_pointer_array.h"
  21 #include "opal/util/output.h"
  22 #include "opal/util/error.h"
  23 
  24 #include "opal/util/bipartite_graph.h"
  25 #include "opal/util/bipartite_graph_internal.h"
  26 
  27 #ifndef container_of
  28 #define container_of(ptr, type, member) ( \
  29         (type *)( ((char *)(ptr)) - offsetof(type,member) ))
  30 #endif
  31 
  32 #define GRAPH_DEBUG 0
  33 #if GRAPH_DEBUG
  34 #  define GRAPH_DEBUG_OUT(args) printf(args)
  35 #else
  36 #  define GRAPH_DEBUG_OUT(args) do {} while(0)
  37 #endif
  38 
  39 #define MAX_COST INT64_MAX
  40 
  41 #ifndef MAX
  42 #  define MAX(a,b) ((a) > (b) ? (a) : (b))
  43 #endif
  44 
  45 #ifndef MIN
  46 #  define MIN(a,b) ((a) < (b) ? (a) : (b))
  47 #endif
  48 
  49 #define f(i,j) flow[n*i + j]
  50 
  51 /* ensure that (a+b<=max) */
  52 static inline void check_add64_overflow(int64_t a, int64_t b)
  53 {
  54     assert(!((b > 0) && (a > (INT64_MAX - b))) &&
  55            !((b < 0) && (a < (INT64_MIN - b))));
  56 }
  57 
  58 static void edge_constructor(opal_bp_graph_edge_t *e)
  59 {
  60     OBJ_CONSTRUCT(&e->outbound_li, opal_list_item_t);
  61     OBJ_CONSTRUCT(&e->inbound_li, opal_list_item_t);
  62 }
  63 
  64 static void edge_destructor(opal_bp_graph_edge_t *e)
  65 {
  66     OBJ_DESTRUCT(&e->outbound_li);
  67     OBJ_DESTRUCT(&e->inbound_li);
  68 }
  69 
  70 OBJ_CLASS_DECLARATION(opal_bp_graph_edge_t);
  71 OBJ_CLASS_INSTANCE(opal_bp_graph_edge_t, opal_object_t,
  72                    edge_constructor, edge_destructor);
  73 
  74 static void dump_vec(const char *name, int *vec, int n)
  75     __opal_attribute_unused__;
  76 
  77 static void dump_vec(const char *name, int *vec, int n)
  78 {
  79     int i;
  80     fprintf(stderr, "%s={", name);
  81     for (i = 0; i < n; ++i) {
  82         fprintf(stderr, "[%d]=%2d, ", i, vec[i]);
  83     }
  84     fprintf(stderr, "}\n");
  85 }
  86 
  87 static void dump_vec64(const char *name, int64_t *vec, int n)
  88     __opal_attribute_unused__;
  89 
  90 static void dump_vec64(const char *name, int64_t *vec, int n)
  91 {
  92     int i;
  93     fprintf(stderr, "%s={", name);
  94     for (i = 0; i < n; ++i) {
  95         fprintf(stderr, "[%d]=%2" PRIi64 ", ", i, vec[i]);
  96     }
  97     fprintf(stderr, "}\n");
  98 }
  99 
 100 
 101 static void dump_flow(int *flow, int n)
 102     __opal_attribute_unused__;
 103 
 104 static void dump_flow(int *flow, int n)
 105 {
 106     int u, v;
 107 
 108     fprintf(stderr, "flow={\n");
 109     for (u = 0; u < n; ++u) {
 110         fprintf(stderr, "u=%d| ", u);
 111         for (v = 0; v < n; ++v) {
 112             fprintf(stderr, "%2d,", f(u,v));
 113         }
 114         fprintf(stderr, "\n");
 115     }
 116     fprintf(stderr, "}\n");
 117 }
 118 
 119 
 120 static int get_capacity(opal_bp_graph_t *g, int source, int target)
 121 {
 122     opal_bp_graph_edge_t *e;
 123 
 124     CHECK_VERTEX_RANGE(g, source);
 125     CHECK_VERTEX_RANGE(g, target);
 126 
 127     FOREACH_OUT_EDGE(g, source, e) {
 128         assert(e->source == source);
 129         if (e->target == target) {
 130             return e->capacity;
 131         }
 132     }
 133 
 134     return 0;
 135 }
 136 
 137 static int
 138 set_capacity(opal_bp_graph_t *g, int source, int target, int cap)
 139 {
 140     opal_bp_graph_edge_t *e;
 141 
 142     CHECK_VERTEX_RANGE(g, source);
 143     CHECK_VERTEX_RANGE(g, target);
 144 
 145     FOREACH_OUT_EDGE(g, source, e) {
 146         assert(e->source == source);
 147         if (e->target == target) {
 148             e->capacity = cap;
 149             return OPAL_SUCCESS;
 150         }
 151     }
 152 
 153     return OPAL_ERR_NOT_FOUND;
 154 }
 155 
 156 static void free_vertex(opal_bp_graph_t *g,
 157                         opal_bp_graph_vertex_t *v)
 158 {
 159     if (NULL != v) {
 160         if (NULL != g->v_data_cleanup_fn && NULL != v->v_data) {
 161             g->v_data_cleanup_fn(v->v_data);
 162         }
 163         free(v);
 164     }
 165 }
 166 
 167 int opal_bp_graph_create(opal_bp_graph_cleanup_fn_t v_data_cleanup_fn,
 168                          opal_bp_graph_cleanup_fn_t e_data_cleanup_fn,
 169                          opal_bp_graph_t **g_out)
 170 {
 171     int err;
 172     opal_bp_graph_t *g = NULL;
 173 
 174     if (NULL == g_out) {
 175         return OPAL_ERR_BAD_PARAM;
 176     }
 177     *g_out = NULL;
 178 
 179     g = calloc(1, sizeof(*g));
 180     if (NULL == g) {
 181         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 182         err = OPAL_ERR_OUT_OF_RESOURCE;
 183         goto out_free_g;
 184     }
 185 
 186     g->source_idx = -1;
 187     g->sink_idx = -1;
 188 
 189     g->v_data_cleanup_fn = v_data_cleanup_fn;
 190     g->e_data_cleanup_fn = e_data_cleanup_fn;
 191 
 192     /* now that we essentially have an empty graph, add vertices to it */
 193     OBJ_CONSTRUCT(&g->vertices, opal_pointer_array_t);
 194     err = opal_pointer_array_init(&g->vertices, 0, INT_MAX, 32);
 195     if (OPAL_SUCCESS != err) {
 196         goto out_free_g;
 197     }
 198 
 199     *g_out = g;
 200     return OPAL_SUCCESS;
 201 
 202  out_free_g:
 203     free(g);
 204     return err;
 205 }
 206 
 207 int opal_bp_graph_free(opal_bp_graph_t *g)
 208 {
 209     int i;
 210     opal_bp_graph_edge_t *e, *next;
 211     opal_bp_graph_vertex_t *v;
 212 
 213     /* remove all edges from all out_edges lists */
 214     for (i = 0; i < NUM_VERTICES(g); ++i) {
 215         v = V_ID_TO_PTR(g, i);
 216         LIST_FOREACH_SAFE_CONTAINED(e, next, &v->out_edges,
 217                                     opal_bp_graph_edge_t, outbound_li) {
 218             opal_list_remove_item(&v->out_edges, &e->outbound_li);
 219             OBJ_RELEASE(e);
 220         }
 221     }
 222     /* now remove from all in_edges lists and free the edge */
 223     for (i = 0; i < NUM_VERTICES(g); ++i) {
 224         v = V_ID_TO_PTR(g, i);
 225         LIST_FOREACH_SAFE_CONTAINED(e, next, &v->in_edges,
 226                                     opal_bp_graph_edge_t, inbound_li) {
 227             opal_list_remove_item(&v->in_edges, &e->inbound_li);
 228 
 229             if (NULL != g->e_data_cleanup_fn && NULL != e->e_data) {
 230                 g->e_data_cleanup_fn(e->e_data);
 231             }
 232             OBJ_RELEASE(e);
 233         }
 234 
 235         free_vertex(g, V_ID_TO_PTR(g, i));
 236         opal_pointer_array_set_item(&g->vertices, i, NULL);
 237     }
 238     g->num_vertices = 0;
 239 
 240     OBJ_DESTRUCT(&g->vertices);
 241     free(g);
 242 
 243     return OPAL_SUCCESS;
 244 }
 245 
 246 int opal_bp_graph_clone(const opal_bp_graph_t *g,
 247                         bool copy_user_data,
 248                         opal_bp_graph_t **g_clone_out)
 249 {
 250     int err;
 251     int i;
 252     int index;
 253     opal_bp_graph_t *gx;
 254     opal_bp_graph_edge_t *e;
 255 
 256     if (NULL == g_clone_out) {
 257         return OPAL_ERR_BAD_PARAM;
 258     }
 259     *g_clone_out = NULL;
 260 
 261     if (copy_user_data) {
 262         opal_output(0, "[%s:%d:%s] user data copy requested but not yet supported", 
 263                     __FILE__, __LINE__, __func__);
 264         abort();
 265         return OPAL_ERR_FATAL;
 266     }
 267 
 268     gx = NULL;
 269     err = opal_bp_graph_create(NULL, NULL, &gx);
 270     if (OPAL_SUCCESS != err) {
 271         return err;
 272     }
 273     assert(NULL != gx);
 274 
 275     /* reconstruct all vertices */
 276     for (i = 0; i < NUM_VERTICES(g); ++i) {
 277         err = opal_bp_graph_add_vertex(gx, NULL, &index);
 278         if (OPAL_SUCCESS != err) {
 279             goto out_free_gx;
 280         }
 281         assert(index == i);
 282     }
 283 
 284     /* now reconstruct all the edges (iterate by source vertex only to avoid
 285      * double-adding) */
 286     for (i = 0; i < NUM_VERTICES(g); ++i) {
 287         FOREACH_OUT_EDGE(g, i, e) {
 288             assert(i == e->source);
 289             err = opal_bp_graph_add_edge(gx, e->source, e->target,
 290                                             e->cost, e->capacity, NULL);
 291             if (OPAL_SUCCESS != err) {
 292                 goto out_free_gx;
 293             }
 294         }
 295     }
 296 
 297     *g_clone_out = gx;
 298     return OPAL_SUCCESS;
 299 
 300  out_free_gx:
 301     /* we don't reach in and manipulate gx's state directly, so it should be
 302      * safe to use the standard free function */
 303     opal_bp_graph_free(gx);
 304     return err;
 305 }
 306 
 307 int opal_bp_graph_indegree(const opal_bp_graph_t *g,
 308                            int vertex)
 309 {
 310     opal_bp_graph_vertex_t *v;
 311 
 312     v = V_ID_TO_PTR(g, vertex);
 313     return opal_list_get_size(&v->in_edges);
 314 }
 315 
 316 int opal_bp_graph_outdegree(const opal_bp_graph_t *g,
 317                                int vertex)
 318 {
 319     opal_bp_graph_vertex_t *v;
 320 
 321     v = V_ID_TO_PTR(g, vertex);
 322     return opal_list_get_size(&v->out_edges);
 323 }
 324 
 325 int opal_bp_graph_add_edge(opal_bp_graph_t *g,
 326                            int from,
 327                            int to,
 328                            int64_t cost,
 329                            int capacity,
 330                            void *e_data)
 331 {
 332     opal_bp_graph_edge_t *e;
 333     opal_bp_graph_vertex_t *v_from, *v_to;
 334 
 335     if (from < 0 || from >= NUM_VERTICES(g)) {
 336         return OPAL_ERR_BAD_PARAM;
 337     }
 338     if (to < 0 || to >= NUM_VERTICES(g)) {
 339         return OPAL_ERR_BAD_PARAM;
 340     }
 341     if (cost == MAX_COST) {
 342         return OPAL_ERR_BAD_PARAM;
 343     }
 344     if (capacity < 0) {
 345         /* negative cost is fine, but negative capacity is not currently
 346          * handled appropriately */
 347         return OPAL_ERR_BAD_PARAM;
 348     }
 349     FOREACH_OUT_EDGE(g, from, e) {
 350         assert(e->source == from);
 351         if (e->target == to) {
 352             return OPAL_EXISTS;
 353         }
 354     }
 355 
 356     /* this reference is owned by the out_edges list */
 357     e = OBJ_NEW(opal_bp_graph_edge_t);
 358     if (NULL == e) {
 359         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 360         return OPAL_ERR_OUT_OF_RESOURCE;
 361     }
 362 
 363     e->source   = from;
 364     e->target   = to;
 365     e->cost     = cost;
 366     e->capacity = capacity;
 367     e->e_data   = e_data;
 368 
 369     v_from = V_ID_TO_PTR(g, from);
 370     opal_list_append(&v_from->out_edges, &e->outbound_li);
 371 
 372     OBJ_RETAIN(e); /* ref owned by in_edges list */
 373     v_to = V_ID_TO_PTR(g, to);
 374     opal_list_append(&v_to->in_edges, &e->inbound_li);
 375 
 376     return OPAL_SUCCESS;
 377 }
 378 
 379 int opal_bp_graph_add_vertex(opal_bp_graph_t *g,
 380                              void *v_data,
 381                              int *index_out)
 382 {
 383     opal_bp_graph_vertex_t *v;
 384 
 385     v = calloc(1, sizeof(*v));
 386     if (NULL == v) {
 387         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 388         return OPAL_ERR_OUT_OF_RESOURCE;
 389     }
 390 
 391     /* add to the ptr array early to simplify cleanup in the incredibly rare
 392      * chance that adding fails */
 393     v->v_index = opal_pointer_array_add(&g->vertices, v);
 394     if (-1 == v->v_index) {
 395         free(v);
 396         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 397         return OPAL_ERR_OUT_OF_RESOURCE;
 398     }
 399     assert(v->v_index == g->num_vertices);
 400 
 401     ++g->num_vertices;
 402 
 403     v->v_data = v_data;
 404     OBJ_CONSTRUCT(&v->out_edges, opal_list_t);
 405     OBJ_CONSTRUCT(&v->in_edges, opal_list_t);
 406 
 407     if (NULL != index_out) {
 408         *index_out = v->v_index;
 409     }
 410 
 411     return OPAL_SUCCESS;
 412 }
 413 
 414 int opal_bp_graph_order(const opal_bp_graph_t *g)
 415 {
 416     return NUM_VERTICES(g);
 417 }
 418 
 419 /**
 420  * shrink a flow matrix for old_n vertices to one works for new_n
 421  *
 422  * Takes a matrix stored in a one-dimensional array of size (old_n*old_n) and
 423  * "truncates" it into a dense array of size (new_n*new_n) that only contain
 424  * the flow values for the first new_n vertices.  E.g., it turns this array
 425  * (old_n=5, new_n=3):
 426  *
 427  *    1  2  3  4  5
 428  *    6  7  8  9 10
 429  *   11 12 13 14 15
 430  *   16 17 18 19 20
 431  *   21 22 23 24 25
 432  *
 433  * into this array;
 434  *
 435  *    1  2  3
 436  *    6  7  8
 437  *   11 12 13
 438  */
 439 static void shrink_flow_matrix(int *flow, int old_n, int new_n)
 440 {
 441     int u, v;
 442 
 443     assert(old_n > new_n);
 444 
 445     for (u = 0; u < new_n; ++u) {
 446         for (v = 0; v < new_n; ++v) {
 447             flow[new_n*u + v] = flow[old_n*u + v];
 448         }
 449     }
 450 }
 451 
 452 /**
 453  * Compute the so-called "bottleneck" capacity value for a path "pred" through
 454  * graph "gx".
 455  */
 456 static int
 457 bottleneck_path(
 458                 opal_bp_graph_t *gx,
 459                 int n,
 460                 int *pred)
 461 {
 462     int u, v;
 463     int min;
 464 
 465     min = INT_MAX;
 466     FOREACH_UV_ON_PATH(pred, gx->source_idx, gx->sink_idx, u, v) {
 467         int cap_f_uv = get_capacity(gx, u, v);
 468         min = MIN(min, cap_f_uv);
 469     }
 470 
 471     return min;
 472 }
 473 
 474 
 475 /**
 476  * This routine implements the Bellman-Ford shortest paths algorithm, slightly
 477  * specialized for our forumlation of flow networks:
 478  * http://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
 479  *
 480  * Specifically, it attempts to find the shortest path from "source" to
 481  * "target".  It returns true if such a path was found, false otherwise.  Any
 482  * found path is returned in "pred" as a predecessor chain (i.e., pred[sink]
 483  * is the start of the path and pred[pred[sink]] is its predecessor, etc.).
 484  *
 485  * The contents of "pred" are only valid if this routine returns true.
 486  */
 487 bool opal_bp_graph_bellman_ford(opal_bp_graph_t *gx,
 488                                 int source,
 489                                 int target,
 490                                 int *pred)
 491 {
 492     int64_t *dist;
 493     int i;
 494     int n;
 495     int u, v;
 496     bool found_target = false;
 497 
 498     if (NULL == gx) {
 499         OPAL_ERROR_LOG(OPAL_ERR_BAD_PARAM);
 500         return false;
 501     }
 502     if (NULL == pred) {
 503         OPAL_ERROR_LOG(OPAL_ERR_BAD_PARAM);
 504         return false;
 505     }
 506     if (source < 0 || source >= NUM_VERTICES(gx)) {
 507         return OPAL_ERR_BAD_PARAM;
 508     }
 509     if (target < 0 || target >= NUM_VERTICES(gx)) {
 510         return OPAL_ERR_BAD_PARAM;
 511     }
 512 
 513     /* initialize */
 514     n = opal_bp_graph_order(gx);
 515     dist = malloc(n * sizeof(*dist));
 516     if (NULL == dist) {
 517         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 518         goto out;
 519     }
 520     for (i = 0; i < n; ++i) {
 521         dist[i] = MAX_COST;
 522         pred[i] = -1;
 523     }
 524     dist[source] = 0;
 525 
 526     /* relax repeatedly */
 527     for (i = 1; i < NUM_VERTICES(gx); ++i) {
 528         bool relaxed = false;
 529 #if GRAPH_DEBUG
 530         dump_vec("pred", pred, NUM_VERTICES(gx));
 531         dump_vec64("dist", dist, NUM_VERTICES(gx));
 532 #endif
 533 
 534         for (u = 0; u < NUM_VERTICES(gx); ++u) {
 535             opal_bp_graph_edge_t *e_ptr;
 536 
 537             FOREACH_OUT_EDGE(gx, u, e_ptr) {
 538                 v = e_ptr->target;
 539 
 540                 /* make sure to only construct paths from edges that actually have
 541                  * non-zero capacity */
 542                 if (e_ptr->capacity > 0 &&
 543                     dist[u] != MAX_COST) { /* avoid signed overflow for "infinity" */
 544                     check_add64_overflow(dist[u], e_ptr->cost);
 545                     if ((dist[u] + e_ptr->cost) < dist[v]) {
 546                         dist[v] = dist[u] + e_ptr->cost;
 547                         pred[v] = u;
 548                         relaxed = true;
 549                     }
 550                 }
 551             }
 552         }
 553         /* optimization: stop if an outer iteration did not succeed in
 554          * changing any dist/pred values (already at optimum) */
 555         if (!relaxed) {
 556             GRAPH_DEBUG_OUT(("relaxed==false, breaking out"));
 557             break;
 558         }
 559     }
 560 
 561     /* check for negative-cost cycles */
 562     for (u = 0; u < NUM_VERTICES(gx); ++u) {
 563         opal_bp_graph_edge_t * e_ptr;
 564 
 565         FOREACH_OUT_EDGE(gx, u, e_ptr) {
 566             v = e_ptr->target;
 567             if (e_ptr->capacity > 0 &&
 568                 dist[u] != MAX_COST && /* avoid signed overflow */
 569                 (dist[u] + e_ptr->cost) < dist[v]) {
 570                 opal_output(0, "[%s:%d:%s] negative-weight cycle detected",
 571                             __FILE__, __LINE__, __func__);
 572                 abort();
 573                 goto out;
 574             }
 575         }
 576     }
 577 
 578     if (dist[target] != MAX_COST) {
 579         found_target = true;
 580     }
 581 
 582  out:
 583 #if GRAPH_DEBUG
 584     dump_vec("pred", pred, NUM_VERTICES(gx));
 585 #endif
 586     assert(pred[source] == -1);
 587     free(dist);
 588     GRAPH_DEBUG_OUT(("bellman_ford: found_target=%s", found_target ? "true" : "false"));
 589     return found_target;
 590 }
 591 
 592 /**
 593  * Transform the given connected, bipartite, acyclic digraph into a flow
 594  * network (i.e., add a source and a sink, with the source connected to vertex
 595  * set V1 and the sink connected to vertex set V2).  This also creates
 596  * residual edges suitable for augmenting-path algorithms.  All "source" nodes
 597  * in the original graph are considered to have an output of 1 and "sink"
 598  * nodes can take an input of 1.  The result is that "forward" edges are all
 599  * created with capacity=1, "backward" (residual) edges are created with
 600  * capacity=0.
 601  *
 602  * After this routine, all capacities are "residual capacities" ($c_f$ in the
 603  * literature).
 604  *
 605  * Initial flow throughout the network is assumed to be 0 at all edges.
 606  *
 607  * The graph will be left in an undefined state if an error occurs (though
 608  * freeing it should still be safe).
 609  */
 610 int opal_bp_graph_bipartite_to_flow(opal_bp_graph_t *g)
 611 {
 612     int err;
 613     int order;
 614     int u, v;
 615     int num_left, num_right;
 616 
 617     /* grab size before adding extra vertices */
 618     order = opal_bp_graph_order(g);
 619 
 620     err = opal_bp_graph_add_vertex(g, NULL, &g->source_idx);
 621     if (OPAL_SUCCESS != err) {
 622         return err;
 623     }
 624     err = opal_bp_graph_add_vertex(g, NULL, &g->sink_idx);
 625     if (OPAL_SUCCESS != err) {
 626         return err;
 627     }
 628 
 629     /* The networks we are interested in are bipartite and have edges only
 630      * from one partition to the other partition (none vice versa).  We
 631      * visualize this conventionally with all of the source vertices on the
 632      * left-hand side of an imaginary rendering of the graph and the target
 633      * vertices on the right-hand side of the rendering.  The direction
 634      * "forward" is considered to be moving from left to right.
 635      */
 636     num_left = 0;
 637     num_right = 0;
 638     for (u = 0; u < order; ++u) {
 639         int inbound = opal_bp_graph_indegree(g, u);
 640         int outbound = opal_bp_graph_outdegree(g, u);
 641 
 642         if (inbound > 0 && outbound > 0) {
 643             opal_output(0, "[%s:%d:%s] graph is not (unidirectionally) bipartite",
 644                         __FILE__, __LINE__, __func__);
 645             abort();
 646         }
 647         else if (inbound > 0) {
 648             /* "right" side of the graph, create edges to the sink */
 649             ++num_right;
 650             err = opal_bp_graph_add_edge(g, u, g->sink_idx,
 651                                             0, /* no cost */
 652                                             /*capacity=*/1,
 653                                             /*e_data=*/NULL);
 654             if (OPAL_SUCCESS != err) {
 655                 GRAPH_DEBUG_OUT(("add_edge failed"));
 656                 return err;
 657             }
 658         }
 659         else if (outbound > 0) {
 660             /* "left" side of the graph, create edges to the source */
 661             ++num_left;
 662             err = opal_bp_graph_add_edge(g, g->source_idx, u,
 663                                             0, /* no cost */
 664                                             /*capacity=*/1,
 665                                             /*e_data=*/NULL);
 666             if (OPAL_SUCCESS != err) {
 667                 GRAPH_DEBUG_OUT(("add_edge failed"));
 668                 return err;
 669             }
 670         }
 671     }
 672 
 673     /* it doesn't make sense to extend this graph with a source and sink
 674      * unless */
 675     if (num_right == 0 || num_left == 0) {
 676         return OPAL_ERR_BAD_PARAM;
 677     }
 678 
 679     /* now run through and create "residual" edges as well (i.e., create edges
 680      * in the reverse direction with 0 initial flow and a residual capacity of
 681      * $c_f(u,v)=c(u,v)-f(u,v)$).  Residual edges can exist where no edges
 682      * exist in the original graph.
 683      */
 684     order = opal_bp_graph_order(g); /* need residuals for newly created
 685                                           source/sink edges too */
 686     for (u = 0; u < order; ++u) {
 687         opal_bp_graph_edge_t * e_ptr;
 688         FOREACH_OUT_EDGE(g, u, e_ptr) {
 689             v = e_ptr->target;
 690 
 691             /* (u,v) exists, add (v,u) if not already present.  Cost is
 692              * negative for these edges because "giving back" flow pays us
 693              * back any cost already incurred. */
 694             err = opal_bp_graph_add_edge(g, v, u,
 695                                             -e_ptr->cost,
 696                                             /*capacity=*/0,
 697                                             /*e_data=*/NULL);
 698             if (OPAL_SUCCESS != err && OPAL_EXISTS != err) {
 699                 return err;
 700             }
 701         }
 702     }
 703 
 704     return OPAL_SUCCESS;
 705 }
 706 
 707 /**
 708  * Implements the "Successive Shortest Path" algorithm for computing the
 709  * minimum cost flow problem.  This is a generalized version of the
 710  * Ford-Fulkerson algorithm.  There are two major changes from F-F:
 711  *  1. In addition to capacities and flows, this algorithm pays attention to
 712  *     costs for traversing an edge.  This particular function leaves the
 713  *     caller's costs alone but sets its own capacities.
 714  *  2. Shortest paths are computed using the cost metric.
 715  *
 716  * The algorithm's sketch looks like:
 717  *  1    Transform network G by adding source and sink, create residual edges
 718  *  2    Initial flow x is zero
 719  *  3    while ( Gx contains a path from s to t ) do
 720  *  4        Find any shortest path P from s to t
 721  *  5        Augment current flow x along P
 722  *  6        update Gx
 723  *
 724  * This function mutates the given graph (adding vertices and edges, changing
 725  * capacties, etc.), so callers may wish to clone the graph before calling
 726  * this routine.
 727  *
 728  * The result is an array of (u,v) vertex pairs, where (u,v) is an edge in the
 729  * original graph which has non-zero flow.
 730  *
 731  * Returns OMPI error codes like OPAL_SUCCESS/OPAL_ERR_OUT_OF_RESOURCE.
 732  *
 733  * This version of the algorithm has a theoretical upper bound on its running
 734  * time of O(|V|^2 * |E| * f), where f is essentially the maximum flow in the
 735  * graph.  In our case, f=min(|V1|,|V2|), where V1 and V2 are the two
 736  * constituent sets of the bipartite graph.
 737  *
 738  * This algorithm's performance could probably be improved by modifying it to
 739  * use vertex potentials and Dijkstra's Algorithm instead of Bellman-Ford.
 740  * Normally vertex potentials are needed in order to use Dijkstra's safely,
 741  * but our graphs are constrained enough that this may not be necessary.
 742  * Switching to Dijkstra's implemented with a heap should yield a reduced
 743  * upper bound of O(|V| * |E| * f * log(|V|)).  Let's consider this a future
 744  * enhancement for the time being, since it's not obvious at this point that
 745  * the faster running time will be worth the additional implementation
 746  * complexity.
 747  */
 748 static int min_cost_flow_ssp(opal_bp_graph_t *gx,
 749                              int **flow_out)
 750 {
 751     int err = OPAL_SUCCESS;
 752     int n;
 753     int *pred = NULL;
 754     int *flow = NULL;
 755     int u, v;
 756     int c;
 757 
 758     GRAPH_DEBUG_OUT(("begin min_cost_flow_ssp()"));
 759 
 760     if (NULL == flow_out) {
 761         return OPAL_ERR_BAD_PARAM;
 762     }
 763     *flow_out = NULL;
 764 
 765     n = opal_bp_graph_order(gx);
 766 
 767     pred = malloc(n*sizeof(*pred));
 768     if (NULL == pred) {
 769         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 770         err = OPAL_ERR_OUT_OF_RESOURCE;
 771         goto out_error;
 772     }
 773 
 774     /* "flow" is a 2d matrix of current flow values, all initialized to zero */
 775     flow = calloc(n*n, sizeof(*flow));
 776     if (NULL == flow) {
 777         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 778         err = OPAL_ERR_OUT_OF_RESOURCE;
 779         goto out_error;
 780     }
 781 
 782     /* loop as long as paths exist from source to sink */
 783     while (opal_bp_graph_bellman_ford(gx, gx->source_idx, gx->sink_idx, pred)) {
 784         int cap_f_path;
 785 
 786         /* find any shortest path P from s to t (already present in pred) */
 787         GRAPH_DEBUG_OUT(("start outer iteration of SSP algorithm"));
 788 #if GRAPH_DEBUG
 789         dump_vec("pred", pred, NUM_VERTICES(gx));
 790         dump_flow(flow, n);
 791 #endif
 792 
 793         cap_f_path = bottleneck_path(gx, n, pred);
 794 
 795         /* augment current flow along P */
 796         FOREACH_UV_ON_PATH(pred, gx->source_idx, gx->sink_idx, u, v) {
 797             assert(u == pred[v]);
 798 
 799             f(u,v) = f(u,v) + cap_f_path; /* "forward" edge */
 800             f(v,u) = f(v,u) - cap_f_path; /* residual network edge */
 801 
 802             assert(f(u,v) == -f(v,u)); /* skew symmetry invariant */
 803 
 804             /* update Gx as we go along: decrease capacity by this new
 805              * augmenting flow */
 806             c = get_capacity(gx, u, v) - cap_f_path;
 807             assert(c >= 0);
 808             err = set_capacity(gx, u, v, c);
 809             if (OPAL_SUCCESS != err) {
 810                 opal_output(0, "[%s:%d:%s] unable to set capacity, missing edge?",
 811                             __FILE__, __LINE__, __func__);
 812                 abort();
 813             }
 814 
 815             c = get_capacity(gx, v, u) + cap_f_path;
 816             assert(c >= 0);
 817             err = set_capacity(gx, v, u, c);
 818             if (OPAL_SUCCESS != err) {
 819                 opal_output(0, "[%s:%d:%s] unable to set capacity, missing edge?",
 820                             __FILE__, __LINE__, __func__);
 821                 abort();
 822             }
 823         }
 824     }
 825 
 826  out:
 827     *flow_out = flow;
 828     free(pred);
 829     return err;
 830 
 831  out_error:
 832     free(*flow_out);
 833     GRAPH_DEBUG_OUT(("returning error %d", err));
 834     goto out;
 835 }
 836 
 837 int opal_bp_graph_solve_bipartite_assignment(const opal_bp_graph_t *g,
 838                                              int *num_match_edges_out,
 839                                              int **match_edges_out)
 840 {
 841     int err;
 842     int i;
 843     int u, v;
 844     int n;
 845     int *flow = NULL;
 846     opal_bp_graph_t *gx = NULL;
 847 
 848     if (NULL == match_edges_out || NULL == num_match_edges_out) {
 849         return OPAL_ERR_BAD_PARAM;
 850     }
 851     *num_match_edges_out = 0;
 852     *match_edges_out = NULL;
 853 
 854     /* don't perturb the caller's data structure */
 855     err = opal_bp_graph_clone(g, false, &gx);
 856     if (OPAL_SUCCESS != err) {
 857         GRAPH_DEBUG_OUT(("opal_bp_graph_clone failed"));
 858         goto out;
 859     }
 860 
 861     /* Transform gx into a residual flow network with capacities, a source, a
 862      * sink, and residual edges.  We track the actual flow separately in the
 863      * "flow" matrix.  Initial capacity for every forward edge is 1.  Initial
 864      * capacity for every backward (residual) edge is 0.
 865      *
 866      * For the remainder of this routine (and the ssp routine) the capacities
 867      * refer to residual capacities ($c_f$) not capacities in the original
 868      * graph.  For convenience we adjust all residual capacities as we go
 869      * along rather than recomputing them from the flow and capacities in the
 870      * original graph.  This allows many other graph operations to have no
 871      * direct knowledge of the flow matrix.
 872      */
 873     err = opal_bp_graph_bipartite_to_flow(gx);
 874     if (OPAL_SUCCESS != err) {
 875         GRAPH_DEBUG_OUT(("bipartite_to_flow failed"));
 876         OPAL_ERROR_LOG(err);
 877         return err;
 878     }
 879 
 880     /* Use the SSP algorithm to compute the min-cost flow over this network.
 881      * Edges with non-zero flow in the result should be part of the matching.
 882      *
 883      * Note that the flow array returned is sized for gx, not for g.  Index
 884      * accordingly later on.
 885      */
 886     err = min_cost_flow_ssp(gx, &flow);
 887     if (OPAL_SUCCESS != err) {
 888         GRAPH_DEBUG_OUT(("min_cost_flow_ssp failed"));
 889         return err;
 890     }
 891     assert(NULL != flow);
 892 
 893     /* don't care about new edges in gx, only old edges in g */
 894     n = opal_bp_graph_order(g);
 895 
 896 #if GRAPH_DEBUG
 897     dump_flow(flow, NUM_VERTICES(gx));
 898 #endif
 899     shrink_flow_matrix(flow, opal_bp_graph_order(gx), n);
 900 #if GRAPH_DEBUG
 901     dump_flow(flow, n);
 902 #endif
 903 
 904     for (u = 0; u < n; ++u) {
 905         for (v = 0; v < n; ++v) {
 906             if (f(u,v) > 0) {
 907                 ++(*num_match_edges_out);
 908             }
 909         }
 910     }
 911 
 912     if (0 == *num_match_edges_out) {
 913         /* avoid attempting to allocate a zero-byte buffer */
 914         goto out;
 915     }
 916 
 917     *match_edges_out = malloc(*num_match_edges_out * 2 * sizeof(int));
 918     if (NULL == *match_edges_out) {
 919         *num_match_edges_out = 0;
 920         OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
 921         err = OPAL_ERR_OUT_OF_RESOURCE;
 922         goto out;
 923     }
 924 
 925     i = 0;
 926     for (u = 0; u < n; ++u) {
 927         for (v = 0; v < n; ++v) {
 928             /* flow exists on this edge so include this edge in the matching */
 929             if (f(u,v) > 0) {
 930                 (*match_edges_out)[i++] = u;
 931                 (*match_edges_out)[i++] = v;
 932             }
 933         }
 934     }
 935 
 936  out:
 937     free(flow);
 938     opal_bp_graph_free(gx);
 939     return err;
 940 }
/* [<][>][^][v][top][bottom][index][help] */
root/opal/util/bipartite_graph.c

DEFINITIONS
