/*
 * Copyright (c) 2014-2022 Cisco Systems, Inc.  All rights reserved
 * Copyright (c) 2017-2019 Amazon.com, Inc. or its affiliates.  All Rights
 *                         reserved.
 * $COPYRIGHT$
 *
 * Additional copyrights may follow
 *
 * $HEADER$
 */

#include "opal_config.h"

#include <stddef.h>
#include <stdlib.h>

#include "opal/class/opal_list.h"
#include "opal/class/opal_pointer_array.h"
#include "opal/constants.h"
#include "opal/util/error.h"
#include "opal/util/misc.h"
#include "opal/util/output.h"
#include "opal_stdint.h"

#include "opal/util/bipartite_graph.h"
#include "opal/util/bipartite_graph_internal.h"

#define GRAPH_DEBUG 0
#if GRAPH_DEBUG
#    define GRAPH_DEBUG_OUT(args) printf(args)
#else
#    define GRAPH_DEBUG_OUT(args) \
        do {                      \
        } while (0)
#endif

#define MAX_COST INT64_MAX

#ifndef MAX
#    define MAX(a, b) ((a) > (b) ? (a) : (b))
#endif

#ifndef MIN
#    define MIN(a, b) ((a) < (b) ? (a) : (b))
#endif

#define f(i, j) flow[n * i + j]

/* ensure that (a+b<=max) */
static inline void check_add64_overflow(int64_t a, int64_t b)
{
    assert(!((b > 0) && (a > (INT64_MAX - b))) && !((b < 0) && (a < (INT64_MIN - b))));
}

static void edge_constructor(opal_bp_graph_edge_t *e)
{
    OBJ_CONSTRUCT(&e->outbound_li, opal_list_item_t);
    OBJ_CONSTRUCT(&e->inbound_li, opal_list_item_t);
}

static void edge_destructor(opal_bp_graph_edge_t *e)
{
    OBJ_DESTRUCT(&e->outbound_li);
    OBJ_DESTRUCT(&e->inbound_li);
}

OBJ_CLASS_DECLARATION(opal_bp_graph_edge_t);
OBJ_CLASS_INSTANCE(opal_bp_graph_edge_t, opal_object_t, edge_constructor, edge_destructor);

static void dump_vec(const char *name, int *vec, int n) __opal_attribute_unused__;

static void dump_vec(const char *name, int *vec, int n)
{
    int i;
    fprintf(stderr, "%s={", name);
    for (i = 0; i < n; ++i) {
        fprintf(stderr, "[%d]=%2d, ", i, vec[i]);
    }
    fprintf(stderr, "}\n");
}

static void dump_vec64(const char *name, int64_t *vec, int n) __opal_attribute_unused__;

static void dump_vec64(const char *name, int64_t *vec, int n)
{
    int i;
    fprintf(stderr, "%s={", name);
    for (i = 0; i < n; ++i) {
        fprintf(stderr, "[%d]=%2" PRIi64 ", ", i, vec[i]);
    }
    fprintf(stderr, "}\n");
}

static void dump_flow(int *flow, int n) __opal_attribute_unused__;

static void dump_flow(int *flow, int n)
{
    int u, v;

    fprintf(stderr, "flow={\n");
    for (u = 0; u < n; ++u) {
        fprintf(stderr, "u=%d| ", u);
        for (v = 0; v < n; ++v) {
            fprintf(stderr, "%2d,", f(u, v));
        }
        fprintf(stderr, "\n");
    }
    fprintf(stderr, "}\n");
}

static int get_capacity(opal_bp_graph_t *g, int source, int target)
{
    opal_bp_graph_edge_t *e;

    CHECK_VERTEX_RANGE(g, source);
    CHECK_VERTEX_RANGE(g, target);

    FOREACH_OUT_EDGE(g, source, e)
    {
        assert(e->source == source);
        if (e->target == target) {
            return e->capacity;
        }
    }

    return 0;
}

static int set_capacity(opal_bp_graph_t *g, int source, int target, int cap)
{
    opal_bp_graph_edge_t *e;

    CHECK_VERTEX_RANGE(g, source);
    CHECK_VERTEX_RANGE(g, target);

    FOREACH_OUT_EDGE(g, source, e)
    {
        assert(e->source == source);
        if (e->target == target) {
            e->capacity = cap;
            return OPAL_SUCCESS;
        }
    }

    return OPAL_ERR_NOT_FOUND;
}

static void free_vertex(opal_bp_graph_t *g, opal_bp_graph_vertex_t *v)
{
    if (NULL != v) {
        if (NULL != g->v_data_cleanup_fn && NULL != v->v_data) {
            g->v_data_cleanup_fn(v->v_data);
        }
        free(v);
    }
}

int opal_bp_graph_create(opal_bp_graph_cleanup_fn_t v_data_cleanup_fn,
                         opal_bp_graph_cleanup_fn_t e_data_cleanup_fn, opal_bp_graph_t **g_out)
{
    int err;
    opal_bp_graph_t *g = NULL;

    if (NULL == g_out) {
        return OPAL_ERR_BAD_PARAM;
    }
    *g_out = NULL;

    g = calloc(1, sizeof(*g));
    if (NULL == g) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        err = OPAL_ERR_OUT_OF_RESOURCE;
        goto out_free_g;
    }

    g->source_idx = -1;
    g->sink_idx = -1;

    g->v_data_cleanup_fn = v_data_cleanup_fn;
    g->e_data_cleanup_fn = e_data_cleanup_fn;

    /* now that we essentially have an empty graph, add vertices to it */
    OBJ_CONSTRUCT(&g->vertices, opal_pointer_array_t);
    err = opal_pointer_array_init(&g->vertices, 0, INT_MAX, 32);
    if (OPAL_SUCCESS != err) {
        goto out_free_g;
    }

    *g_out = g;
    return OPAL_SUCCESS;

out_free_g:
    free(g);
    return err;
}

int opal_bp_graph_free(opal_bp_graph_t *g)
{
    int i;
    opal_bp_graph_edge_t *e, *next;
    opal_bp_graph_vertex_t *v;

    /* remove all edges from all out_edges lists */
    for (i = 0; i < NUM_VERTICES(g); ++i) {
        v = V_ID_TO_PTR(g, i);
        LIST_FOREACH_SAFE_CONTAINED(e, next, &v->out_edges, opal_bp_graph_edge_t, outbound_li)
        {
            opal_list_remove_item(&v->out_edges, &e->outbound_li);
            OBJ_RELEASE(e);
        }
    }
    /* now remove from all in_edges lists and free the edge */
    for (i = 0; i < NUM_VERTICES(g); ++i) {
        v = V_ID_TO_PTR(g, i);
        LIST_FOREACH_SAFE_CONTAINED(e, next, &v->in_edges, opal_bp_graph_edge_t, inbound_li)
        {
            opal_list_remove_item(&v->in_edges, &e->inbound_li);

            if (NULL != g->e_data_cleanup_fn && NULL != e->e_data) {
                g->e_data_cleanup_fn(e->e_data);
            }
            OBJ_RELEASE(e);
        }

        free_vertex(g, V_ID_TO_PTR(g, i));
        opal_pointer_array_set_item(&g->vertices, i, NULL);
    }
    g->num_vertices = 0;

    OBJ_DESTRUCT(&g->vertices);
    free(g);

    return OPAL_SUCCESS;
}

int opal_bp_graph_clone(const opal_bp_graph_t *g, bool copy_user_data,
                        opal_bp_graph_t **g_clone_out)
{
    int err;
    int i;
    int index;
    opal_bp_graph_t *gx;
    opal_bp_graph_edge_t *e;

    if (NULL == g_clone_out) {
        return OPAL_ERR_BAD_PARAM;
    }
    *g_clone_out = NULL;

    if (copy_user_data) {
        opal_output(0, "[%s:%d:%s] user data copy requested but not yet supported", __FILE__,
                    __LINE__, __func__);
        abort();
        return OPAL_ERR_FATAL;
    }

    gx = NULL;
    err = opal_bp_graph_create(NULL, NULL, &gx);
    if (OPAL_SUCCESS != err) {
        return err;
    }
    assert(NULL != gx);

    /* reconstruct all vertices */
    for (i = 0; i < NUM_VERTICES(g); ++i) {
        err = opal_bp_graph_add_vertex(gx, NULL, &index);
        if (OPAL_SUCCESS != err) {
            goto out_free_gx;
        }
        assert(index == i);
    }

    /* now reconstruct all the edges (iterate by source vertex only to avoid
     * double-adding) */
    for (i = 0; i < NUM_VERTICES(g); ++i) {
        FOREACH_OUT_EDGE(g, i, e)
        {
            assert(i == e->source);
            err = opal_bp_graph_add_edge(gx, e->source, e->target, e->cost, e->capacity, NULL);
            if (OPAL_SUCCESS != err) {
                goto out_free_gx;
            }
        }
    }

    *g_clone_out = gx;
    return OPAL_SUCCESS;

out_free_gx:
    /* we don't reach in and manipulate gx's state directly, so it should be
     * safe to use the standard free function */
    opal_bp_graph_free(gx);
    return err;
}

int opal_bp_graph_indegree(const opal_bp_graph_t *g, int vertex)
{
    opal_bp_graph_vertex_t *v;

    v = V_ID_TO_PTR(g, vertex);
    return opal_list_get_size(&v->in_edges);
}

int opal_bp_graph_outdegree(const opal_bp_graph_t *g, int vertex)
{
    opal_bp_graph_vertex_t *v;

    v = V_ID_TO_PTR(g, vertex);
    return opal_list_get_size(&v->out_edges);
}

int opal_bp_graph_add_edge(opal_bp_graph_t *g, int from, int to, int64_t cost, int capacity,
                           void *e_data)
{
    opal_bp_graph_edge_t *e;
    opal_bp_graph_vertex_t *v_from, *v_to;

    if (from < 0 || from >= NUM_VERTICES(g)) {
        return OPAL_ERR_BAD_PARAM;
    }
    if (to < 0 || to >= NUM_VERTICES(g)) {
        return OPAL_ERR_BAD_PARAM;
    }
    if (cost == MAX_COST) {
        return OPAL_ERR_BAD_PARAM;
    }
    if (capacity < 0) {
        /* negative cost is fine, but negative capacity is not currently
         * handled appropriately */
        return OPAL_ERR_BAD_PARAM;
    }
    FOREACH_OUT_EDGE(g, from, e)
    {
        assert(e->source == from);
        if (e->target == to) {
            return OPAL_EXISTS;
        }
    }

    /* this reference is owned by the out_edges list */
    e = OBJ_NEW(opal_bp_graph_edge_t);
    if (NULL == e) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        return OPAL_ERR_OUT_OF_RESOURCE;
    }

    e->source = from;
    e->target = to;
    e->cost = cost;
    e->capacity = capacity;
    e->e_data = e_data;

    v_from = V_ID_TO_PTR(g, from);
    opal_list_append(&v_from->out_edges, &e->outbound_li);

    OBJ_RETAIN(e); /* ref owned by in_edges list */
    v_to = V_ID_TO_PTR(g, to);
    opal_list_append(&v_to->in_edges, &e->inbound_li);

    return OPAL_SUCCESS;
}

int opal_bp_graph_add_vertex(opal_bp_graph_t *g, void *v_data, int *index_out)
{
    opal_bp_graph_vertex_t *v;

    v = calloc(1, sizeof(*v));
    if (NULL == v) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        return OPAL_ERR_OUT_OF_RESOURCE;
    }

    /* add to the ptr array early to simplify cleanup in the incredibly rare
     * chance that adding fails */
    v->v_index = opal_pointer_array_add(&g->vertices, v);
    if (-1 == v->v_index) {
        free(v);
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        return OPAL_ERR_OUT_OF_RESOURCE;
    }
    assert(v->v_index == g->num_vertices);

    ++g->num_vertices;

    v->v_data = v_data;
    OBJ_CONSTRUCT(&v->out_edges, opal_list_t);
    OBJ_CONSTRUCT(&v->in_edges, opal_list_t);

    if (NULL != index_out) {
        *index_out = v->v_index;
    }

    return OPAL_SUCCESS;
}

int opal_bp_graph_get_vertex_data(opal_bp_graph_t *g, int v_index, void **v_data_out)
{
    opal_bp_graph_vertex_t *v;

    v = V_ID_TO_PTR(g, v_index);
    if (NULL == v) {
        return OPAL_ERR_BAD_PARAM;
    }
    *v_data_out = v->v_data;

    return OPAL_SUCCESS;
}

int opal_bp_graph_order(const opal_bp_graph_t *g)
{
    return NUM_VERTICES(g);
}

/**
 * shrink a flow matrix for old_n vertices to one works for new_n
 *
 * Takes a matrix stored in a one-dimensional array of size (old_n*old_n) and
 * "truncates" it into a dense array of size (new_n*new_n) that only contain
 * the flow values for the first new_n vertices.  E.g., it turns this array
 * (old_n=5, new_n=3):
 *
 *    1  2  3  4  5
 *    6  7  8  9 10
 *   11 12 13 14 15
 *   16 17 18 19 20
 *   21 22 23 24 25
 *
 * into this array;
 *
 *    1  2  3
 *    6  7  8
 *   11 12 13
 */
static void shrink_flow_matrix(int *flow, int old_n, int new_n)
{
    int u, v;

    assert(old_n > new_n);

    for (u = 0; u < new_n; ++u) {
        for (v = 0; v < new_n; ++v) {
            flow[new_n * u + v] = flow[old_n * u + v];
        }
    }
}

/**
 * Compute the so-called "bottleneck" capacity value for a path "pred" through
 * graph "gx".
 */
static int bottleneck_path(opal_bp_graph_t *gx, int n, int *pred)
{
    int u, v;
    int min;

    min = INT_MAX;
    FOREACH_UV_ON_PATH(pred, gx->source_idx, gx->sink_idx, u, v)
    {
        int cap_f_uv = get_capacity(gx, u, v);
        min = MIN(min, cap_f_uv);
    }

    return min;
}

/**
 * This routine implements the Bellman-Ford shortest paths algorithm, slightly
 * specialized for our forumlation of flow networks:
 * https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm
 *
 * Specifically, it attempts to find the shortest path from "source" to
 * "target".  It returns true if such a path was found, false otherwise.  Any
 * found path is returned in "pred" as a predecessor chain (i.e., pred[sink]
 * is the start of the path and pred[pred[sink]] is its predecessor, etc.).
 *
 * The contents of "pred" are only valid if this routine returns true.
 */
bool opal_bp_graph_bellman_ford(opal_bp_graph_t *gx, int source, int target, int *pred)
{
    int64_t *dist;
    int i;
    int n;
    int u, v;
    bool found_target = false;

    if (NULL == gx) {
        OPAL_ERROR_LOG(OPAL_ERR_BAD_PARAM);
        return false;
    }
    if (NULL == pred) {
        OPAL_ERROR_LOG(OPAL_ERR_BAD_PARAM);
        return false;
    }
    if (source < 0 || source >= NUM_VERTICES(gx)) {
        return OPAL_ERR_BAD_PARAM;
    }
    if (target < 0 || target >= NUM_VERTICES(gx)) {
        return OPAL_ERR_BAD_PARAM;
    }

    /* initialize */
    n = opal_bp_graph_order(gx);
    dist = malloc(n * sizeof(*dist));
    if (NULL == dist) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        goto out;
    }
    for (i = 0; i < n; ++i) {
        dist[i] = MAX_COST;
        pred[i] = -1;
    }
    dist[source] = 0;

    /* relax repeatedly */
    for (i = 1; i < NUM_VERTICES(gx); ++i) {
        bool relaxed = false;
#if GRAPH_DEBUG
        dump_vec("pred", pred, NUM_VERTICES(gx));
        dump_vec64("dist", dist, NUM_VERTICES(gx));
#endif

        for (u = 0; u < NUM_VERTICES(gx); ++u) {
            opal_bp_graph_edge_t *e_ptr;

            FOREACH_OUT_EDGE(gx, u, e_ptr)
            {
                v = e_ptr->target;

                /* make sure to only construct paths from edges that actually have
                 * non-zero capacity */
                if (e_ptr->capacity > 0
                    && dist[u] != MAX_COST) { /* avoid signed overflow for "infinity" */
                    check_add64_overflow(dist[u], e_ptr->cost);
                    if ((dist[u] + e_ptr->cost) < dist[v]) {
                        dist[v] = dist[u] + e_ptr->cost;
                        pred[v] = u;
                        relaxed = true;
                    }
                }
            }
        }
        /* optimization: stop if an outer iteration did not succeed in
         * changing any dist/pred values (already at optimum) */
        if (!relaxed) {
            GRAPH_DEBUG_OUT(("relaxed==false, breaking out"));
            break;
        }
    }

    /* check for negative-cost cycles */
    for (u = 0; u < NUM_VERTICES(gx); ++u) {
        opal_bp_graph_edge_t *e_ptr;

        FOREACH_OUT_EDGE(gx, u, e_ptr)
        {
            v = e_ptr->target;
            if (e_ptr->capacity > 0 && dist[u] != MAX_COST && /* avoid signed overflow */
                (dist[u] + e_ptr->cost) < dist[v]) {
                opal_output(0, "[%s:%d:%s] negative-weight cycle detected", __FILE__, __LINE__,
                            __func__);
                abort();
                goto out;
            }
        }
    }

    if (dist[target] != MAX_COST) {
        found_target = true;
    }

out:
#if GRAPH_DEBUG
    dump_vec("pred", pred, NUM_VERTICES(gx));
#endif
    assert(pred[source] == -1);
    free(dist);
    GRAPH_DEBUG_OUT(("bellman_ford: found_target=%s", found_target ? "true" : "false"));
    return found_target;
}

/**
 * Transform the given connected, bipartite, acyclic digraph into a flow
 * network (i.e., add a source and a sink, with the source connected to vertex
 * set V1 and the sink connected to vertex set V2).  This also creates
 * residual edges suitable for augmenting-path algorithms.  All "source" nodes
 * in the original graph are considered to have an output of 1 and "sink"
 * nodes can take an input of 1.  The result is that "forward" edges are all
 * created with capacity=1, "backward" (residual) edges are created with
 * capacity=0.
 *
 * After this routine, all capacities are "residual capacities" ($c_f$ in the
 * literature).
 *
 * Initial flow throughout the network is assumed to be 0 at all edges.
 *
 * The graph will be left in an undefined state if an error occurs (though
 * freeing it should still be safe).
 */
int opal_bp_graph_bipartite_to_flow(opal_bp_graph_t *g)
{
    int err;
    int order;
    int u, v;
    int num_left, num_right;

    /* grab size before adding extra vertices */
    order = opal_bp_graph_order(g);

    err = opal_bp_graph_add_vertex(g, NULL, &g->source_idx);
    if (OPAL_SUCCESS != err) {
        return err;
    }
    err = opal_bp_graph_add_vertex(g, NULL, &g->sink_idx);
    if (OPAL_SUCCESS != err) {
        return err;
    }

    /* The networks we are interested in are bipartite and have edges only
     * from one partition to the other partition (none vice versa).  We
     * visualize this conventionally with all of the source vertices on the
     * left-hand side of an imaginary rendering of the graph and the target
     * vertices on the right-hand side of the rendering.  The direction
     * "forward" is considered to be moving from left to right.
     */
    num_left = 0;
    num_right = 0;
    for (u = 0; u < order; ++u) {
        int inbound = opal_bp_graph_indegree(g, u);
        int outbound = opal_bp_graph_outdegree(g, u);

        if (inbound > 0 && outbound > 0) {
            opal_output(0, "[%s:%d:%s] graph is not (unidirectionally) bipartite", __FILE__,
                        __LINE__, __func__);
            abort();
        } else if (inbound > 0) {
            /* "right" side of the graph, create edges to the sink */
            ++num_right;
            err = opal_bp_graph_add_edge(g, u, g->sink_idx, 0, /* no cost */
                                         /*capacity=*/1,
                                         /*e_data=*/NULL);
            if (OPAL_SUCCESS != err) {
                GRAPH_DEBUG_OUT(("add_edge failed"));
                return err;
            }
        } else if (outbound > 0) {
            /* "left" side of the graph, create edges to the source */
            ++num_left;
            err = opal_bp_graph_add_edge(g, g->source_idx, u, 0, /* no cost */
                                         /*capacity=*/1,
                                         /*e_data=*/NULL);
            if (OPAL_SUCCESS != err) {
                GRAPH_DEBUG_OUT(("add_edge failed"));
                return err;
            }
        }
    }

    /* it doesn't make sense to extend this graph with a source and sink
     * unless */
    if (num_right == 0 || num_left == 0) {
        return OPAL_ERR_BAD_PARAM;
    }

    /* now run through and create "residual" edges as well (i.e., create edges
     * in the reverse direction with 0 initial flow and a residual capacity of
     * $c_f(u,v)=c(u,v)-f(u,v)$).  Residual edges can exist where no edges
     * exist in the original graph.
     */
    order = opal_bp_graph_order(g); /* need residuals for newly created
                                          source/sink edges too */
    for (u = 0; u < order; ++u) {
        opal_bp_graph_edge_t *e_ptr;
        FOREACH_OUT_EDGE(g, u, e_ptr)
        {
            v = e_ptr->target;

            /* (u,v) exists, add (v,u) if not already present.  Cost is
             * negative for these edges because "giving back" flow pays us
             * back any cost already incurred. */
            err = opal_bp_graph_add_edge(g, v, u, -e_ptr->cost,
                                         /*capacity=*/0,
                                         /*e_data=*/NULL);
            if (OPAL_SUCCESS != err && OPAL_EXISTS != err) {
                return err;
            }
        }
    }

    return OPAL_SUCCESS;
}

/**
 * Implements the "Successive Shortest Path" algorithm for computing the
 * minimum cost flow problem.  This is a generalized version of the
 * Ford-Fulkerson algorithm.  There are two major changes from F-F:
 *  1. In addition to capacities and flows, this algorithm pays attention to
 *     costs for traversing an edge.  This particular function leaves the
 *     caller's costs alone but sets its own capacities.
 *  2. Shortest paths are computed using the cost metric.
 *
 * The algorithm's sketch looks like:
 *  1    Transform network G by adding source and sink, create residual edges
 *  2    Initial flow x is zero
 *  3    while ( Gx contains a path from s to t ) do
 *  4        Find any shortest path P from s to t
 *  5        Augment current flow x along P
 *  6        update Gx
 *
 * This function mutates the given graph (adding vertices and edges, changing
 * capacties, etc.), so callers may wish to clone the graph before calling
 * this routine.
 *
 * The result is an array of (u,v) vertex pairs, where (u,v) is an edge in the
 * original graph which has non-zero flow.
 *
 * Returns OMPI error codes like OPAL_SUCCESS/OPAL_ERR_OUT_OF_RESOURCE.
 *
 * This version of the algorithm has a theoretical upper bound on its running
 * time of O(|V|^2 * |E| * f), where f is essentially the maximum flow in the
 * graph.  In our case, f=min(|V1|,|V2|), where V1 and V2 are the two
 * constituent sets of the bipartite graph.
 *
 * This algorithm's performance could probably be improved by modifying it to
 * use vertex potentials and Dijkstra's Algorithm instead of Bellman-Ford.
 * Normally vertex potentials are needed in order to use Dijkstra's safely,
 * but our graphs are constrained enough that this may not be necessary.
 * Switching to Dijkstra's implemented with a heap should yield a reduced
 * upper bound of O(|V| * |E| * f * log(|V|)).  Let's consider this a future
 * enhancement for the time being, since it's not obvious at this point that
 * the faster running time will be worth the additional implementation
 * complexity.
 */
static int min_cost_flow_ssp(opal_bp_graph_t *gx, int **flow_out)
{
    int err = OPAL_SUCCESS;
    int n;
    int *pred = NULL;
    int *flow = NULL;
    int u, v;
    int c;

    GRAPH_DEBUG_OUT(("begin min_cost_flow_ssp()"));

    if (NULL == flow_out) {
        return OPAL_ERR_BAD_PARAM;
    }
    *flow_out = NULL;

    n = opal_bp_graph_order(gx);

    pred = malloc(n * sizeof(*pred));
    if (NULL == pred) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        err = OPAL_ERR_OUT_OF_RESOURCE;
        goto out_error;
    }

    /* "flow" is a 2d matrix of current flow values, all initialized to zero */
    flow = calloc(n * n, sizeof(*flow));
    if (NULL == flow) {
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        err = OPAL_ERR_OUT_OF_RESOURCE;
        goto out_error;
    }

    /* loop as long as paths exist from source to sink */
    while (opal_bp_graph_bellman_ford(gx, gx->source_idx, gx->sink_idx, pred)) {
        int cap_f_path;

        /* find any shortest path P from s to t (already present in pred) */
        GRAPH_DEBUG_OUT(("start outer iteration of SSP algorithm"));
#if GRAPH_DEBUG
        dump_vec("pred", pred, NUM_VERTICES(gx));
        dump_flow(flow, n);
#endif

        cap_f_path = bottleneck_path(gx, n, pred);

        /* augment current flow along P */
        FOREACH_UV_ON_PATH(pred, gx->source_idx, gx->sink_idx, u, v)
        {
            assert(u == pred[v]);

            f(u, v) = f(u, v) + cap_f_path; /* "forward" edge */
            f(v, u) = f(v, u) - cap_f_path; /* residual network edge */

            assert(f(u, v) == -f(v, u)); /* skew symmetry invariant */

            /* update Gx as we go along: decrease capacity by this new
             * augmenting flow */
            c = get_capacity(gx, u, v) - cap_f_path;
            assert(c >= 0);
            err = set_capacity(gx, u, v, c);
            if (OPAL_SUCCESS != err) {
                opal_output(0, "[%s:%d:%s] unable to set capacity, missing edge?", __FILE__,
                            __LINE__, __func__);
                abort();
            }

            c = get_capacity(gx, v, u) + cap_f_path;
            assert(c >= 0);
            err = set_capacity(gx, v, u, c);
            if (OPAL_SUCCESS != err) {
                opal_output(0, "[%s:%d:%s] unable to set capacity, missing edge?", __FILE__,
                            __LINE__, __func__);
                abort();
            }
        }
    }

out:
    *flow_out = flow;
    free(pred);
    return err;

out_error:
    free(*flow_out);
    GRAPH_DEBUG_OUT(("returning error %d", err));
    goto out;
}

int opal_bp_graph_solve_bipartite_assignment(const opal_bp_graph_t *g, int *num_match_edges_out,
                                             int **match_edges_out)
{
    int err;
    int i;
    int u, v;
    int n;
    int *flow = NULL;
    opal_bp_graph_t *gx = NULL;

    if (NULL == match_edges_out || NULL == num_match_edges_out) {
        return OPAL_ERR_BAD_PARAM;
    }
    *num_match_edges_out = 0;
    *match_edges_out = NULL;

    /* don't perturb the caller's data structure */
    err = opal_bp_graph_clone(g, false, &gx);
    if (OPAL_SUCCESS != err) {
        GRAPH_DEBUG_OUT(("opal_bp_graph_clone failed"));
        goto out;
    }

    /* Transform gx into a residual flow network with capacities, a source, a
     * sink, and residual edges.  We track the actual flow separately in the
     * "flow" matrix.  Initial capacity for every forward edge is 1.  Initial
     * capacity for every backward (residual) edge is 0.
     *
     * For the remainder of this routine (and the ssp routine) the capacities
     * refer to residual capacities ($c_f$) not capacities in the original
     * graph.  For convenience we adjust all residual capacities as we go
     * along rather than recomputing them from the flow and capacities in the
     * original graph.  This allows many other graph operations to have no
     * direct knowledge of the flow matrix.
     */
    err = opal_bp_graph_bipartite_to_flow(gx);
    if (OPAL_SUCCESS != err) {
        GRAPH_DEBUG_OUT(("bipartite_to_flow failed"));
        OPAL_ERROR_LOG(err);
        return err;
    }

    /* Use the SSP algorithm to compute the min-cost flow over this network.
     * Edges with non-zero flow in the result should be part of the matching.
     *
     * Note that the flow array returned is sized for gx, not for g.  Index
     * accordingly later on.
     */
    err = min_cost_flow_ssp(gx, &flow);
    if (OPAL_SUCCESS != err) {
        GRAPH_DEBUG_OUT(("min_cost_flow_ssp failed"));
        return err;
    }
    assert(NULL != flow);

    /* don't care about new edges in gx, only old edges in g */
    n = opal_bp_graph_order(g);

#if GRAPH_DEBUG
    dump_flow(flow, NUM_VERTICES(gx));
#endif
    shrink_flow_matrix(flow, opal_bp_graph_order(gx), n);
#if GRAPH_DEBUG
    dump_flow(flow, n);
#endif

    for (u = 0; u < n; ++u) {
        for (v = 0; v < n; ++v) {
            if (f(u, v) > 0) {
                ++(*num_match_edges_out);
            }
        }
    }

    if (0 == *num_match_edges_out) {
        /* avoid attempting to allocate a zero-byte buffer */
        goto out;
    }

    *match_edges_out = malloc(*num_match_edges_out * 2 * sizeof(int));
    if (NULL == *match_edges_out) {
        *num_match_edges_out = 0;
        OPAL_ERROR_LOG(OPAL_ERR_OUT_OF_RESOURCE);
        err = OPAL_ERR_OUT_OF_RESOURCE;
        goto out;
    }

    i = 0;
    for (u = 0; u < n; ++u) {
        for (v = 0; v < n; ++v) {
            /* flow exists on this edge so include this edge in the matching */
            if (f(u, v) > 0) {
                (*match_edges_out)[i++] = u;
                (*match_edges_out)[i++] = v;
            }
        }
    }

out:
    free(flow);
    opal_bp_graph_free(gx);
    return err;
}