2 * Copyright (C) 2005-2011 University of Karlsruhe. All right reserved.
4 * This file is part of libFirm.
6 * This file may be distributed and/or modified under the terms of the
7 * GNU General Public License version 2 as published by the Free Software
8 * Foundation and appearing in the file LICENSE.GPL included in the
9 * packaging of this file.
11 * Licensees holding valid libFirm Professional Edition licenses may use
12 * this file in accordance with the libFirm Commercial License.
13 * Agreement provided with the Software.
15 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
16 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * @brief Sparse matrix storage with linked lists for rows and cols.
23 * @author Daniel Grund
25 * Sparse matrix storage with linked lists for rows and cols.
26 * Matrix is optimized for left-to-right and top-to-bottom access.
27 * Complexity is O(1) then.
28 * Random access or right-to-left and bottom-to-top is O(m*n).
37 #include "sp_matrix.h"
43 typedef enum iter_direction_t {
48 * Embedded list pointer.
50 typedef struct sp_matrix_list_head_t {
51 struct sp_matrix_list_head_t *next;
52 } sp_matrix_list_head_t;
57 typedef struct entry_t {
58 sp_matrix_list_head_t col_chain; /**< points to next element in same column */
59 sp_matrix_list_head_t row_chain; /**< points to next element in same row */
60 matrix_elem_t e; /**< The actual element */
64 /* These specify the dimensions of the matrix.
65 * They equal the largest values ever used in matrix_set */
67 /* These are the dimensions of allocated arrays below.
68 * rowc >= maxrow and colc >= maxcol hold. */
70 /* number of entries */
72 /* arrays of sp_matrix_list_head* as entry-points to rows and cols */
73 sp_matrix_list_head_t **rows, **cols;
74 /* for iteration: first is to remember start-point;
75 * last was returned just before
76 * next is used in case last was removed from list */
78 sp_matrix_list_head_t *first, *last, *next;
79 int iter_row; /* used for iteration over all elements */
80 /* for each column the last inserted element in col list */
81 sp_matrix_list_head_t **last_col_el;
82 /* for each row the last inserted element in row list */
83 sp_matrix_list_head_t **last_row_el;
86 #define SP_MATRIX_INIT_LIST_HEAD(ptr) do { (ptr)->next = NULL; } while (0)
88 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
89 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
91 #define is_empty_row(row) (row > m->maxrow || (m->rows[row])->next == NULL)
92 #define is_empty_col(col) (col > m->maxcol || (m->cols[col])->next == NULL)
94 #define list_entry_by_col(h) (&_container_of(h, entry_t, col_chain)->e)
95 #define list_entry_by_row(h) (&_container_of(h, entry_t, row_chain)->e)
98 * Returns the size of a single matrix element.
100 unsigned matrix_get_elem_size(void)
102 return sizeof(entry_t);
106 * Returns the new size for an array of size old_size,
107 * which must at least store an entry at position min.
109 static inline int m_new_size(int old_size, int min)
112 assert(min >= old_size);
117 assert(bits < sizeof(min) * 8 - 1);
122 * Allocates space for @p count entries in the rows array and
123 * initializes all entries from @p start to the end.
125 static inline void m_alloc_row(sp_matrix_t *m, int start, int count)
130 m->rows = XREALLOC(m->rows, sp_matrix_list_head_t *, m->rowc);
131 m->last_row_el = XREALLOC(m->last_row_el, sp_matrix_list_head_t *, m->rowc);
133 for (p = start; p < m->rowc; ++p) {
134 m->rows[p] = XMALLOC(sp_matrix_list_head_t);
135 SP_MATRIX_INIT_LIST_HEAD(m->rows[p]);
136 m->last_row_el[p] = m->rows[p];
141 * Allocates space for @p count entries in the cols array and
142 * initializes all entries from @p start to the end.
144 static inline void m_alloc_col(sp_matrix_t *m, int start, int count)
149 m->cols = XREALLOC(m->cols, sp_matrix_list_head_t*, m->colc);
150 m->last_col_el = XREALLOC(m->last_col_el, sp_matrix_list_head_t*, m->colc);
152 for (p = start; p < m->colc; ++p) {
153 m->cols[p] = XMALLOC(sp_matrix_list_head_t);
154 SP_MATRIX_INIT_LIST_HEAD(m->cols[p]);
155 m->last_col_el[p] = m->cols[p];
160 * Searches in row @p row for the matrix element m[row, col], starting at element @p start.
161 * @return If the element exists:
162 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
163 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
164 * @p prev_prev always points to the previous element of @p prev
166 static inline matrix_elem_t *m_search_in_row_from(const sp_matrix_t *m,
167 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
169 sp_matrix_list_head_t *row_start;
170 matrix_elem_t *res = NULL;
172 row_start = m->rows[row];
175 while ((*prev)->next != NULL && list_entry_by_row((*prev)->next)->col <= col) {
176 (*prev_prev) = (*prev);
177 *prev = (*prev)->next;
180 if (*prev != row_start) {
181 matrix_elem_t *me = list_entry_by_row(*prev);
183 if (me->row == row && me->col == col)
188 m->last_row_el[row] = *prev;
195 * Searches in row @p row for the matrix element m[row, col].
196 * @return If the element exists:
197 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
198 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
199 * @p prev_prev always points to the previous element of @p prev
201 static inline matrix_elem_t *m_search_in_row(const sp_matrix_t *m,
202 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
204 sp_matrix_list_head_t *start = m->rows[row];
208 if (m->last_row_el[row] != start) {
209 matrix_elem_t *el = list_entry_by_row(m->last_row_el[row]);
211 *prev_prev = start = m->last_row_el[row];
215 return m_search_in_row_from(m, row, col, start, prev, prev_prev);
219 * Searches in col @p col for the matrix element m[row, col], starting at @p start.
220 * @return If the element exists:
221 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
222 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
223 * @p prev_prev always points to the previous element of @p prev
225 static inline matrix_elem_t *m_search_in_col_from(const sp_matrix_t *m,
226 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
228 sp_matrix_list_head_t *col_start;
229 matrix_elem_t *res = NULL;
231 col_start = m->cols[col];
234 while ((*prev)->next != NULL && list_entry_by_col((*prev)->next)->row <= row) {
235 *prev_prev = (*prev);
236 *prev = (*prev)->next;
239 if (*prev != col_start) {
240 matrix_elem_t *me = list_entry_by_col(*prev);
242 if (me->row == row && me->col == col)
247 m->last_col_el[col] = *prev;
254 * Searches in col @p col for the matrix element m[row, col].
255 * @return If the element exists:
256 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
257 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
258 * @p prev_prev always points to the previous element of @p prev
260 static inline matrix_elem_t *m_search_in_col(const sp_matrix_t *m,
261 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
263 sp_matrix_list_head_t *start = m->cols[col];
267 if (m->last_col_el[col] != start) {
268 matrix_elem_t *el = list_entry_by_col(m->last_col_el[col]);
270 *prev_prev = start = m->last_col_el[col];
274 return m_search_in_col_from(m, row, col, start, prev, prev_prev);
277 sp_matrix_t *new_matrix(int row_init, int col_init)
279 sp_matrix_t *res = XMALLOCZ(sp_matrix_t);
282 m_alloc_row(res, 0, MAX(0, row_init));
283 m_alloc_col(res, 0, MAX(0, col_init));
287 void del_matrix(sp_matrix_t *m)
291 for (i = 0; i < m->rowc; ++i) {
292 if (! is_empty_row(i)) {
294 sp_matrix_list_head_t *n;
296 n = m->rows[i]->next;
298 /* get current matrix element */
299 e = _container_of(n, entry_t, row_chain);
307 for (i = 0; i < m->colc; ++i)
309 xfree(m->last_col_el);
310 xfree(m->last_row_el);
316 void matrix_set(sp_matrix_t *m, int row, int col, double val)
318 matrix_elem_t *me = NULL;
320 sp_matrix_list_head_t *leftof, *above;
321 sp_matrix_list_head_t *prev_leftof, *prev_above;
323 /* if necessary enlarge the matrix */
324 if (row > m->maxrow) {
327 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
329 if (col > m->maxcol) {
332 m_alloc_col(m, m->colc, m_new_size(m->colc, col));
335 /* search for existing entry */
336 if (m->maxrow < m->maxcol)
337 me = m_search_in_col(m, row, col, &above, &prev_above);
339 me = m_search_in_row(m, row, col, &leftof, &prev_leftof);
341 /* if it exists, set the value and return */
344 me->val = (float)val;
346 entr = _container_of(me, entry_t, e);
348 /* remove row_chain entry */
350 prev_leftof->next = entr->row_chain.next;
352 m->rows[row]->next = entr->row_chain.next;
354 /* remove col_chain entry */
356 prev_above->next = entr->col_chain.next;
358 m->cols[col]->next = entr->col_chain.next;
360 entr->row_chain.next = NULL;
361 entr->col_chain.next = NULL;
363 /* set the last pointer to the "previous" element */
364 if (m->last_col_el[col] == &entr->col_chain ||
365 m->last_row_el[row] == &entr->row_chain)
367 m->last_col_el[col] = prev_above ? prev_above : m->cols[col];
368 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
377 /* if it does not exist and 0 should be set just quit */
381 /* if it does not exist and val != 0 search the other direction */
382 if (m->maxrow >= m->maxcol)
383 m_search_in_col(m, row, col, &above, &prev_above);
385 m_search_in_row(m, row, col, &leftof, &prev_leftof);
386 /* now leftof and above are the entry_t's prior the new one in each direction */
388 /* insert new entry */
389 entr = XMALLOC(entry_t);
392 entr->e.val = (float)val;
394 /* add row_chain entry */
395 entr->row_chain.next = leftof->next;
396 leftof->next = &entr->row_chain;
398 /* add col_chain entry */
399 entr->col_chain.next = above->next;
400 above->next = &entr->col_chain;
402 m->last_col_el[col] = &entr->col_chain;
403 m->last_row_el[row] = &entr->row_chain;
408 void matrix_set_row_bulk(sp_matrix_t *m, int row, int *cols, int num_cols, double val)
410 matrix_elem_t *me = NULL;
413 sp_matrix_list_head_t *leftof, *above;
414 sp_matrix_list_head_t *prev_leftof, *prev_above;
416 /* if necessary enlarge the matrix */
417 if (row > m->maxrow) {
420 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
422 if (cols[num_cols - 1] > m->maxcol) {
423 m->maxcol = cols[num_cols - 1];
424 if (cols[num_cols - 1] >= m->colc)
425 m_alloc_col(m, m->colc, m_new_size(m->colc, cols[num_cols - 1]));
428 /* set start values */
432 for (i = 0; i < num_cols; ++i) {
433 /* search for existing entry */
434 me = m_search_in_row(m, row, cols[i], &leftof, &prev_leftof);
436 /* if it exists, set the value and return */
439 me->val = (float)val;
441 entr = _container_of(me, entry_t, e);
443 /* remove row_chain entry */
445 prev_leftof->next = entr->row_chain.next;
447 m->rows[row]->next = entr->row_chain.next;
449 /* remove col_chain entry */
451 prev_above->next = entr->col_chain.next;
453 m->cols[cols[i]]->next = entr->col_chain.next;
455 entr->row_chain.next = NULL;
456 entr->col_chain.next = NULL;
458 /* set the last pointer to the "previous" element */
459 if (m->last_col_el[cols[i]] == &entr->col_chain ||
460 m->last_row_el[row] == &entr->row_chain)
462 m->last_col_el[cols[i]] = prev_above ? prev_above : m->cols[cols[i]];
463 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
473 /* if it does not exist and 0 should be set just quit */
477 /* we have to search the col list as well, to get the above pointer */
478 m_search_in_col(m, row, cols[i], &above, &prev_above);
480 /* now leftof and above are the entry_t's prior the new one in each direction */
482 /* insert new entry */
483 entr = XMALLOC(entry_t);
485 entr->e.col = cols[i];
486 entr->e.val = (float)val;
488 m->last_col_el[cols[i]] = &entr->col_chain;
489 m->last_row_el[row] = &entr->row_chain;
491 /* add row_chain entry */
492 entr->row_chain.next = leftof->next;
493 leftof->next = &entr->row_chain;
495 /* add col_chain entry */
496 entr->col_chain.next = above->next;
497 above->next = &entr->col_chain;
503 double matrix_get(const sp_matrix_t *m, int row, int col)
505 sp_matrix_list_head_t *dummy, *dummy2;
508 if (is_empty_row(row) || is_empty_col(col))
511 if (m->maxrow < m->maxcol)
512 me = m_search_in_col(m, row, col, &dummy, &dummy2);
514 me = m_search_in_row(m, row, col, &dummy, &dummy2);
517 assert(me->col == col && me->row == row);
519 return me ? me->val : 0.0;
522 int matrix_get_entries(const sp_matrix_t *m)
527 int matrix_get_rowcount(const sp_matrix_t *m)
529 return m->maxrow + 1;
532 int matrix_get_colcount(const sp_matrix_t *m)
534 return m->maxcol + 1;
537 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row)
539 if (is_empty_row(row))
543 m->first = m->rows[row];
544 m->last = m->first->next;
545 m->next = m->last ? m->last->next : NULL;
547 assert (list_entry_by_row(m->last)->row == row);
549 return list_entry_by_row(m->last);
552 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col)
554 if (is_empty_col(col))
558 m->first = m->cols[col];
559 m->last = m->first->next;
560 m->next = m->last ? m->last->next : NULL;
562 assert (list_entry_by_col(m->last)->col == col);
564 return list_entry_by_col(m->last);
567 static inline const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow)
569 const matrix_elem_t *res;
572 for (i = startrow; i <= m->maxrow; ++i) {
573 res = matrix_row_first(m, i);
584 const matrix_elem_t *matrix_first(sp_matrix_t *m)
586 return matrix_first_from(m, 0);
589 const matrix_elem_t *matrix_next(sp_matrix_t *m)
591 assert(m->first && "Start iteration with matrix_???_first, before calling me!");
593 if (m->next == NULL) {
595 return matrix_first_from(m, ++m->iter_row);
601 m->next = m->next->next;
604 return list_entry_by_col(m->last);
605 else /* right or all */
606 return list_entry_by_row(m->last);
609 static int cmp_count(const void *e1, const void *e2)
611 return (int *)e2 - (int *)e1;
614 static inline void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow)
616 bitset_set(fullrow, row);
617 matrix_foreach_in_col(m, row, e) {
618 if (! bitset_is_set(fullrow, e->row)) {
619 assert(0.0 == matrix_get(m, e->col, e->row));
620 matrix_set(m, e->col, e->row, e->val);
621 matrix_set(m, e->row, e->col, 0.0);
626 void matrix_optimize(sp_matrix_t *m)
632 size = MAX(m->maxcol, m->maxrow)+1;
634 /* kill all double-entries (Mij and Mji are set) */
635 matrix_foreach(m, e) {
638 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will always have the same color as root");
639 t_val = matrix_get(m, e->col, e->row);
640 if (fabs(t_val) > 1e-10) {
641 matrix_set(m, e->col, e->row, 0);
642 matrix_set(m, e->row, e->col, e->val + t_val);
646 c = ALLOCAN(int, size);
648 fullrow = bitset_alloca(size);
650 /* kill 'all' rows containing only 1 entry */
653 /* count elements in rows */
654 memset(c, 0, size * sizeof(*c));
659 for (i = 0; i<size; ++i)
660 if (c[i] == 1 && ! bitset_is_set(fullrow, i)) {
662 /* if the other row isn't empty move the e in there, else fill e's row */
663 matrix_elem_t const *const e = matrix_row_first(m, i);
666 matrix_fill_row(m, e->col, fullrow);
668 matrix_fill_row(m, e->row, fullrow);
674 memset(c, 0, size * sizeof(*c));
678 qsort(c, size, sizeof(*c), cmp_count);
680 for (i = 0; i < size; ++i) {
681 if (! bitset_is_set(fullrow, i))
682 matrix_fill_row(m, i, fullrow);
686 void matrix_dump(sp_matrix_t *m, FILE *out, int factor)
690 for (i = 0; i <= m->maxrow; ++i) {
692 matrix_foreach_in_row(m, i, e) {
693 for (o = last_idx + 1; o < e->col; ++o)
694 fprintf(out, " %4.1f" , 0.0);
696 fprintf(out, " %4.1f", factor * e->val);
700 for (o = last_idx + 1; o <= m->maxcol; ++o)
701 fprintf(out, " %4.1f" , 0.0);
707 void matrix_self_test(int d)
710 sp_matrix_t *m = new_matrix(10, 10);
712 for (i = 0; i < d; ++i)
713 for (o = 0; o < d; ++o)
714 matrix_set(m, i, o, i*o);
716 for (i = 0; i < d; ++i)
717 for (o = 0; o<d; ++o)
718 assert(matrix_get(m, i, o) == i*o);
721 matrix_foreach_in_row(m,1,e) {
725 assert(!matrix_next(m)); /*iter must finish */
728 matrix_foreach_in_col(m,d-1,e) {
732 assert(!matrix_next(m));
734 m = new_matrix(16,16);
735 matrix_set(m, 1,1,9);
736 matrix_set(m, 1,2,8);
737 matrix_set(m, 1,3,7);
739 matrix_set(m, 1,3,6);
740 matrix_set(m, 1,2,5);
741 matrix_set(m, 1,1,4);
743 matrix_foreach_in_row(m, 1, e) {
744 assert(e->row == 1 && e->col == i && e->val == i+3);
751 matrix_set(m, 1,1,1);
752 matrix_set(m, 2,2,2);
753 matrix_set(m, 3,3,3);
754 matrix_set(m, 3,5,4);
755 matrix_set(m, 4,4,5);
756 matrix_set(m, 5,5,6);
759 assert(e->val == ++i);
761 matrix_set(m, 1,1,0);
762 assert(5 == matrix_get_entries(m));