Bugfix: Do not add the cost vector of a R1 reduced node to pbqp matrix row/column...
[libfirm] / heuristical.c
1 #include "adt/array.h"
2 #include "assert.h"
3 #include "error.h"
4
5 #include "bucket.h"
6 #include "heuristical.h"
7 #include "html_dumper.h"
8 #include "kaps.h"
9 #include "matrix.h"
10 #include "pbqp_edge.h"
11 #include "pbqp_edge_t.h"
12 #include "pbqp_node.h"
13 #include "pbqp_node_t.h"
14 #include "vector.h"
15
16 static pbqp_edge **edge_bucket;
17 static pbqp_node **node_buckets[4];
18 static pbqp_node **reduced_bucket = NULL;
19 static int         buckets_filled = 0;
20
21 #if KAPS_STATISTIC
22 static int dump = 0;
23 #endif
24
25 /* Forward declarations. */
26 static void apply_Brute_Force(pbqp *pbqp);
27
28 static void insert_into_edge_bucket(pbqp_edge *edge)
29 {
30         if (edge_bucket_contains(edge_bucket, edge)) {
31                 /* Edge is already inserted. */
32                 return;
33         }
34
35         edge_bucket_insert(&edge_bucket, edge);
36 }
37
38 static void init_buckets(void)
39 {
40         int i;
41
42         edge_bucket_init(&edge_bucket);
43         node_bucket_init(&reduced_bucket);
44
45         for (i = 0; i < 4; ++i) {
46                 node_bucket_init(&node_buckets[i]);
47         }
48 }
49
50 static void free_buckets(void)
51 {
52         int i;
53
54         for (i = 0; i < 4; ++i) {
55                 node_bucket_free(&node_buckets[i]);
56         }
57
58         edge_bucket_free(&edge_bucket);
59         node_bucket_free(&reduced_bucket);
60
61         buckets_filled = 0;
62 }
63
64 static void fill_node_buckets(pbqp *pbqp)
65 {
66         unsigned node_index;
67         unsigned node_len;
68
69         assert(pbqp);
70         node_len = pbqp->num_nodes;
71
72         for (node_index = 0; node_index < node_len; ++node_index) {
73                 unsigned   degree;
74                 pbqp_node *node = get_node(pbqp, node_index);
75
76                 if (!node) continue;
77
78                 degree = pbqp_node_get_degree(node);
79
80                 /* We have only one bucket for nodes with arity >= 3. */
81                 if (degree > 3) {
82                         degree = 3;
83                 }
84
85                 node_bucket_insert(&node_buckets[degree], node);
86         }
87
88         buckets_filled = 1;
89 }
90
91 static void normalize_towards_source(pbqp *pbqp, pbqp_edge *edge)
92 {
93         pbqp_matrix    *mat;
94         pbqp_node      *src_node;
95         pbqp_node      *tgt_node;
96         vector         *src_vec;
97         vector         *tgt_vec;
98         int             src_len;
99         int             tgt_len;
100         int             src_index;
101
102         assert(pbqp);
103         assert(edge);
104
105         src_node = edge->src;
106         tgt_node = edge->tgt;
107         assert(src_node);
108         assert(tgt_node);
109
110         src_vec = src_node->costs;
111         tgt_vec = tgt_node->costs;
112         assert(src_vec);
113         assert(tgt_vec);
114
115         src_len = src_vec->len;
116         tgt_len = tgt_vec->len;
117         assert(src_len > 0);
118         assert(tgt_len > 0);
119
120         mat = edge->costs;
121         assert(mat);
122
123         /* Normalize towards source node. */
124         for (src_index = 0; src_index < src_len; ++src_index) {
125                 num min = pbqp_matrix_get_row_min(mat, src_index, tgt_vec);
126
127                 if (min != 0) {
128                         if (src_vec->entries[src_index].data == INF_COSTS) {
129                                 pbqp_matrix_set_row_value(mat, src_index, 0);
130                         } else {
131                                 pbqp_matrix_sub_row_value(mat, src_index, tgt_vec, min);
132                         }
133                         src_vec->entries[src_index].data = pbqp_add(
134                                         src_vec->entries[src_index].data, min);
135
136                         if (min == INF_COSTS) {
137                                 unsigned edge_index;
138                                 unsigned edge_len = pbqp_node_get_degree(src_node);
139
140                                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
141                                         pbqp_edge *edge_candidate = src_node->edges[edge_index];
142                                         if (edge_candidate != edge) {
143                                                 insert_into_edge_bucket(edge_candidate);
144                                         }
145                                 }
146                         }
147                 }
148         }
149 }
150
151 static void normalize_towards_target(pbqp *pbqp, pbqp_edge *edge)
152 {
153         pbqp_matrix    *mat;
154         pbqp_node      *src_node;
155         pbqp_node      *tgt_node;
156         vector         *src_vec;
157         vector         *tgt_vec;
158         int             src_len;
159         int             tgt_len;
160         int             tgt_index;
161
162         assert(pbqp);
163         assert(edge);
164
165         src_node = edge->src;
166         tgt_node = edge->tgt;
167         assert(src_node);
168         assert(tgt_node);
169
170         src_vec = src_node->costs;
171         tgt_vec = tgt_node->costs;
172         assert(src_vec);
173         assert(tgt_vec);
174
175         src_len = src_vec->len;
176         tgt_len = tgt_vec->len;
177         assert(src_len > 0);
178         assert(tgt_len > 0);
179
180         mat = edge->costs;
181         assert(mat);
182
183         for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
184                 num min = pbqp_matrix_get_col_min(mat, tgt_index, src_vec);
185
186                 if (min != 0) {
187                         if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
188                                 pbqp_matrix_set_col_value(mat, tgt_index, 0);
189                         } else {
190                                 pbqp_matrix_sub_col_value(mat, tgt_index, src_vec, min);
191                         }
192                         tgt_vec->entries[tgt_index].data = pbqp_add(
193                                         tgt_vec->entries[tgt_index].data, min);
194
195                         if (min == INF_COSTS) {
196                                 unsigned edge_index;
197                                 unsigned edge_len = pbqp_node_get_degree(tgt_node);
198
199                                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
200                                         pbqp_edge *edge_candidate = tgt_node->edges[edge_index];
201                                         if (edge_candidate != edge) {
202                                                 insert_into_edge_bucket(edge_candidate);
203                                         }
204                                 }
205                         }
206                 }
207         }
208 }
209
210 static void reorder_node(pbqp_node *node)
211 {
212         unsigned    degree     = pbqp_node_get_degree(node);
213         /* Assume node lost one incident edge. */
214         unsigned    old_degree = degree + 1;
215
216         if (!buckets_filled) return;
217
218         /* Same bucket as before */
219         if (degree > 2) return;
220
221         if (!node_bucket_contains(node_buckets[old_degree], node)) {
222                 /* Old arity is new arity, so we have nothing to do. */
223                 assert(node_bucket_contains(node_buckets[degree], node));
224                 return;
225         }
226
227         /* Delete node from old bucket... */
228         node_bucket_remove(&node_buckets[old_degree], node);
229
230         /* ..and add to new one. */
231         node_bucket_insert(&node_buckets[degree], node);
232 }
233
234 static void check_melting_possibility(pbqp *pbqp, pbqp_edge *edge)
235 {
236         pbqp_matrix    *mat;
237         pbqp_node      *src_node;
238         pbqp_node      *tgt_node;
239         vector         *src_vec;
240         vector         *tgt_vec;
241         int             src_len;
242         int             tgt_len;
243         int             src_index;
244         int             tgt_index;
245
246         assert(pbqp);
247         assert(edge);
248
249         src_node = edge->src;
250         tgt_node = edge->tgt;
251         assert(src_node);
252         assert(tgt_node);
253
254         src_vec = src_node->costs;
255         tgt_vec = tgt_node->costs;
256         assert(src_vec);
257         assert(tgt_vec);
258
259         src_len = src_vec->len;
260         tgt_len = tgt_vec->len;
261         assert(src_len > 0);
262         assert(tgt_len > 0);
263
264         mat = edge->costs;
265         assert(mat);
266
267         if (src_len == 1 && tgt_len == 1) {
268                 //panic("Something is wrong");
269         }
270
271         int allRowsOk = 1;
272         for (src_index = 0; src_index < src_len; ++src_index) {
273                 int onlyOneZero = 0;
274                 if (src_vec->entries[src_index].data == INF_COSTS) {
275                         continue;
276                 }
277                 for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
278                         if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
279                                 continue;
280                         }
281                         if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
282                                 if (onlyOneZero) {
283                                         onlyOneZero = 0;
284                                         break;
285                                 } else {
286                                         onlyOneZero = 1;
287                                         continue;
288                                 }
289                         }
290                         if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
291                                 continue;
292                         }
293                         onlyOneZero = 0;
294                         break;
295                 }
296                 allRowsOk &= onlyOneZero;
297         }
298
299         int allColsOk = 1;
300         for (tgt_index = 0; tgt_index < tgt_len; ++tgt_index) {
301                 int onlyOneZero = 0;
302                 if (tgt_vec->entries[tgt_index].data == INF_COSTS) {
303                         continue;
304                 }
305                 for (src_index = 0; src_index < src_len; ++src_index) {
306                         if (src_vec->entries[src_index].data == INF_COSTS) {
307                                 continue;
308                         }
309                         if (mat->entries[src_index * tgt_len + tgt_index] == 0) {
310                                 if (onlyOneZero) {
311                                         onlyOneZero = 0;
312                                         break;
313                                 } else {
314                                         onlyOneZero = 1;
315                                         continue;
316                                 }
317                         }
318                         if (mat->entries[src_index * tgt_len + tgt_index] == INF_COSTS) {
319                                 continue;
320                         }
321                         onlyOneZero = 0;
322                         break;
323                 }
324                 allColsOk &= onlyOneZero;
325         }
326
327         if (allRowsOk && allColsOk) {
328                 panic("Hurray");
329         }
330 }
331
332 static void simplify_edge(pbqp *pbqp, pbqp_edge *edge)
333 {
334         pbqp_matrix    *mat;
335         pbqp_node      *src_node;
336         pbqp_node      *tgt_node;
337         vector         *src_vec;
338         vector         *tgt_vec;
339         int             src_len;
340         int             tgt_len;
341
342         assert(pbqp);
343         assert(edge);
344
345         src_node = edge->src;
346         tgt_node = edge->tgt;
347         assert(src_node);
348         assert(tgt_node);
349
350         /* If edge are already deleted, we have nothing to do. */
351         if (!is_connected(src_node, edge) || !is_connected(tgt_node, edge))
352                 return;
353
354         if (pbqp->dump_file) {
355                 char txt[100];
356                 sprintf(txt, "Simplification of Edge n%d-n%d", src_node->index, tgt_node->index);
357                 dump_section(pbqp->dump_file, 3, txt);
358         }
359
360         src_vec = src_node->costs;
361         tgt_vec = tgt_node->costs;
362         assert(src_vec);
363         assert(tgt_vec);
364
365         src_len = src_vec->len;
366         tgt_len = tgt_vec->len;
367         assert(src_len > 0);
368         assert(tgt_len > 0);
369
370         mat = edge->costs;
371         assert(mat);
372
373         if (pbqp->dump_file) {
374                 fputs("Input:<br>\n", pbqp->dump_file);
375                 dump_simplifyedge(pbqp, edge);
376         }
377
378         normalize_towards_source(pbqp, edge);
379         normalize_towards_target(pbqp, edge);
380
381         if (pbqp->dump_file) {
382                 fputs("<br>\nOutput:<br>\n", pbqp->dump_file);
383                 dump_simplifyedge(pbqp, edge);
384         }
385
386         if (pbqp_matrix_is_zero(mat, src_vec, tgt_vec)) {
387                 if (pbqp->dump_file) {
388                         fputs("edge has been eliminated<br>\n", pbqp->dump_file);
389                 }
390
391 #if KAPS_STATISTIC
392                 if (dump == 0) {
393                         pbqp->num_edges++;
394                 }
395 #endif
396
397                 delete_edge(edge);
398                 reorder_node(src_node);
399                 reorder_node(tgt_node);
400         } else {
401                 //check_melting_possibility(pbqp, edge);
402         }
403 }
404
405 static void initial_simplify_edges(pbqp *pbqp)
406 {
407         unsigned node_index;
408         unsigned node_len;
409
410         assert(pbqp);
411
412         if (pbqp->dump_file) {
413                 pbqp_dump_input(pbqp);
414                 dump_section(pbqp->dump_file, 1, "2. Simplification of Cost Matrices");
415         }
416
417         node_len = pbqp->num_nodes;
418
419         init_buckets();
420
421         /* First simplify all edges. */
422         for (node_index = 0; node_index < node_len; ++node_index) {
423                 unsigned    edge_index;
424                 pbqp_node  *node = get_node(pbqp, node_index);
425                 pbqp_edge **edges;
426                 unsigned    edge_len;
427
428                 if (!node) continue;
429
430                 edges = node->edges;
431                 edge_len = pbqp_node_get_degree(node);
432
433                 for (edge_index = 0; edge_index < edge_len; ++edge_index) {
434                         pbqp_edge *edge = edges[edge_index];
435
436                         /* Simplify only once per edge. */
437                         if (node != edge->src) continue;
438
439                         simplify_edge(pbqp, edge);
440                 }
441         }
442 }
443
444 static num determine_solution(pbqp *pbqp)
445 {
446         unsigned node_index;
447         unsigned node_len;
448         num      solution   = 0;
449         FILE     *file;
450
451         assert(pbqp);
452         file = pbqp->dump_file;
453
454         if (file) {
455                 dump_section(file, 1, "4. Determine Solution/Minimum");
456                 dump_section(file, 2, "4.1. Trivial Solution");
457         }
458
459         /* Solve trivial nodes and calculate solution. */
460         node_len = node_bucket_get_length(node_buckets[0]);
461
462 #if KAPS_STATISTIC
463         if (dump == 0) {
464                 pbqp->num_r0 = node_len;
465         }
466 #endif
467
468         for (node_index = 0; node_index < node_len; ++node_index) {
469                 pbqp_node *node = node_buckets[0][node_index];
470                 assert(node);
471
472                 node->solution = vector_get_min_index(node->costs);
473                 solution       = pbqp_add(solution,
474                                 node->costs->entries[node->solution].data);
475                 if (file) {
476                         fprintf(file, "node n%d is set to %d<br>\n", node->index, node->solution);
477                         dump_node(file, node);
478                 }
479         }
480
481         if (file) {
482                 dump_section(file, 2, "Minimum");
483                 fprintf(file, "Minimum is equal to %lld.", solution);
484         }
485
486         return solution;
487 }
488
489 static void back_propagate(pbqp *pbqp)
490 {
491         unsigned node_index;
492         unsigned node_len   = node_bucket_get_length(reduced_bucket);
493
494         assert(pbqp);
495         if (pbqp->dump_file) {
496                 dump_section(pbqp->dump_file, 2, "Back Propagation");
497         }
498
499         for (node_index = node_len; node_index > 0; --node_index) {
500                 pbqp_node *node = reduced_bucket[node_index - 1];
501
502                 switch (pbqp_node_get_degree(node)) {
503                         case 1:
504                                 back_propagate_RI(pbqp, node);
505                                 break;
506                         case 2:
507                                 back_propagate_RII(pbqp, node);
508                                 break;
509                         default:
510                                 panic("Only nodes with degree one or two should be in this bucket");
511                                 break;
512                 }
513         }
514 }
515
516 static void apply_heuristic_reductions(pbqp *pbqp)
517 {
518         for (;;) {
519                 if (edge_bucket_get_length(edge_bucket) > 0) {
520                         apply_edge(pbqp);
521                 } else if (node_bucket_get_length(node_buckets[1]) > 0) {
522                         apply_RI(pbqp);
523                 } else if (node_bucket_get_length(node_buckets[2]) > 0) {
524                         apply_RII(pbqp);
525                 } else if (node_bucket_get_length(node_buckets[3]) > 0) {
526                         apply_RN(pbqp);
527                 } else {
528                         return;
529                 }
530         }
531 }
532
533 void solve_pbqp_heuristical(pbqp *pbqp)
534 {
535         /* Reduce nodes degree ... */
536         initial_simplify_edges(pbqp);
537
538         /* ... and put node into bucket representing their degree. */
539         fill_node_buckets(pbqp);
540
541 #if KAPS_STATISTIC
542         FILE *fh = fopen("solutions.pb", "a");
543         fprintf(fh, "Solution");
544         fclose(fh);
545 #endif
546
547         apply_heuristic_reductions(pbqp);
548
549         pbqp->solution = determine_solution(pbqp);
550
551 #if KAPS_STATISTIC
552         fh = fopen("solutions.pb", "a");
553         fprintf(fh, ": %lld RE:%u R0:%u R1:%u R2:%u RN/BF:%u\n", pbqp->solution,
554                                 pbqp->num_edges, pbqp->num_r0, pbqp->num_r1, pbqp->num_r2,
555                                 pbqp->num_rn);
556         fclose(fh);
557 #endif
558
559         /* Solve reduced nodes. */
560         back_propagate(pbqp);
561
562         free_buckets();
563 }
564
565 void apply_edge(pbqp *pbqp)
566 {
567         pbqp_edge *edge = edge_bucket_pop(&edge_bucket);
568
569         simplify_edge(pbqp, edge);
570 }
571
572 void apply_RI(pbqp *pbqp)
573 {
574         pbqp_node   *node       = node_bucket_pop(&node_buckets[1]);
575         pbqp_edge   *edge       = node->edges[0];
576         pbqp_matrix *mat        = edge->costs;
577         int          is_src     = edge->src == node;
578         pbqp_node   *other_node;
579
580         assert(pbqp_node_get_degree(node) == 1);
581
582         if (is_src) {
583                 other_node = edge->tgt;
584         } else {
585                 other_node = edge->src;
586         }
587
588         if (pbqp->dump_file) {
589                 char     txt[100];
590                 sprintf(txt, "RI-Reduction of Node n%d", node->index);
591                 dump_section(pbqp->dump_file, 2, txt);
592                 pbqp_dump_graph(pbqp);
593                 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
594                 dump_node(pbqp->dump_file, node);
595                 dump_node(pbqp->dump_file, other_node);
596                 dump_edge(pbqp->dump_file, edge);
597         }
598
599         if (is_src) {
600                 pbqp_matrix_add_to_all_cols(mat, node->costs);
601                 normalize_towards_target(pbqp, edge);
602         } else {
603                 pbqp_matrix_add_to_all_rows(mat, node->costs);
604                 normalize_towards_source(pbqp, edge);
605         }
606         disconnect_edge(other_node, edge);
607
608         if (pbqp->dump_file) {
609                 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
610                 dump_node(pbqp->dump_file, other_node);
611         }
612
613         reorder_node(other_node);
614
615 #if KAPS_STATISTIC
616         if (dump == 0) {
617                 pbqp->num_r1++;
618         }
619 #endif
620
621         /* Add node to back propagation list. */
622         node_bucket_insert(&reduced_bucket, node);
623 }
624
625 void apply_RII(pbqp *pbqp)
626 {
627         pbqp_node   *node       = node_bucket_pop(&node_buckets[2]);
628         pbqp_edge   *src_edge   = node->edges[0];
629         pbqp_edge   *tgt_edge   = node->edges[1];
630         int          src_is_src = src_edge->src == node;
631         int          tgt_is_src = tgt_edge->src == node;
632         pbqp_matrix *src_mat;
633         pbqp_matrix *tgt_mat;
634         pbqp_node   *src_node;
635         pbqp_node   *tgt_node;
636         pbqp_matrix *mat;
637         vector      *vec;
638         vector      *node_vec;
639         vector      *src_vec;
640         vector      *tgt_vec;
641         unsigned     col_index;
642         unsigned     col_len;
643         unsigned     row_index;
644         unsigned     row_len;
645         unsigned     node_len;
646
647         assert(pbqp);
648         assert(pbqp_node_get_degree(node) == 2);
649
650         if (src_is_src) {
651                 src_node = src_edge->tgt;
652         } else {
653                 src_node = src_edge->src;
654         }
655
656         if (tgt_is_src) {
657                 tgt_node = tgt_edge->tgt;
658         } else {
659                 tgt_node = tgt_edge->src;
660         }
661
662         /* Swap nodes if necessary. */
663         if (tgt_node->index < src_node->index) {
664                 pbqp_node *tmp_node;
665                 pbqp_edge *tmp_edge;
666
667                 tmp_node = src_node;
668                 src_node = tgt_node;
669                 tgt_node = tmp_node;
670
671                 tmp_edge = src_edge;
672                 src_edge = tgt_edge;
673                 tgt_edge = tmp_edge;
674
675                 src_is_src = src_edge->src == node;
676                 tgt_is_src = tgt_edge->src == node;
677         }
678
679         if (pbqp->dump_file) {
680                 char     txt[100];
681                 sprintf(txt, "RII-Reduction of Node n%d", node->index);
682                 dump_section(pbqp->dump_file, 2, txt);
683                 pbqp_dump_graph(pbqp);
684                 fputs("<br>\nBefore reduction:<br>\n", pbqp->dump_file);
685                 dump_node(pbqp->dump_file, src_node);
686                 dump_edge(pbqp->dump_file, src_edge);
687                 dump_node(pbqp->dump_file, node);
688                 dump_edge(pbqp->dump_file, tgt_edge);
689                 dump_node(pbqp->dump_file, tgt_node);
690         }
691
692         src_mat = src_edge->costs;
693         tgt_mat = tgt_edge->costs;
694
695         src_vec  = src_node->costs;
696         tgt_vec  = tgt_node->costs;
697         node_vec = node->costs;
698
699         row_len  = src_vec->len;
700         col_len  = tgt_vec->len;
701         node_len = node_vec->len;
702
703         mat = pbqp_matrix_alloc(pbqp, row_len, col_len);
704
705         for (row_index = 0; row_index < row_len; ++row_index) {
706                 for (col_index = 0; col_index < col_len; ++col_index) {
707                         vec = vector_copy(pbqp, node_vec);
708
709                         if (src_is_src) {
710                                 vector_add_matrix_col(vec, src_mat, row_index);
711                         } else {
712                                 vector_add_matrix_row(vec, src_mat, row_index);
713                         }
714
715                         if (tgt_is_src) {
716                                 vector_add_matrix_col(vec, tgt_mat, col_index);
717                         } else {
718                                 vector_add_matrix_row(vec, tgt_mat, col_index);
719                         }
720
721                         mat->entries[row_index * col_len + col_index] = vector_get_min(vec);
722
723                         obstack_free(&pbqp->obstack, vec);
724                 }
725         }
726
727         pbqp_edge *edge = get_edge(pbqp, src_node->index, tgt_node->index);
728
729         /* Disconnect node. */
730         disconnect_edge(src_node, src_edge);
731         disconnect_edge(tgt_node, tgt_edge);
732
733 #if KAPS_STATISTIC
734         if (dump == 0) {
735                 pbqp->num_r2++;
736         }
737 #endif
738
739         /* Add node to back propagation list. */
740         node_bucket_insert(&reduced_bucket, node);
741
742         if (edge == NULL) {
743                 edge = alloc_edge(pbqp, src_node->index, tgt_node->index, mat);
744         } else {
745                 pbqp_matrix_add(edge->costs, mat);
746
747                 /* Free local matrix. */
748                 obstack_free(&pbqp->obstack, mat);
749
750                 reorder_node(src_node);
751                 reorder_node(tgt_node);
752         }
753
754         if (pbqp->dump_file) {
755                 fputs("<br>\nAfter reduction:<br>\n", pbqp->dump_file);
756                 dump_edge(pbqp->dump_file, edge);
757         }
758
759         /* Edge has changed so we simplify it. */
760         simplify_edge(pbqp, edge);
761 }
762
763 static void select_alternative(pbqp_node *node, unsigned selected_index)
764 {
765         unsigned  edge_index;
766         unsigned  node_index;
767         unsigned  node_len;
768         vector   *node_vec;
769         unsigned  max_degree = pbqp_node_get_degree(node);
770
771         assert(node);
772         node->solution = selected_index;
773         node_vec = node->costs;
774         node_len = node_vec->len;
775         assert(selected_index < node_len);
776
777         /* Set all other costs to infinity. */
778         for (node_index = 0; node_index < node_len; ++node_index) {
779                 if (node_index != selected_index) {
780                         node_vec->entries[node_index].data = INF_COSTS;
781                 }
782         }
783
784         /* Add all incident edges to edge bucket, since they are now independent. */
785         for (edge_index = 0; edge_index < max_degree; ++edge_index) {
786                 insert_into_edge_bucket(node->edges[edge_index]);
787         }
788 }
789
790 static pbqp_node *get_node_with_max_degree(void)
791 {
792         pbqp_node  **bucket       = node_buckets[3];
793         unsigned     bucket_len   = node_bucket_get_length(bucket);
794         unsigned     bucket_index;
795         unsigned     max_degree   = 0;
796         pbqp_node   *result       = NULL;
797
798         for (bucket_index = 0; bucket_index < bucket_len; ++bucket_index) {
799                 pbqp_node *candidate = bucket[bucket_index];
800                 unsigned   degree    = pbqp_node_get_degree(candidate);
801
802                 if (degree > max_degree) {
803                         result = candidate;
804                         max_degree = degree;
805                 }
806         }
807
808         return result;
809 }
810
811 static unsigned get_local_minimal_alternative(pbqp *pbqp, pbqp_node *node)
812 {
813         pbqp_edge   *edge;
814         vector      *node_vec;
815         vector      *vec;
816         pbqp_matrix *mat;
817         unsigned     edge_index;
818         unsigned     max_degree   = 0;
819         unsigned     node_index;
820         unsigned     node_len;
821         unsigned     min_index    = 0;
822         num          min          = INF_COSTS;
823         int          is_src;
824
825         assert(pbqp);
826         assert(node);
827         node_vec = node->costs;
828         node_len = node_vec->len;
829
830         for (node_index = 0; node_index < node_len; ++node_index) {
831                 num value = node_vec->entries[node_index].data;
832
833                 for (edge_index = 0; edge_index < max_degree; ++edge_index) {
834                         edge   = node->edges[edge_index];
835                         mat    = edge->costs;
836                         is_src = edge->src == node;
837
838                         if (is_src) {
839                                 vec = vector_copy(pbqp, edge->tgt->costs);
840                                 vector_add_matrix_row(vec, mat, node_index);
841                         } else {
842                                 vec = vector_copy(pbqp, edge->src->costs);
843                                 vector_add_matrix_col(vec, mat, node_index);
844                         }
845
846                         value = pbqp_add(value, vector_get_min(vec));
847
848                         obstack_free(&pbqp->obstack, vec);
849                 }
850
851                 if (value < min) {
852                         min = value;
853                         min_index = node_index;
854                 }
855         }
856
857         return min_index;
858 }
859
860 void apply_RN(pbqp *pbqp)
861 {
862         pbqp_node   *node         = NULL;
863         unsigned     min_index    = 0;
864
865         assert(pbqp);
866
867         /* We want to reduce a node with maximum degree. */
868         node = get_node_with_max_degree();
869         assert(node);
870         assert(pbqp_node_get_degree(node) > 2);
871
872         if (pbqp->dump_file) {
873                 char     txt[100];
874                 sprintf(txt, "RN-Reduction of Node n%d", node->index);
875                 dump_section(pbqp->dump_file, 2, txt);
876                 pbqp_dump_graph(pbqp);
877         }
878
879         min_index = get_local_minimal_alternative(pbqp, node);
880
881         if (pbqp->dump_file) {
882                 fprintf(pbqp->dump_file, "node n%d is set to %d<br><br>\n",
883                                         node->index, min_index);
884         }
885
886 #if KAPS_STATISTIC
887         if (dump == 0) {
888                 FILE *fh = fopen("solutions.pb", "a");
889                 fprintf(fh, "[%u]", min_index);
890                 fclose(fh);
891                 pbqp->num_rn++;
892         }
893 #endif
894
895         /* Now that we found the local minimum set all other costs to infinity. */
896         select_alternative(node, min_index);
897 }
898
899 static void apply_brute_force_reductions(pbqp *pbqp)
900 {
901         for (;;) {
902                 if (edge_bucket_get_length(edge_bucket) > 0) {
903                         apply_edge(pbqp);
904                 } else if (node_bucket_get_length(node_buckets[1]) > 0) {
905                         apply_RI(pbqp);
906                 } else if (node_bucket_get_length(node_buckets[2]) > 0) {
907                         apply_RII(pbqp);
908                 } else if (node_bucket_get_length(node_buckets[3]) > 0) {
909                         apply_Brute_Force(pbqp);
910                 } else {
911                         return;
912                 }
913         }
914 }
915
916 static unsigned get_minimal_alternative(pbqp *pbqp, pbqp_node *node)
917 {
918         vector      *node_vec;
919         unsigned     node_index;
920         unsigned     node_len;
921         unsigned     min_index    = 0;
922         num          min          = INF_COSTS;
923         unsigned     bucket_index;
924
925         assert(pbqp);
926         assert(node);
927         node_vec     = node->costs;
928         node_len     = node_vec->len;
929         bucket_index = node->bucket_index;
930
931         for (node_index = 0; node_index < node_len; ++node_index) {
932                 pbqp_node_bucket bucket_deg3;
933                 num              value;
934                 unsigned         bucket_0_length;
935                 unsigned         bucket_red_length;
936
937                 char *tmp = obstack_finish(&pbqp->obstack);
938
939                 node_bucket_init(&bucket_deg3);
940
941                 /* Some node buckets and the edge bucket should be empty. */
942                 assert(node_bucket_get_length(node_buckets[1]) == 0);
943                 assert(node_bucket_get_length(node_buckets[2]) == 0);
944                 assert(edge_bucket_get_length(edge_bucket)     == 0);
945
946                 /* char *tmp = obstack_finish(&pbqp->obstack); */
947
948                 /* Save current PBQP state. */
949                 node_bucket_copy(&bucket_deg3, node_buckets[3]);
950                 node_bucket_shrink(&node_buckets[3], 0);
951                 node_bucket_deep_copy(pbqp, &node_buckets[3], bucket_deg3);
952                 node_bucket_update(pbqp, node_buckets[3]);
953                 bucket_0_length   = node_bucket_get_length(node_buckets[0]);
954                 bucket_red_length = node_bucket_get_length(reduced_bucket);
955
956                 /* Select alternative and solve PBQP recursively. */
957                 select_alternative(node_buckets[3][bucket_index], node_index);
958                 apply_brute_force_reductions(pbqp);
959
960                 value = determine_solution(pbqp);
961
962                 if (value < min) {
963                         min = value;
964                         min_index = node_index;
965                 }
966
967                 /* Some node buckets and the edge bucket should still be empty. */
968                 assert(node_bucket_get_length(node_buckets[1]) == 0);
969                 assert(node_bucket_get_length(node_buckets[2]) == 0);
970                 assert(edge_bucket_get_length(edge_bucket)     == 0);
971
972                 /* Clear modified buckets... */
973                 node_bucket_shrink(&node_buckets[3], 0);
974
975                 /* ... and restore old PBQP state. */
976                 node_bucket_shrink(&node_buckets[0], bucket_0_length);
977                 node_bucket_shrink(&reduced_bucket, bucket_red_length);
978                 node_bucket_copy(&node_buckets[3], bucket_deg3);
979                 node_bucket_update(pbqp, node_buckets[3]);
980
981                 /* Free copies. */
982                 /* obstack_free(&pbqp->obstack, tmp); */
983                 node_bucket_free(&bucket_deg3);
984
985                 obstack_free(&pbqp->obstack, tmp);
986         }
987
988         return min_index;
989 }
990
991 void apply_Brute_Force(pbqp *pbqp)
992 {
993         pbqp_node   *node         = NULL;
994         unsigned     min_index    = 0;
995
996         assert(pbqp);
997
998         /* We want to reduce a node with maximum degree. */
999         node = get_node_with_max_degree();
1000         assert(node);
1001         assert(pbqp_node_get_degree(node) > 2);
1002
1003         if (pbqp->dump_file) {
1004                 char     txt[100];
1005                 sprintf(txt, "BF-Reduction of Node n%d", node->index);
1006                 dump_section(pbqp->dump_file, 2, txt);
1007                 pbqp_dump_graph(pbqp);
1008         }
1009
1010 #if KAPS_STATISTIC
1011         dump++;
1012 #endif
1013
1014         min_index = get_minimal_alternative(pbqp, node);
1015         node = pbqp->nodes[node->index];
1016
1017         if (pbqp->dump_file) {
1018                 fprintf(pbqp->dump_file, "node n%d is set to %d<br><br>\n",
1019                                         node->index, min_index);
1020         }
1021
1022 #if KAPS_STATISTIC
1023         dump--;
1024         if (dump == 0) {
1025                 FILE *fh = fopen("solutions.pb", "a");
1026                 fprintf(fh, "[%u]", min_index);
1027                 fclose(fh);
1028                 pbqp->num_bf++;
1029         }
1030 #endif
1031
1032         /* Now that we found the minimum set all other costs to infinity. */
1033         select_alternative(node, min_index);
1034 }
1035
1036 void solve_pbqp_brute_force(pbqp *pbqp)
1037 {
1038         /* Reduce nodes degree ... */
1039         initial_simplify_edges(pbqp);
1040
1041         /* ... and put node into bucket representing their degree. */
1042         fill_node_buckets(pbqp);
1043
1044 #if KAPS_STATISTIC
1045         FILE *fh = fopen("solutions.pb", "a");
1046         fprintf(fh, "Solution");
1047         fclose(fh);
1048 #endif
1049
1050         apply_brute_force_reductions(pbqp);
1051
1052         pbqp->solution = determine_solution(pbqp);
1053
1054 #if KAPS_STATISTIC
1055         fh = fopen("solutions.pb", "a");
1056         fprintf(fh, ": %lld RE:%u R0:%u R1:%u R2:%u RN/BF:%u\n", pbqp->solution,
1057                         pbqp->num_edges, pbqp->num_r0, pbqp->num_r1, pbqp->num_r2,
1058                         pbqp->num_bf);
1059         fclose(fh);
1060 #endif
1061
1062         /* Solve reduced nodes. */
1063         back_propagate(pbqp);
1064
1065         free_buckets();
1066 }
1067
1068 void back_propagate_RI(pbqp *pbqp, pbqp_node *node)
1069 {
1070         pbqp_edge   *edge;
1071         pbqp_node   *other;
1072         pbqp_matrix *mat;
1073         vector      *vec;
1074         int          is_src;
1075
1076         assert(pbqp);
1077         assert(node);
1078
1079         edge = node->edges[0];
1080         mat = edge->costs;
1081         is_src = edge->src == node;
1082         vec = node->costs;
1083
1084         if (is_src) {
1085                 other = edge->tgt;
1086                 assert(other);
1087
1088                 /* Update pointer for brute force solver. */
1089                 other = pbqp->nodes[other->index];
1090
1091                 node->solution = pbqp_matrix_get_col_min_index(mat, other->solution, vec);
1092         } else {
1093                 other = edge->src;
1094                 assert(other);
1095
1096                 /* Update pointer for brute force solver. */
1097                 other = pbqp->nodes[other->index];
1098
1099                 node->solution = pbqp_matrix_get_row_min_index(mat, other->solution, vec);
1100         }
1101
1102         if (pbqp->dump_file) {
1103                 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
1104         }
1105 }
1106
1107 void back_propagate_RII(pbqp *pbqp, pbqp_node *node)
1108 {
1109         pbqp_edge   *src_edge   = node->edges[0];
1110         pbqp_edge   *tgt_edge   = node->edges[1];
1111         int          src_is_src = src_edge->src == node;
1112         int          tgt_is_src = tgt_edge->src == node;
1113         pbqp_matrix *src_mat;
1114         pbqp_matrix *tgt_mat;
1115         pbqp_node   *src_node;
1116         pbqp_node   *tgt_node;
1117         vector      *vec;
1118         vector      *node_vec;
1119         unsigned     col_index;
1120         unsigned     row_index;
1121
1122         assert(pbqp);
1123
1124         if (src_is_src) {
1125                 src_node = src_edge->tgt;
1126         } else {
1127                 src_node = src_edge->src;
1128         }
1129
1130         if (tgt_is_src) {
1131                 tgt_node = tgt_edge->tgt;
1132         } else {
1133                 tgt_node = tgt_edge->src;
1134         }
1135
1136         /* Swap nodes if necessary. */
1137         if (tgt_node->index < src_node->index) {
1138                 pbqp_node *tmp_node;
1139                 pbqp_edge *tmp_edge;
1140
1141                 tmp_node = src_node;
1142                 src_node = tgt_node;
1143                 tgt_node = tmp_node;
1144
1145                 tmp_edge = src_edge;
1146                 src_edge = tgt_edge;
1147                 tgt_edge = tmp_edge;
1148
1149                 src_is_src = src_edge->src == node;
1150                 tgt_is_src = tgt_edge->src == node;
1151         }
1152
1153         /* Update pointer for brute force solver. */
1154         src_node = pbqp->nodes[src_node->index];
1155         tgt_node = pbqp->nodes[tgt_node->index];
1156
1157         src_mat = src_edge->costs;
1158         tgt_mat = tgt_edge->costs;
1159
1160         node_vec = node->costs;
1161
1162         row_index = src_node->solution;
1163         col_index = tgt_node->solution;
1164
1165         vec = vector_copy(pbqp, node_vec);
1166
1167         if (src_is_src) {
1168                 vector_add_matrix_col(vec, src_mat, row_index);
1169         } else {
1170                 vector_add_matrix_row(vec, src_mat, row_index);
1171         }
1172
1173         if (tgt_is_src) {
1174                 vector_add_matrix_col(vec, tgt_mat, col_index);
1175         } else {
1176                 vector_add_matrix_row(vec, tgt_mat, col_index);
1177         }
1178
1179         node->solution = vector_get_min_index(vec);
1180         if (pbqp->dump_file) {
1181                 fprintf(pbqp->dump_file, "node n%d is set to %d<br>\n", node->index, node->solution);
1182         }
1183
1184         obstack_free(&pbqp->obstack, vec);
1185 }
1186
1187 int node_is_reduced(pbqp_node *node)
1188 {
1189         if (!reduced_bucket) return 0;
1190
1191         if (pbqp_node_get_degree(node) == 0) return 1;
1192
1193         return node_bucket_contains(reduced_bucket, node);
1194 }