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