2 * Author: Daniel Grund, Christian Wuerdig
4 * Copyright: (c) Universitaet Karlsruhe
5 * Licence: This file protected by GPL - GNU GENERAL PUBLIC LICENSE.
6 * CVS-Id: $Id: sp_matrix.c 24123 2008-11-28 15:08:27Z mallon $
8 * Sparse matrix storage with linked lists for rows and cols.
9 * Matrix is optimized for left-to-right and top-to-bottom access.
10 * Complexity is O(1) then.
11 * Random access or right-to-left and bottom-to-top is O(m*n).
20 #include "sp_matrix.h"
26 typedef enum iter_direction_t {
31 * Embedded list pointer.
33 typedef struct sp_matrix_list_head_t {
34 struct sp_matrix_list_head_t *next;
35 } sp_matrix_list_head_t;
40 typedef struct entry_t {
41 sp_matrix_list_head_t col_chain; /**< points to next element in same column */
42 sp_matrix_list_head_t row_chain; /**< points to next element in same row */
43 matrix_elem_t e; /**< The actual element */
47 /* These specify the dimensions of the matrix.
48 * They equal the largest values ever used in matrix_set */
50 /* These are the dimensions of allocated arrays below.
51 * rowc >= maxrow and colc >= maxcol hold. */
53 /* number of entries */
55 /* arrays of sp_matrix_list_head* as entry-points to rows and cols */
56 sp_matrix_list_head_t **rows, **cols;
57 /* for iteration: first is to remember start-point;
58 * last was returned just before
59 * next is used in case last was removed from list */
61 sp_matrix_list_head_t *first, *last, *next;
62 int iter_row; /* used for iteration over all elements */
63 /* for each column the last inserted element in col list */
64 sp_matrix_list_head_t **last_col_el;
65 /* for each row the last inserted element in row list */
66 sp_matrix_list_head_t **last_row_el;
69 #define SP_MATRIX_INIT_LIST_HEAD(ptr) do { (ptr)->next = NULL; } while (0)
71 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
72 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
74 #define is_empty_row(row) (row > m->maxrow || (m->rows[row])->next == NULL)
75 #define is_empty_col(col) (col > m->maxcol || (m->cols[col])->next == NULL)
77 #define list_entry_by_col(h) (&_container_of(h, entry_t, col_chain)->e)
78 #define list_entry_by_row(h) (&_container_of(h, entry_t, row_chain)->e)
81 * Returns the size of a single matrix element.
83 unsigned matrix_get_elem_size(void)
85 return sizeof(entry_t);
89 * Returns the new size for an array of size old_size,
90 * which must at least store an entry at position min.
92 static inline int m_new_size(int old_size, int min)
95 assert(min >= old_size);
100 assert(bits < sizeof(min) * 8 - 1);
105 * Allocates space for @p count entries in the rows array and
106 * initializes all entries from @p start to the end.
108 static inline void m_alloc_row(sp_matrix_t *m, int start, int count)
113 m->rows = XREALLOC(m->rows, sp_matrix_list_head_t *, m->rowc);
114 m->last_row_el = XREALLOC(m->last_row_el, sp_matrix_list_head_t *, m->rowc);
116 for (p = start; p < m->rowc; ++p) {
117 m->rows[p] = XMALLOC(sp_matrix_list_head_t);
118 SP_MATRIX_INIT_LIST_HEAD(m->rows[p]);
119 m->last_row_el[p] = m->rows[p];
124 * Allocates space for @p count entries in the cols array and
125 * initializes all entries from @p start to the end.
127 static inline void m_alloc_col(sp_matrix_t *m, int start, int count)
132 m->cols = XREALLOC(m->cols, sp_matrix_list_head_t*, m->colc);
133 m->last_col_el = XREALLOC(m->last_col_el, sp_matrix_list_head_t*, m->colc);
135 for (p = start; p < m->colc; ++p) {
136 m->cols[p] = XMALLOC(sp_matrix_list_head_t);
137 SP_MATRIX_INIT_LIST_HEAD(m->cols[p]);
138 m->last_col_el[p] = m->cols[p];
143 * Searches in row @p row for the matrix element m[row, col], starting at element @p start.
144 * @return If the element exists:
145 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
146 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
147 * @p prev_prev always points to the previous element of @p prev
149 static inline matrix_elem_t *m_search_in_row_from(const sp_matrix_t *m,
150 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
152 sp_matrix_list_head_t *row_start;
153 matrix_elem_t *res = NULL;
155 row_start = m->rows[row];
158 while ((*prev)->next != NULL && list_entry_by_row((*prev)->next)->col <= col) {
159 (*prev_prev) = (*prev);
160 *prev = (*prev)->next;
163 if (*prev != row_start) {
164 matrix_elem_t *me = list_entry_by_row(*prev);
166 if (me->row == row && me->col == col)
171 m->last_row_el[row] = *prev;
178 * Searches in row @p row for the matrix element m[row, col].
179 * @return If the element exists:
180 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
181 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
182 * @p prev_prev always points to the previous element of @p prev
184 static inline matrix_elem_t *m_search_in_row(const sp_matrix_t *m,
185 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
187 sp_matrix_list_head_t *start = m->rows[row];
191 if (m->last_row_el[row] != start) {
192 matrix_elem_t *el = list_entry_by_row(m->last_row_el[row]);
194 *prev_prev = start = m->last_row_el[row];
198 return m_search_in_row_from(m, row, col, start, prev, prev_prev);
202 * Searches in col @p col for the matrix element m[row, col], starting at @p start.
203 * @return If the element exists:
204 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
205 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
206 * @p prev_prev always points to the previous element of @p prev
208 static inline matrix_elem_t *m_search_in_col_from(const sp_matrix_t *m,
209 int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
211 sp_matrix_list_head_t *col_start;
212 matrix_elem_t *res = NULL;
214 col_start = m->cols[col];
217 while ((*prev)->next != NULL && list_entry_by_col((*prev)->next)->row <= row) {
218 *prev_prev = (*prev);
219 *prev = (*prev)->next;
222 if (*prev != col_start) {
223 matrix_elem_t *me = list_entry_by_col(*prev);
225 if (me->row == row && me->col == col)
230 m->last_col_el[col] = *prev;
237 * Searches in col @p col for the matrix element m[row, col].
238 * @return If the element exists:
239 * Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
240 * Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
241 * @p prev_prev always points to the previous element of @p prev
243 static inline matrix_elem_t *m_search_in_col(const sp_matrix_t *m,
244 int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
246 sp_matrix_list_head_t *start = m->cols[col];
250 if (m->last_col_el[col] != start) {
251 matrix_elem_t *el = list_entry_by_col(m->last_col_el[col]);
253 *prev_prev = start = m->last_col_el[col];
257 return m_search_in_col_from(m, row, col, start, prev, prev_prev);
260 sp_matrix_t *new_matrix(int row_init, int col_init)
262 sp_matrix_t *res = XMALLOCZ(sp_matrix_t);
265 m_alloc_row(res, 0, MAX(0, row_init));
266 m_alloc_col(res, 0, MAX(0, col_init));
270 void del_matrix(sp_matrix_t *m)
274 for (i = 0; i < m->rowc; ++i) {
275 if (! is_empty_row(i)) {
277 sp_matrix_list_head_t *n;
279 n = m->rows[i]->next;
281 /* get current matrix element */
282 e = _container_of(n, entry_t, row_chain);
290 for (i = 0; i < m->colc; ++i)
292 xfree(m->last_col_el);
293 xfree(m->last_row_el);
299 void matrix_set(sp_matrix_t *m, int row, int col, double val)
301 matrix_elem_t *me = NULL;
303 sp_matrix_list_head_t *leftof, *above;
304 sp_matrix_list_head_t *prev_leftof, *prev_above;
306 /* if necessary enlarge the matrix */
307 if (row > m->maxrow) {
310 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
312 if (col > m->maxcol) {
315 m_alloc_col(m, m->colc, m_new_size(m->colc, col));
318 /* search for existing entry */
319 if (m->maxrow < m->maxcol)
320 me = m_search_in_col(m, row, col, &above, &prev_above);
322 me = m_search_in_row(m, row, col, &leftof, &prev_leftof);
324 /* if it exists, set the value and return */
327 me->val = (float)val;
329 entr = _container_of(me, entry_t, e);
331 /* remove row_chain entry */
333 prev_leftof->next = entr->row_chain.next;
335 m->rows[row]->next = entr->row_chain.next;
337 /* remove col_chain entry */
339 prev_above->next = entr->col_chain.next;
341 m->cols[col]->next = entr->col_chain.next;
343 entr->row_chain.next = NULL;
344 entr->col_chain.next = NULL;
346 /* set the last pointer to the "previous" element */
347 if (m->last_col_el[col] == &entr->col_chain ||
348 m->last_row_el[row] == &entr->row_chain)
350 m->last_col_el[col] = prev_above ? prev_above : m->cols[col];
351 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
360 /* if it does not exist and 0 should be set just quit */
364 /* if it does not exist and val != 0 search the other direction */
365 if (m->maxrow >= m->maxcol)
366 m_search_in_col(m, row, col, &above, &prev_above);
368 m_search_in_row(m, row, col, &leftof, &prev_leftof);
369 /* now leftof and above are the entry_t's prior the new one in each direction */
371 /* insert new entry */
372 entr = XMALLOC(entry_t);
375 entr->e.val = (float)val;
377 /* add row_chain entry */
378 entr->row_chain.next = leftof->next;
379 leftof->next = &entr->row_chain;
381 /* add col_chain entry */
382 entr->col_chain.next = above->next;
383 above->next = &entr->col_chain;
385 m->last_col_el[col] = &entr->col_chain;
386 m->last_row_el[row] = &entr->row_chain;
391 void matrix_set_row_bulk(sp_matrix_t *m, int row, int *cols, int num_cols, double val)
393 matrix_elem_t *me = NULL;
396 sp_matrix_list_head_t *leftof, *above;
397 sp_matrix_list_head_t *prev_leftof, *prev_above;
399 /* if necessary enlarge the matrix */
400 if (row > m->maxrow) {
403 m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
405 if (cols[num_cols - 1] > m->maxcol) {
406 m->maxcol = cols[num_cols - 1];
407 if (cols[num_cols - 1] >= m->colc)
408 m_alloc_col(m, m->colc, m_new_size(m->colc, cols[num_cols - 1]));
411 /* set start values */
415 for (i = 0; i < num_cols; ++i) {
416 /* search for existing entry */
417 me = m_search_in_row(m, row, cols[i], &leftof, &prev_leftof);
419 /* if it exists, set the value and return */
422 me->val = (float)val;
424 entr = _container_of(me, entry_t, e);
426 /* remove row_chain entry */
428 prev_leftof->next = entr->row_chain.next;
430 m->rows[row]->next = entr->row_chain.next;
432 /* remove col_chain entry */
434 prev_above->next = entr->col_chain.next;
436 m->cols[cols[i]]->next = entr->col_chain.next;
438 entr->row_chain.next = NULL;
439 entr->col_chain.next = NULL;
441 /* set the last pointer to the "previous" element */
442 if (m->last_col_el[cols[i]] == &entr->col_chain ||
443 m->last_row_el[row] == &entr->row_chain)
445 m->last_col_el[cols[i]] = prev_above ? prev_above : m->cols[cols[i]];
446 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
456 /* if it does not exist and 0 should be set just quit */
460 /* we have to search the col list as well, to get the above pointer */
461 m_search_in_col(m, row, cols[i], &above, &prev_above);
463 /* now leftof and above are the entry_t's prior the new one in each direction */
465 /* insert new entry */
466 entr = XMALLOC(entry_t);
468 entr->e.col = cols[i];
469 entr->e.val = (float)val;
471 m->last_col_el[cols[i]] = &entr->col_chain;
472 m->last_row_el[row] = &entr->row_chain;
474 /* add row_chain entry */
475 entr->row_chain.next = leftof->next;
476 leftof->next = &entr->row_chain;
478 /* add col_chain entry */
479 entr->col_chain.next = above->next;
480 above->next = &entr->col_chain;
486 double matrix_get(const sp_matrix_t *m, int row, int col)
488 sp_matrix_list_head_t *dummy, *dummy2;
491 if (is_empty_row(row) || is_empty_col(col))
494 if (m->maxrow < m->maxcol)
495 me = m_search_in_col(m, row, col, &dummy, &dummy2);
497 me = m_search_in_row(m, row, col, &dummy, &dummy2);
500 assert(me->col == col && me->row == row);
502 return me ? me->val : 0.0;
505 int matrix_get_entries(const sp_matrix_t *m)
510 int matrix_get_rowcount(const sp_matrix_t *m)
512 return m->maxrow + 1;
515 int matrix_get_colcount(const sp_matrix_t *m)
517 return m->maxcol + 1;
520 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row)
522 if (is_empty_row(row))
526 m->first = m->rows[row];
527 m->last = m->first->next;
528 m->next = m->last ? m->last->next : NULL;
530 assert (list_entry_by_row(m->last)->row == row);
532 return list_entry_by_row(m->last);
535 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col)
537 if (is_empty_col(col))
541 m->first = m->cols[col];
542 m->last = m->first->next;
543 m->next = m->last ? m->last->next : NULL;
545 assert (list_entry_by_col(m->last)->col == col);
547 return list_entry_by_col(m->last);
550 static inline const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow)
552 const matrix_elem_t *res;
555 for (i = startrow; i <= m->maxrow; ++i) {
556 res = matrix_row_first(m, i);
567 const matrix_elem_t *matrix_first(sp_matrix_t *m)
569 return matrix_first_from(m, 0);
572 const matrix_elem_t *matrix_next(sp_matrix_t *m)
574 assert(m->first && "Start iteration with matrix_???_first, before calling me!");
576 if (m->next == NULL) {
578 return matrix_first_from(m, ++m->iter_row);
584 m->next = m->next->next;
587 return list_entry_by_col(m->last);
588 else /* right or all */
589 return list_entry_by_row(m->last);
592 static int cmp_count(const void *e1, const void *e2)
594 return (int *)e2 - (int *)e1;
597 static inline void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow)
599 const matrix_elem_t *e;
600 bitset_set(fullrow, row);
601 matrix_foreach_in_col(m, row, e) {
602 if (! bitset_is_set(fullrow, e->row)) {
603 assert(0.0 == matrix_get(m, e->col, e->row));
604 matrix_set(m, e->col, e->row, e->val);
605 matrix_set(m, e->row, e->col, 0.0);
610 void matrix_optimize(sp_matrix_t *m)
614 const matrix_elem_t *e;
617 size = MAX(m->maxcol, m->maxrow)+1;
619 /* kill all double-entries (Mij and Mji are set) */
620 matrix_foreach(m, e) {
623 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will always have the same color as root");
624 t_val = matrix_get(m, e->col, e->row);
625 if (fabs(t_val) > 1e-10) {
626 matrix_set(m, e->col, e->row, 0);
627 matrix_set(m, e->row, e->col, e->val + t_val);
631 c = alloca(size * sizeof(*c));
633 fullrow = bitset_alloca(size);
635 /* kill 'all' rows containing only 1 entry */
638 /* count elements in rows */
639 memset(c, 0, size * sizeof(*c));
644 for (i = 0; i<size; ++i)
645 if (c[i] == 1 && ! bitset_is_set(fullrow, i)) {
647 /* if the other row isn't empty move the e in there, else fill e's row */
648 if (e = matrix_row_first(m, i), e) {
650 matrix_fill_row(m, e->col, fullrow);
652 matrix_fill_row(m, e->row, fullrow);
658 memset(c, 0, size * sizeof(*c));
662 qsort(c, size, sizeof(*c), cmp_count);
664 for (i = 0; i < size; ++i) {
665 if (! bitset_is_set(fullrow, i))
666 matrix_fill_row(m, i, fullrow);
670 void matrix_dump(sp_matrix_t *m, FILE *out, int factor)
673 const matrix_elem_t *e;
675 for (i = 0; i <= m->maxrow; ++i) {
677 matrix_foreach_in_row(m, i, e) {
678 for (o = last_idx + 1; o < e->col; ++o)
679 fprintf(out, " %4.1f" , 0.0);
681 fprintf(out, " %4.1f", factor * e->val);
685 for (o = last_idx + 1; o <= m->maxcol; ++o)
686 fprintf(out, " %4.1f" , 0.0);
692 void matrix_self_test(int d)
695 const matrix_elem_t *e;
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);
743 for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
746 matrix_set(m, 1,1,0);
747 assert(5 == matrix_get_entries(m));