2 * This file is part of libFirm.
3 * Copyright (C) 2012 University of Karlsruhe.
8 * @brief Sparse matrix storage with linked lists for rows and cols.
11 * Sparse matrix storage with linked lists for rows and cols.
12 * Matrix is optimized for left-to-right and top-to-bottom access.
13 * Complexity is O(1) then.
14 * Random access or right-to-left and bottom-to-top is O(m*n).
23 #include "sp_matrix.h"
29 typedef enum iter_direction_t {
34 * Embedded list pointer.
36 typedef struct sp_matrix_list_head_t {
37 struct sp_matrix_list_head_t *next;
38 } sp_matrix_list_head_t;
43 typedef struct entry_t {
44 sp_matrix_list_head_t col_chain; /**< points to next element in same column */
45 sp_matrix_list_head_t row_chain; /**< points to next element in same row */
46 matrix_elem_t e; /**< The actual element */
50 /* These specify the dimensions of the matrix.
51 * They equal the largest values ever used in matrix_set */
53 /* These are the dimensions of allocated arrays below.
54 * rowc >= maxrow and colc >= maxcol hold. */
56 /* number of entries */
58 /* arrays of sp_matrix_list_head* as entry-points to rows and cols */
59 sp_matrix_list_head_t **rows, **cols;
60 /* for iteration: first is to remember start-point;
61 * last was returned just before
62 * next is used in case last was removed from list */
64 sp_matrix_list_head_t *first, *last, *next;
65 int iter_row; /* used for iteration over all elements */
66 /* for each column the last inserted element in col list */
67 sp_matrix_list_head_t **last_col_el;
68 /* for each row the last inserted element in row list */
69 sp_matrix_list_head_t **last_row_el;
72 #define SP_MATRIX_INIT_LIST_HEAD(ptr) do { (ptr)->next = NULL; } while (0)
74 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
75 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
77 #define is_empty_row(row) (row > m->maxrow || (m->rows[row])->next == NULL)
78 #define is_empty_col(col) (col > m->maxcol || (m->cols[col])->next == NULL)
80 #define list_entry_by_col(h) (&_container_of(h, entry_t, col_chain)->e)
81 #define list_entry_by_row(h) (&_container_of(h, entry_t, row_chain)->e)
84 * Returns the size of a single matrix element.
86 unsigned matrix_get_elem_size(void)
88 return sizeof(entry_t);
92 * Returns the new size for an array of size old_size,
93 * which must at least store an entry at position min.
95 static inline int m_new_size(int old_size, int min)
98 assert(min >= old_size);
103 assert(bits < sizeof(min) * 8 - 1);
108 * Allocates space for @p count entries in the rows array and
109 * initializes all entries from @p start to the end.
111 static inline void m_alloc_row(sp_matrix_t *m, int start, int count)
116 m->rows = XREALLOC(m->rows, sp_matrix_list_head_t *, m->rowc);
117 m->last_row_el = XREALLOC(m->last_row_el, sp_matrix_list_head_t *, m->rowc);
119 for (p = start; p < m->rowc; ++p) {
120 m->rows[p] = XMALLOC(sp_matrix_list_head_t);
121 SP_MATRIX_INIT_LIST_HEAD(m->rows[p]);
122 m->last_row_el[p] = m->rows[p];
127 * Allocates space for @p count entries in the cols array and
128 * initializes all entries from @p start to the end.
130 static inline void m_alloc_col(sp_matrix_t *m, int start, int count)
135 m->cols = XREALLOC(m->cols, sp_matrix_list_head_t*, m->colc);
136 m->last_col_el = XREALLOC(m->last_col_el, sp_matrix_list_head_t*, m->colc);
138 for (p = start; p < m->colc; ++p) {
139 m->cols[p] = XMALLOC(sp_matrix_list_head_t);
140 SP_MATRIX_INIT_LIST_HEAD(m->cols[p]);
141 m->last_col_el[p] = m->cols[p];
146 * Searches in row @p row for the matrix element m[row, col], starting at element @p start.
147 * @return If the element exists:
148 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
149 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
150 * @p prev_prev always points to the previous element of @p prev
152 static inline matrix_elem_t *m_search_in_row_from(const sp_matrix_t *m,
153 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
155 sp_matrix_list_head_t *row_start;
156 matrix_elem_t *res = NULL;
158 row_start = m->rows[row];
161 while ((*prev)->next != NULL && list_entry_by_row((*prev)->next)->col <= col) {
162 (*prev_prev) = (*prev);
163 *prev = (*prev)->next;
166 if (*prev != row_start) {
167 matrix_elem_t *me = list_entry_by_row(*prev);
169 if (me->row == row && me->col == col)
174 m->last_row_el[row] = *prev;
181 * Searches in row @p row for the matrix element m[row, col].
182 * @return If the element exists:
183 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
184 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
185 * @p prev_prev always points to the previous element of @p prev
187 static inline matrix_elem_t *m_search_in_row(const sp_matrix_t *m,
188 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
190 sp_matrix_list_head_t *start = m->rows[row];
194 if (m->last_row_el[row] != start) {
195 matrix_elem_t *el = list_entry_by_row(m->last_row_el[row]);
197 *prev_prev = start = m->last_row_el[row];
201 return m_search_in_row_from(m, row, col, start, prev, prev_prev);
205 * Searches in col @p col for the matrix element m[row, col], starting at @p start.
206 * @return If the element exists:
207 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
208 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
209 * @p prev_prev always points to the previous element of @p prev
211 static inline matrix_elem_t *m_search_in_col_from(const sp_matrix_t *m,
212 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
214 sp_matrix_list_head_t *col_start;
215 matrix_elem_t *res = NULL;
217 col_start = m->cols[col];
220 while ((*prev)->next != NULL && list_entry_by_col((*prev)->next)->row <= row) {
221 *prev_prev = (*prev);
222 *prev = (*prev)->next;
225 if (*prev != col_start) {
226 matrix_elem_t *me = list_entry_by_col(*prev);
228 if (me->row == row && me->col == col)
233 m->last_col_el[col] = *prev;
240 * Searches in col @p col for the matrix element m[row, col].
241 * @return If the element exists:
242 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
243 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
244 * @p prev_prev always points to the previous element of @p prev
246 static inline matrix_elem_t *m_search_in_col(const sp_matrix_t *m,
247 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
249 sp_matrix_list_head_t *start = m->cols[col];
253 if (m->last_col_el[col] != start) {
254 matrix_elem_t *el = list_entry_by_col(m->last_col_el[col]);
256 *prev_prev = start = m->last_col_el[col];
260 return m_search_in_col_from(m, row, col, start, prev, prev_prev);
263 sp_matrix_t *new_matrix(int row_init, int col_init)
265 sp_matrix_t *res = XMALLOCZ(sp_matrix_t);
268 m_alloc_row(res, 0, MAX(0, row_init));
269 m_alloc_col(res, 0, MAX(0, col_init));
273 void del_matrix(sp_matrix_t *m)
277 for (i = 0; i < m->rowc; ++i) {
278 if (! is_empty_row(i)) {
280 sp_matrix_list_head_t *n;
282 n = m->rows[i]->next;
284 /* get current matrix element */
285 e = _container_of(n, entry_t, row_chain);
293 for (i = 0; i < m->colc; ++i)
295 xfree(m->last_col_el);
296 xfree(m->last_row_el);
302 void matrix_set(sp_matrix_t *m, int row, int col, double val)
304 matrix_elem_t *me = NULL;
306 sp_matrix_list_head_t *leftof, *above;
307 sp_matrix_list_head_t *prev_leftof, *prev_above;
309 /* if necessary enlarge the matrix */
310 if (row > m->maxrow) {
313 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
315 if (col > m->maxcol) {
318 m_alloc_col(m, m->colc, m_new_size(m->colc, col));
321 /* search for existing entry */
322 if (m->maxrow < m->maxcol)
323 me = m_search_in_col(m, row, col, &above, &prev_above);
325 me = m_search_in_row(m, row, col, &leftof, &prev_leftof);
327 /* if it exists, set the value and return */
330 me->val = (float)val;
332 entr = _container_of(me, entry_t, e);
334 /* remove row_chain entry */
336 prev_leftof->next = entr->row_chain.next;
338 m->rows[row]->next = entr->row_chain.next;
340 /* remove col_chain entry */
342 prev_above->next = entr->col_chain.next;
344 m->cols[col]->next = entr->col_chain.next;
346 entr->row_chain.next = NULL;
347 entr->col_chain.next = NULL;
349 /* set the last pointer to the "previous" element */
350 if (m->last_col_el[col] == &entr->col_chain ||
351 m->last_row_el[row] == &entr->row_chain)
353 m->last_col_el[col] = prev_above ? prev_above : m->cols[col];
354 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
363 /* if it does not exist and 0 should be set just quit */
367 /* if it does not exist and val != 0 search the other direction */
368 if (m->maxrow >= m->maxcol)
369 m_search_in_col(m, row, col, &above, &prev_above);
371 m_search_in_row(m, row, col, &leftof, &prev_leftof);
372 /* now leftof and above are the entry_t's prior the new one in each direction */
374 /* insert new entry */
375 entr = XMALLOC(entry_t);
378 entr->e.val = (float)val;
380 /* add row_chain entry */
381 entr->row_chain.next = leftof->next;
382 leftof->next = &entr->row_chain;
384 /* add col_chain entry */
385 entr->col_chain.next = above->next;
386 above->next = &entr->col_chain;
388 m->last_col_el[col] = &entr->col_chain;
389 m->last_row_el[row] = &entr->row_chain;
394 void matrix_set_row_bulk(sp_matrix_t *m, int row, int *cols, int num_cols, double val)
396 matrix_elem_t *me = NULL;
399 sp_matrix_list_head_t *leftof, *above;
400 sp_matrix_list_head_t *prev_leftof, *prev_above;
402 /* if necessary enlarge the matrix */
403 if (row > m->maxrow) {
406 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
408 if (cols[num_cols - 1] > m->maxcol) {
409 m->maxcol = cols[num_cols - 1];
410 if (cols[num_cols - 1] >= m->colc)
411 m_alloc_col(m, m->colc, m_new_size(m->colc, cols[num_cols - 1]));
414 /* set start values */
418 for (i = 0; i < num_cols; ++i) {
419 /* search for existing entry */
420 me = m_search_in_row(m, row, cols[i], &leftof, &prev_leftof);
422 /* if it exists, set the value and return */
425 me->val = (float)val;
427 entr = _container_of(me, entry_t, e);
429 /* remove row_chain entry */
431 prev_leftof->next = entr->row_chain.next;
433 m->rows[row]->next = entr->row_chain.next;
435 /* remove col_chain entry */
437 prev_above->next = entr->col_chain.next;
439 m->cols[cols[i]]->next = entr->col_chain.next;
441 entr->row_chain.next = NULL;
442 entr->col_chain.next = NULL;
444 /* set the last pointer to the "previous" element */
445 if (m->last_col_el[cols[i]] == &entr->col_chain ||
446 m->last_row_el[row] == &entr->row_chain)
448 m->last_col_el[cols[i]] = prev_above ? prev_above : m->cols[cols[i]];
449 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
459 /* if it does not exist and 0 should be set just quit */
463 /* we have to search the col list as well, to get the above pointer */
464 m_search_in_col(m, row, cols[i], &above, &prev_above);
466 /* now leftof and above are the entry_t's prior the new one in each direction */
468 /* insert new entry */
469 entr = XMALLOC(entry_t);
471 entr->e.col = cols[i];
472 entr->e.val = (float)val;
474 m->last_col_el[cols[i]] = &entr->col_chain;
475 m->last_row_el[row] = &entr->row_chain;
477 /* add row_chain entry */
478 entr->row_chain.next = leftof->next;
479 leftof->next = &entr->row_chain;
481 /* add col_chain entry */
482 entr->col_chain.next = above->next;
483 above->next = &entr->col_chain;
489 double matrix_get(const sp_matrix_t *m, int row, int col)
491 sp_matrix_list_head_t *dummy, *dummy2;
494 if (is_empty_row(row) || is_empty_col(col))
497 if (m->maxrow < m->maxcol)
498 me = m_search_in_col(m, row, col, &dummy, &dummy2);
500 me = m_search_in_row(m, row, col, &dummy, &dummy2);
503 assert(me->col == col && me->row == row);
505 return me ? me->val : 0.0;
508 int matrix_get_entries(const sp_matrix_t *m)
513 int matrix_get_rowcount(const sp_matrix_t *m)
515 return m->maxrow + 1;
518 int matrix_get_colcount(const sp_matrix_t *m)
520 return m->maxcol + 1;
523 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row)
525 if (is_empty_row(row))
529 m->first = m->rows[row];
530 m->last = m->first->next;
531 m->next = m->last ? m->last->next : NULL;
533 assert (list_entry_by_row(m->last)->row == row);
535 return list_entry_by_row(m->last);
538 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col)
540 if (is_empty_col(col))
544 m->first = m->cols[col];
545 m->last = m->first->next;
546 m->next = m->last ? m->last->next : NULL;
548 assert (list_entry_by_col(m->last)->col == col);
550 return list_entry_by_col(m->last);
553 static inline const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow)
555 const matrix_elem_t *res;
558 for (i = startrow; i <= m->maxrow; ++i) {
559 res = matrix_row_first(m, i);
570 const matrix_elem_t *matrix_first(sp_matrix_t *m)
572 return matrix_first_from(m, 0);
575 const matrix_elem_t *matrix_next(sp_matrix_t *m)
577 assert(m->first && "Start iteration with matrix_???_first, before calling me!");
579 if (m->next == NULL) {
581 return matrix_first_from(m, ++m->iter_row);
587 m->next = m->next->next;
590 return list_entry_by_col(m->last);
591 else /* right or all */
592 return list_entry_by_row(m->last);
595 static int cmp_count(const void *e1, const void *e2)
597 return (int *)e2 - (int *)e1;
600 static inline void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow)
602 bitset_set(fullrow, row);
603 matrix_foreach_in_col(m, row, e) {
604 if (! bitset_is_set(fullrow, e->row)) {
605 assert(0.0 == matrix_get(m, e->col, e->row));
606 matrix_set(m, e->col, e->row, e->val);
607 matrix_set(m, e->row, e->col, 0.0);
612 void matrix_optimize(sp_matrix_t *m)
618 size = MAX(m->maxcol, m->maxrow)+1;
620 /* kill all double-entries (Mij and Mji are set) */
621 matrix_foreach(m, e) {
624 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will always have the same color as root");
625 t_val = matrix_get(m, e->col, e->row);
626 if (fabs(t_val) > 1e-10) {
627 matrix_set(m, e->col, e->row, 0);
628 matrix_set(m, e->row, e->col, e->val + t_val);
632 c = ALLOCAN(int, size);
634 fullrow = bitset_alloca(size);
636 /* kill 'all' rows containing only 1 entry */
639 /* count elements in rows */
640 memset(c, 0, size * sizeof(*c));
645 for (i = 0; i<size; ++i)
646 if (c[i] == 1 && ! bitset_is_set(fullrow, i)) {
648 /* if the other row isn't empty move the e in there, else fill e's row */
649 matrix_elem_t const *const e = matrix_row_first(m, i);
652 matrix_fill_row(m, e->col, fullrow);
654 matrix_fill_row(m, e->row, fullrow);
660 memset(c, 0, size * sizeof(*c));
664 qsort(c, size, sizeof(*c), cmp_count);
666 for (i = 0; i < size; ++i) {
667 if (! bitset_is_set(fullrow, i))
668 matrix_fill_row(m, i, fullrow);
672 void matrix_dump(sp_matrix_t *m, FILE *out, int factor)
676 for (i = 0; i <= m->maxrow; ++i) {
678 matrix_foreach_in_row(m, i, e) {
679 for (o = last_idx + 1; o < e->col; ++o)
680 fprintf(out, " %4.1f" , 0.0);
682 fprintf(out, " %4.1f", factor * e->val);
686 for (o = last_idx + 1; o <= m->maxcol; ++o)
687 fprintf(out, " %4.1f" , 0.0);
693 void matrix_self_test(int d)
696 sp_matrix_t *m = new_matrix(10, 10);
698 for (i = 0; i < d; ++i)
699 for (o = 0; o < d; ++o)
700 matrix_set(m, i, o, i*o);
702 for (i = 0; i < d; ++i)
703 for (o = 0; o<d; ++o)
704 assert(matrix_get(m, i, o) == i*o);
707 matrix_foreach_in_row(m,1,e) {
711 assert(!matrix_next(m)); /*iter must finish */
714 matrix_foreach_in_col(m,d-1,e) {
718 assert(!matrix_next(m));
720 m = new_matrix(16,16);
721 matrix_set(m, 1,1,9);
722 matrix_set(m, 1,2,8);
723 matrix_set(m, 1,3,7);
725 matrix_set(m, 1,3,6);
726 matrix_set(m, 1,2,5);
727 matrix_set(m, 1,1,4);
729 matrix_foreach_in_row(m, 1, e) {
730 assert(e->row == 1 && e->col == i && e->val == i+3);
737 matrix_set(m, 1,1,1);
738 matrix_set(m, 2,2,2);
739 matrix_set(m, 3,3,3);
740 matrix_set(m, 3,5,4);
741 matrix_set(m, 4,4,5);
742 matrix_set(m, 5,5,6);
745 assert(e->val == ++i);
747 matrix_set(m, 1,1,0);
748 assert(5 == matrix_get_entries(m));