Do not count an RM as R1, too.
[libfirm] / optimal.c
index 56d165b..7f17182 100644 (file)
--- a/optimal.c
+++ b/optimal.c
 #include "timing.h"
 
 pbqp_edge **edge_bucket;
+pbqp_edge **rm_bucket;
 pbqp_node **node_buckets[4];
 pbqp_node **reduced_bucket = NULL;
 static int         buckets_filled = 0;
 
-#if KAPS_STATISTIC
-static int dump = 0;
-#endif
-
 static void insert_into_edge_bucket(pbqp_edge *edge)
 {
        if (edge_bucket_contains(edge_bucket, edge)) {
@@ -65,11 +62,22 @@ static void insert_into_edge_bucket(pbqp_edge *edge)
        edge_bucket_insert(&edge_bucket, edge);
 }
 
+static void insert_into_rm_bucket(pbqp_edge *edge)
+{
+       if (edge_bucket_contains(rm_bucket, edge)) {
+               /* Edge is already inserted. */
+               return;
+       }
+
+       edge_bucket_insert(&rm_bucket, edge);
+}
+
 static void init_buckets(void)
 {
        int i;
 
        edge_bucket_init(&edge_bucket);
+       edge_bucket_init(&rm_bucket);
        node_bucket_init(&reduced_bucket);
 
        for (i = 0; i < 4; ++i) {
@@ -86,6 +94,7 @@ void free_buckets(void)
        }
 
        edge_bucket_free(&edge_bucket);
+       edge_bucket_free(&rm_bucket);
        node_bucket_free(&reduced_bucket);
 
        buckets_filled = 0;
@@ -100,8 +109,8 @@ void fill_node_buckets(pbqp *pbqp)
        node_len = pbqp->num_nodes;
 
        #if KAPS_TIMING
-               ir_timer_t *t_fill_buckets = ir_timer_register("be_pbqp_fill_buckets", "PBQP Fill Nodes into buckets");
-               ir_timer_reset_and_start(t_fill_buckets);
+               ir_timer_t *t_fill_buckets = ir_timer_new();
+               ir_timer_start(t_fill_buckets);
        #endif
 
        for (node_index = 0; node_index < node_len; ++node_index) {
@@ -124,7 +133,7 @@ void fill_node_buckets(pbqp *pbqp)
 
        #if KAPS_TIMING
                ir_timer_stop(t_fill_buckets);
-               printf("%-20s: %8.3lf msec\n", ir_timer_get_description(t_fill_buckets), (double)ir_timer_elapsed_usec(t_fill_buckets) / 1000.0);
+               printf("PBQP Fill Nodes into buckets: %8.3lf msec\n", (double)ir_timer_elapsed_usec(t_fill_buckets) / 1000.0);
        #endif
 }
 
@@ -135,9 +144,10 @@ static void normalize_towards_source(pbqp_edge *edge)
        pbqp_node      *tgt_node;
        vector         *src_vec;
        vector         *tgt_vec;
-       int             src_len;
-       int             tgt_len;
-       int             src_index;
+       unsigned        src_len;
+       unsigned        tgt_len;
+       unsigned        src_index;
+       unsigned        new_infinity = 0;
 
        assert(edge);
 
@@ -166,22 +176,28 @@ static void normalize_towards_source(pbqp_edge *edge)
                if (min != 0) {
                        if (src_vec->entries[src_index].data == INF_COSTS) {
                                pbqp_matrix_set_row_value(mat, src_index, 0);
-                       } else {
-                               pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
+                               continue;
                        }
+
+                       pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
                        src_vec->entries[src_index].data = pbqp_add(
                                        src_vec->entries[src_index].data, min);
 
                        if (min == INF_COSTS) {
-                               unsigned edge_index;
-                               unsigned edge_len = pbqp_node_get_degree(src_node);
-
-                               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
-                                       pbqp_edge *edge_candidate = src_node->edges[edge_index];
-                                       if (edge_candidate != edge) {
-                                               insert_into_edge_bucket(edge_candidate);
-                                       }
-                               }
+                               new_infinity = 1;
+                       }
+               }
+       }
+
+       if (new_infinity) {
+               unsigned edge_index;
+               unsigned edge_len = pbqp_node_get_degree(src_node);
+
+               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+                       pbqp_edge *edge_candidate = src_node->edges[edge_index];
+
+                       if (edge_candidate != edge) {
+                               insert_into_edge_bucket(edge_candidate);
                        }
                }
        }
@@ -194,9 +210,10 @@ static void normalize_towards_target(pbqp_edge *edge)
        pbqp_node      *tgt_node;
        vector         *src_vec;
        vector         *tgt_vec;
-       int             src_len;
-       int             tgt_len;
-       int             tgt_index;
+       unsigned        src_len;
+       unsigned        tgt_len;
+       unsigned        tgt_index;
+       unsigned        new_infinity = 0;
 
        assert(edge);
 
@@ -218,69 +235,222 @@ static void normalize_towards_target(pbqp_edge *edge)
        mat = edge->costs;
        assert(mat);
 
+       /* Normalize towards target node. */
        for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
                num min = pbqp_matrix_get_col_min(mat, tgt_index, src_vec);
 
                if (min != 0) {
                        if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
                                pbqp_matrix_set_col_value(mat, tgt_index, 0);
-                       } else {
-                               pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
+                               continue;
                        }
+
+                       pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
                        tgt_vec->entries[tgt_index].data = pbqp_add(
                                        tgt_vec->entries[tgt_index].data, min);
 
                        if (min == INF_COSTS) {
-                               unsigned edge_index;
-                               unsigned edge_len = pbqp_node_get_degree(tgt_node);
-
-                               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
-                                       pbqp_edge *edge_candidate = tgt_node->edges[edge_index];
-                                       if (edge_candidate != edge) {
-                                               insert_into_edge_bucket(edge_candidate);
-                                       }
-                               }
+                               new_infinity = 1;
+                       }
+               }
+       }
+
+       if (new_infinity) {
+               unsigned edge_index;
+               unsigned edge_len = pbqp_node_get_degree(tgt_node);
+
+               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+                       pbqp_edge *edge_candidate = tgt_node->edges[edge_index];
+
+                       if (edge_candidate != edge) {
+                               insert_into_edge_bucket(edge_candidate);
                        }
                }
        }
 }
 
-static void reorder_node(pbqp_node *node)
+/**
+ * Tries to apply RM for the source node of the given edge.
+ *
+ * Checks whether the source node of edge can be merged into the target node of
+ * edge, and performs the merge, if possible.
+ */
+static void merge_source_into_target(pbqp *pbqp, pbqp_edge *edge)
 {
-       unsigned    degree     = pbqp_node_get_degree(node);
-       /* Assume node lost one incident edge. */
-       unsigned    old_degree = degree + 1;
+       pbqp_matrix    *mat;
+       pbqp_node      *src_node;
+       pbqp_node      *tgt_node;
+       vector         *src_vec;
+       vector         *tgt_vec;
+       unsigned       *mapping;
+       unsigned        src_len;
+       unsigned        tgt_len;
+       unsigned        src_index;
+       unsigned        tgt_index;
+       unsigned        edge_index;
+       unsigned        edge_len;
 
-       if (!buckets_filled) return;
+       assert(pbqp);
+       assert(edge);
 
-       /* Same bucket as before */
-       if (degree > 2) return;
+       src_node = edge->src;
+       tgt_node = edge->tgt;
+       assert(src_node);
+       assert(tgt_node);
 
-       if (!node_bucket_contains(node_buckets[old_degree], node)) {
-               /* Old arity is new arity, so we have nothing to do. */
-               assert(node_bucket_contains(node_buckets[degree], node));
-               return;
+       src_vec = src_node->costs;
+       tgt_vec = tgt_node->costs;
+       assert(src_vec);
+       assert(tgt_vec);
+
+       src_len = src_vec->len;
+       tgt_len = tgt_vec->len;
+
+       /* Matrizes are normalized. */
+       assert(src_len > 1);
+       assert(tgt_len > 1);
+
+       mat = edge->costs;
+       assert(mat);
+
+       mapping = NEW_ARR_F(unsigned, tgt_len);
+
+       /* Check that each column has at most one zero entry. */
+       for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
+               unsigned onlyOneZero = 0;
+
+               if (tgt_vec->entries[tgt_index].data == INF_COSTS)
+                       continue;
+
+               for (src_index = 0; src_index < src_len; ++src_index) {
+                       if (src_vec->entries[src_index].data == INF_COSTS)
+                               continue;
+
+                       if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS)
+                               continue;
+
+                       /* Matrix entry is finite. */
+                       if (onlyOneZero) {
+                               DEL_ARR_F(mapping);
+                               return;
+                       }
+
+                       onlyOneZero = 1;
+                       mapping[tgt_index] = src_index;
+               }
        }
 
-       /* Delete node from old bucket... */
-       node_bucket_remove(&node_buckets[old_degree], node);
+       /* We know that we can merge the source node into the target node. */
+       edge_len = pbqp_node_get_degree(src_node);
 
-       /* ..and add to new one. */
-       node_bucket_insert(&node_buckets[degree], node);
+#if KAPS_STATISTIC
+       pbqp->num_rm++;
+#endif
+
+       /* Reconnect the source's edges with the target node. */
+       for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+               pbqp_edge   *old_edge = src_node->edges[edge_index];
+               pbqp_edge   *new_edge;
+               pbqp_matrix *old_matrix;
+               pbqp_matrix *new_matrix;
+               pbqp_node   *other_node;
+               vector      *other_vec;
+               unsigned     other_len;
+               unsigned     other_index;
+               unsigned     tgt_index;
+
+               assert(old_edge);
+
+               if (old_edge == edge)
+                       continue;
+
+               old_matrix = old_edge->costs;
+               assert(old_matrix);
+
+               if (old_edge->tgt == src_node) {
+                       other_node = old_edge->src;
+                       other_len  = old_matrix->rows;
+               }
+               else {
+                       other_node = old_edge->tgt;
+                       other_len = old_matrix->cols;
+               }
+               assert(other_node);
+               other_vec = other_node->costs;
+
+               new_matrix = pbqp_matrix_alloc(pbqp, tgt_len, other_len);
+
+               /* Source node selects the column of the old_matrix. */
+               if (old_edge->tgt == src_node) {
+                       for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
+                               unsigned src_index = mapping[tgt_index];
+
+                               if (tgt_vec->entries[tgt_index].data == INF_COSTS)
+                                       continue;
+
+                               for (other_index = 0; other_index < other_len; ++other_index) {
+                                       if (other_vec->entries[other_index].data == INF_COSTS)
+                                               continue;
+
+                                       new_matrix->entries[tgt_index*other_len+other_index] = old_matrix->entries[other_index*src_len+src_index];
+                               }
+                       }
+               }
+               /* Source node selects the row of the old_matrix. */
+               else {
+                       for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
+                               unsigned src_index = mapping[tgt_index];
+
+                               if (tgt_vec->entries[tgt_index].data == INF_COSTS)
+                                       continue;
+
+                               for (other_index = 0; other_index < other_len; ++other_index) {
+                                       if (other_vec->entries[other_index].data == INF_COSTS)
+                                               continue;
+
+                                       new_matrix->entries[tgt_index*other_len+other_index] = old_matrix->entries[src_index*other_len+other_index];
+                               }
+                       }
+               }
+
+               add_edge_costs(pbqp, tgt_node->index, other_node->index, new_matrix);
+
+               delete_edge(old_edge);
+               reorder_node(src_node);
+               reorder_node(other_node);
+
+               new_edge = get_edge(pbqp, tgt_node->index, other_node->index);
+               insert_into_rm_bucket(new_edge);
+       }
+
+       /* Reduce the remaining source node via RI. */
+       apply_RI(pbqp);
+
+#if KAPS_STATISTIC
+       pbqp->num_r1--;
+#endif
 }
 
-#if 0
-static void check_melting_possibility(pbqp *pbqp, pbqp_edge *edge)
+/**
+ * Tries to apply RM for the target node of the given edge.
+ *
+ * Checks whether the target node of edge can be merged into the source node of
+ * edge, and performs the merge, if possible.
+ */
+static void merge_target_into_source(pbqp *pbqp, pbqp_edge *edge)
 {
        pbqp_matrix    *mat;
        pbqp_node      *src_node;
        pbqp_node      *tgt_node;
        vector         *src_vec;
        vector         *tgt_vec;
-       int             src_len;
-       int             tgt_len;
-       int             src_index;
-       int             tgt_index;
+       unsigned       *mapping;
+       unsigned        src_len;
+       unsigned        tgt_len;
+       unsigned        src_index;
+       unsigned        tgt_index;
+       unsigned        edge_index;
+       unsigned        edge_len;
 
        assert(pbqp);
        assert(edge);
@@ -297,77 +467,189 @@ static void check_melting_possibility(pbqp *pbqp, pbqp_edge *edge)
 
        src_len = src_vec->len;
        tgt_len = tgt_vec->len;
-       assert(src_len > 0);
-       assert(tgt_len > 0);
+
+       /* Matrizes are normalized. */
+       assert(src_len > 1);
+       assert(tgt_len > 1);
 
        mat = edge->costs;
        assert(mat);
 
-       if (src_len == 1 && tgt_len == 1) {
-               //panic("Something is wrong");
-       }
+       mapping = NEW_ARR_F(unsigned, src_len);
 
-       int allRowsOk = 1;
+       /* Check that each row has at most one zero entry. */
        for (src_index = 0; src_index < src_len; ++src_index) {
-               int onlyOneZero = 0;
-               if (src_vec->entries[src_index].data == INF_COSTS) {
+               unsigned onlyOneZero = 0;
+
+               if (src_vec->entries[src_index].data == INF_COSTS)
                        continue;
-               }
+
                for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
-                       if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
+                       if (tgt_vec->entries[tgt_index].data == INF_COSTS)
                                continue;
-                       }
-                       if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
-                               if (onlyOneZero) {
-                                       onlyOneZero = 0;
-                                       break;
-                               } else {
-                                       onlyOneZero = 1;
-                                       continue;
-                               }
-                       }
-                       if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
+
+                       if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS)
                                continue;
+
+                       /* Matrix entry is finite. */
+                       if (onlyOneZero) {
+                               DEL_ARR_F(mapping);
+                               return;
                        }
-                       onlyOneZero = 0;
-                       break;
+
+                       onlyOneZero = 1;
+                       mapping[src_index] = tgt_index;
                }
-               allRowsOk &= onlyOneZero;
        }
 
-       int allColsOk = 1;
-       for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
-               int onlyOneZero = 0;
-               if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
+       /* We know that we can merge the target node into the source node. */
+       edge_len = pbqp_node_get_degree(tgt_node);
+
+#if KAPS_STATISTIC
+       pbqp->num_rm++;
+#endif
+
+       /* Reconnect the target's edges with the source node. */
+       for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+               pbqp_edge   *old_edge = tgt_node->edges[edge_index];
+               pbqp_edge   *new_edge;
+               pbqp_matrix *old_matrix;
+               pbqp_matrix *new_matrix;
+               pbqp_node   *other_node;
+               vector      *other_vec;
+               unsigned     other_len;
+               unsigned     other_index;
+               unsigned     src_index;
+
+               assert(old_edge);
+
+               if (old_edge == edge)
                        continue;
+
+               old_matrix = old_edge->costs;
+               assert(old_matrix);
+
+               if (old_edge->tgt == tgt_node) {
+                       other_node = old_edge->src;
+                       other_len  = old_matrix->rows;
                }
-               for (src_index = 0; src_index < src_len; ++src_index) {
-                       if (src_vec->entries[src_index].data == INF_COSTS) {
-                               continue;
-                       }
-                       if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
-                               if (onlyOneZero) {
-                                       onlyOneZero = 0;
-                                       break;
-                               } else {
-                                       onlyOneZero = 1;
+               else {
+                       other_node = old_edge->tgt;
+                       other_len = old_matrix->cols;
+               }
+               assert(other_node);
+               other_vec = other_node->costs;
+
+               new_matrix = pbqp_matrix_alloc(pbqp, src_len, other_len);
+
+               /* Target node selects the column of the old_matrix. */
+               if (old_edge->tgt == tgt_node) {
+                       for (src_index = 0; src_index < src_len; ++src_index) {
+                               unsigned tgt_index = mapping[src_index];
+
+                               if (src_vec->entries[src_index].data == INF_COSTS)
                                        continue;
+
+                               for (other_index = 0; other_index < other_len; ++other_index) {
+                                       if (other_vec->entries[other_index].data == INF_COSTS)
+                                               continue;
+
+                                       new_matrix->entries[src_index*other_len+other_index] = old_matrix->entries[other_index*tgt_len+tgt_index];
                                }
                        }
-                       if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
-                               continue;
+               }
+               /* Source node selects the row of the old_matrix. */
+               else {
+                       for (src_index = 0; src_index < src_len; ++src_index) {
+                               unsigned tgt_index = mapping[src_index];
+
+                               if (src_vec->entries[src_index].data == INF_COSTS)
+                                       continue;
+
+                               for (other_index = 0; other_index < other_len; ++other_index) {
+                                       if (other_vec->entries[other_index].data == INF_COSTS)
+                                               continue;
+
+                                       new_matrix->entries[src_index*other_len+other_index] = old_matrix->entries[tgt_index*other_len+other_index];
+                               }
                        }
-                       onlyOneZero = 0;
-                       break;
                }
-               allColsOk &= onlyOneZero;
+
+               add_edge_costs(pbqp, src_node->index, other_node->index, new_matrix);
+
+               delete_edge(old_edge);
+               reorder_node(tgt_node);
+               reorder_node(other_node);
+
+               new_edge = get_edge(pbqp, src_node->index, other_node->index);
+               insert_into_rm_bucket(new_edge);
        }
 
-       if (allRowsOk && allColsOk) {
-               panic("Hurray");
+       /* Reduce the remaining source node via RI. */
+       apply_RI(pbqp);
+
+#if KAPS_STATISTIC
+       pbqp->num_r1--;
+#endif
+}
+
+/**
+ * Merge neighbors into the given node.
+ */
+void apply_RM(pbqp *pbqp, pbqp_node *node)
+{
+       pbqp_edge **edges;
+       unsigned    edge_index;
+       unsigned    edge_len;
+
+       assert(node);
+       assert(pbqp);
+
+       edges    = node->edges;
+       edge_len = pbqp_node_get_degree(node);
+
+       /* Check all incident edges. */
+       for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+               pbqp_edge *edge = edges[edge_index];
+
+               insert_into_rm_bucket(edge);
+       }
+
+       /* ALAP: Merge neighbors into given node. */
+       while(edge_bucket_get_length(rm_bucket) > 0) {
+               pbqp_edge *edge = edge_bucket_pop(&rm_bucket);
+               assert(edge);
+
+               if (edge->src == node)
+                       merge_target_into_source(pbqp, edge);
+               else
+                       merge_source_into_target(pbqp, edge);
        }
 }
-#endif
+
+void reorder_node(pbqp_node *node)
+{
+       unsigned    degree     = pbqp_node_get_degree(node);
+       /* Assume node lost one incident edge. */
+       unsigned    old_degree = degree + 1;
+
+       if (!buckets_filled) return;
+
+       /* Same bucket as before */
+       if (degree > 2) return;
+
+       if (!node_bucket_contains(node_buckets[old_degree], node)) {
+               /* Old arity is new arity, so we have nothing to do. */
+               assert(node_bucket_contains(node_buckets[degree], node));
+               return;
+       }
+
+       /* Delete node from old bucket... */
+       node_bucket_remove(&node_buckets[old_degree], node);
+
+       /* ..and add to new one. */
+       node_bucket_insert(&node_buckets[degree], node);
+}
 
 void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
 {
@@ -437,9 +719,7 @@ void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
 #endif
 
 #if KAPS_STATISTIC
-               if (dump == 0) {
-                       pbqp->num_edges++;
-               }
+               pbqp->num_edges++;
 #endif
 
                delete_edge(edge);
@@ -456,8 +736,8 @@ void initial_simplify_edges(pbqp *pbqp)
        assert(pbqp);
 
        #if KAPS_TIMING
-               ir_timer_t *t_int_simpl = ir_timer_register("be_pbqp_init_simp", "PBQP Initial simplify edges");
-               ir_timer_reset_and_start(t_int_simpl);
+               ir_timer_t *t_int_simpl = ir_timer_new();
+               ir_timer_start(t_int_simpl);
        #endif
 
 #if    KAPS_DUMP
@@ -495,7 +775,7 @@ void initial_simplify_edges(pbqp *pbqp)
 
        #if KAPS_TIMING
                ir_timer_stop(t_int_simpl);
-               printf("%-20s: %8.3lf msec\n", ir_timer_get_description(t_int_simpl), (double)ir_timer_elapsed_usec(t_int_simpl) / 1000.0);
+               printf("PBQP Initial simplify edges: %8.3lf msec\n", (double)ir_timer_elapsed_usec(t_int_simpl) / 1000.0);
        #endif
 }
 
@@ -506,7 +786,7 @@ num determine_solution(pbqp *pbqp)
        num      solution   = 0;
 
        #if KAPS_TIMING
-               ir_timer_t *t_det_solution = ir_timer_register("be_det_solution", "PBQP Determine Solution");
+               ir_timer_t *t_det_solution = ir_timer_new();
                ir_timer_reset_and_start(t_det_solution);
        #endif
 
@@ -529,9 +809,7 @@ num determine_solution(pbqp *pbqp)
        node_len = node_bucket_get_length(node_buckets[0]);
 
 #if KAPS_STATISTIC
-       if (dump == 0) {
-               pbqp->num_r0 = node_len;
-       }
+       pbqp->num_r0 = node_len;
 #endif
 
        for (node_index = 0; node_index < node_len; ++node_index) {
@@ -563,7 +841,7 @@ num determine_solution(pbqp *pbqp)
 
        #if KAPS_TIMING
                ir_timer_stop(t_det_solution);
-               printf("%-20s: %8.3lf msec\n", ir_timer_get_description(t_det_solution), (double)ir_timer_elapsed_usec(t_det_solution) / 1000.0);
+               printf("PBQP Determine Solution: %8.3lf msec\n", (double)ir_timer_elapsed_usec(t_det_solution) / 1000.0);
        #endif
 
        return solution;
@@ -769,9 +1047,7 @@ void apply_RI(pbqp *pbqp)
        reorder_node(other_node);
 
 #if KAPS_STATISTIC
-       if (dump == 0) {
-               pbqp->num_r1++;
-       }
+       pbqp->num_r1++;
 #endif
 
        /* Add node to back propagation list. */
@@ -889,9 +1165,7 @@ void apply_RII(pbqp *pbqp)
        disconnect_edge(tgt_node, tgt_edge);
 
 #if KAPS_STATISTIC
-       if (dump == 0) {
-               pbqp->num_r2++;
-       }
+       pbqp->num_r2++;
 #endif
 
        /* Add node to back propagation list. */
@@ -921,6 +1195,123 @@ void apply_RII(pbqp *pbqp)
        simplify_edge(pbqp, edge);
 }
 
+static void select_column(pbqp_edge *edge, unsigned col_index)
+{
+       pbqp_matrix    *mat;
+       pbqp_node      *src_node;
+       pbqp_node      *tgt_node;
+       vector         *src_vec;
+       vector         *tgt_vec;
+       unsigned        src_len;
+       unsigned        tgt_len;
+       unsigned        src_index;
+       unsigned        new_infinity = 0;
+
+       assert(edge);
+
+       src_node = edge->src;
+       tgt_node = edge->tgt;
+       assert(src_node);
+       assert(tgt_node);
+
+       src_vec = src_node->costs;
+       tgt_vec = tgt_node->costs;
+       assert(src_vec);
+       assert(tgt_vec);
+
+       src_len = src_vec->len;
+       tgt_len = tgt_vec->len;
+       assert(src_len > 0);
+       assert(tgt_len > 0);
+
+       mat = edge->costs;
+       assert(mat);
+
+       for (src_index = 0; src_index < src_len; ++src_index) {
+               num elem = mat->entries[src_index * tgt_len + col_index];
+
+               if (elem != 0) {
+                       if (elem == INF_COSTS && src_vec->entries[src_index].data != INF_COSTS)
+                               new_infinity = 1;
+
+                       src_vec->entries[src_index].data = pbqp_add(
+                                       src_vec->entries[src_index].data, elem);
+               }
+       }
+
+       if (new_infinity) {
+               unsigned edge_index;
+               unsigned edge_len = pbqp_node_get_degree(src_node);
+
+               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+                       pbqp_edge *edge_candidate = src_node->edges[edge_index];
+
+                       if (edge_candidate != edge) {
+                               insert_into_edge_bucket(edge_candidate);
+                       }
+               }
+       }
+
+       delete_edge(edge);
+       reorder_node(src_node);
+       reorder_node(tgt_node);
+}
+
+static void select_row(pbqp_edge *edge, unsigned row_index)
+{
+       pbqp_matrix    *mat;
+       pbqp_node      *src_node;
+       pbqp_node      *tgt_node;
+       vector         *tgt_vec;
+       unsigned        tgt_len;
+       unsigned        tgt_index;
+       unsigned        new_infinity = 0;
+
+       assert(edge);
+
+       src_node = edge->src;
+       tgt_node = edge->tgt;
+       assert(tgt_node);
+
+       tgt_vec = tgt_node->costs;
+       assert(tgt_vec);
+
+       tgt_len = tgt_vec->len;
+       assert(tgt_len > 0);
+
+       mat = edge->costs;
+       assert(mat);
+
+       for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
+               num elem = mat->entries[row_index * tgt_len + tgt_index];
+
+               if (elem != 0) {
+                       if (elem == INF_COSTS && tgt_vec->entries[tgt_index].data != INF_COSTS)
+                               new_infinity = 1;
+
+                       tgt_vec->entries[tgt_index].data = pbqp_add(
+                                       tgt_vec->entries[tgt_index].data, elem);
+               }
+       }
+
+       if (new_infinity) {
+               unsigned edge_index;
+               unsigned edge_len = pbqp_node_get_degree(tgt_node);
+
+               for (edge_index = 0; edge_index < edge_len; ++edge_index) {
+                       pbqp_edge *edge_candidate = tgt_node->edges[edge_index];
+
+                       if (edge_candidate != edge) {
+                               insert_into_edge_bucket(edge_candidate);
+                       }
+               }
+       }
+
+       delete_edge(edge);
+       reorder_node(src_node);
+       reorder_node(tgt_node);
+}
+
 void select_alternative(pbqp_node *node, unsigned selected_index)
 {
        unsigned  edge_index;
@@ -942,9 +1333,14 @@ void select_alternative(pbqp_node *node, unsigned selected_index)
                }
        }
 
-       /* Add all incident edges to edge bucket, since they are now independent. */
+       /* Select corresponding row/column for incident edges. */
        for (edge_index = 0; edge_index < max_degree; ++edge_index) {
-               insert_into_edge_bucket(node->edges[edge_index]);
+               pbqp_edge *edge = node->edges[edge_index];
+
+               if (edge->src == node)
+                       select_row(edge, selected_index);
+               else
+                       select_column(edge, selected_index);
        }
 }
 
@@ -976,7 +1372,7 @@ unsigned get_local_minimal_alternative(pbqp *pbqp, pbqp_node *node)
        vector      *vec;
        pbqp_matrix *mat;
        unsigned     edge_index;
-       unsigned     max_degree   = 0;
+       unsigned     max_degree;
        unsigned     node_index;
        unsigned     node_len;
        unsigned     min_index    = 0;
@@ -985,8 +1381,9 @@ unsigned get_local_minimal_alternative(pbqp *pbqp, pbqp_node *node)
 
        assert(pbqp);
        assert(node);
-       node_vec = node->costs;
-       node_len = node_vec->len;
+       node_vec   = node->costs;
+       node_len   = node_vec->len;
+       max_degree = pbqp_node_get_degree(node);
 
        for (node_index = 0; node_index < node_len; ++node_index) {
                num value = node_vec->entries[node_index].data;