nsz Git - libfirm/blob - ir/lpp/sp_matrix.c

   1 /**
   2  * Author:      Daniel Grund, Christian Wuerdig
   3  * Date:        07.04.2005
   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 $
   7  *
   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).
  12  */
  13
  14 #ifdef _WIN32
  15 #include <malloc.h>
  16 #endif
  17 #ifdef HAVE_ALLOCA_H
  18 #include <alloca.h>
  19 #endif
  20
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <assert.h>
  24 #include <math.h>
  25
  26 #include "sp_matrix.h"
  27
  28 #include "irtools.h"
  29 #include "bitset.h"
  30 #include "xmalloc.h"
  31
  32 typedef enum iter_direction_t {
  33         down, right, all
  34 } iter_direction_t;
  35
  36 /**
  37  * Embedded list pointer.
  38  */
  39 typedef struct sp_matrix_list_head_t {
  40         struct sp_matrix_list_head_t *next;
  41 } sp_matrix_list_head_t;
  42
  43 /**
  44  * A matrix entry.
  45  */
  46 typedef struct entry_t {
  47         sp_matrix_list_head_t col_chain; /**< points to next element in same column */
  48         sp_matrix_list_head_t row_chain; /**< points to next element in same row */
  49         matrix_elem_t         e;         /**< The actual element */
  50 } entry_t;
  51
  52 struct sp_matrix_t {
  53         /* These specify the dimensions of the matrix.
  54          * They equal the largest values ever used in matrix_set */
  55         int maxrow, maxcol;
  56         /* These are the dimensions of allocated arrays below.
  57          * rowc >= maxrow and colc >= maxcol hold. */
  58         int rowc, colc;
  59         /* number of entries */
  60         int entries;
  61         /* arrays of sp_matrix_list_head* as entry-points to rows and cols */
  62         sp_matrix_list_head_t **rows, **cols;
  63         /* for iteration: first is to remember start-point;
  64          *                last was returned just before
  65          *                next is used in case last was removed from list */
  66         iter_direction_t dir;
  67         sp_matrix_list_head_t *first, *last, *next;
  68         int iter_row; /* used for iteration over all elements */
  69         /* for each column the last inserted element in col list */
  70         sp_matrix_list_head_t **last_col_el;
  71         /* for each row the last inserted element in row list */
  72         sp_matrix_list_head_t **last_row_el;
  73 };
  74
  75 #define SP_MATRIX_INIT_LIST_HEAD(ptr) do { (ptr)->next = NULL; } while (0)
  76
  77 #define _offsetof(type,member) ((char *) &(((type *) 0)->member) - (char *) 0)
  78 #define _container_of(ptr,type,member) ((type *) ((char *) (ptr) - _offsetof(type, member)))
  79
  80 #define is_empty_row(row) (row > m->maxrow || (m->rows[row])->next == NULL)
  81 #define is_empty_col(col) (col > m->maxcol || (m->cols[col])->next == NULL)
  82
  83 #define list_entry_by_col(h) (&_container_of(h, entry_t, col_chain)->e)
  84 #define list_entry_by_row(h) (&_container_of(h, entry_t, row_chain)->e)
  85
  86 /**
  87  * Returns the size of a single matrix element.
  88  */
  89 unsigned matrix_get_elem_size(void)
  90 {
  91         return sizeof(entry_t);
  92 }
  93
  94 /**
  95  * Returns the new size for an array of size old_size,
  96  *  which must at least store an entry at position min.
  97  */
  98 static inline int m_new_size(int old_size, int min)
  99 {
 100         unsigned bits = 0;
 101         assert(min >= old_size);
 102         while (min > 0) {
 103                 min >>= 1;
 104                 bits++;
 105         }
 106         assert(bits < sizeof(min) * 8 - 1);
 107         return 1 << bits;
 108 }
 109
 110 /**
 111  * Allocates space for @p count entries in the rows array and
 112  * initializes all entries from @p start to the end.
 113  */
 114 static inline void m_alloc_row(sp_matrix_t *m, int start, int count)
 115 {
 116         int p;
 117
 118         m->rowc        = count;
 119         m->rows        = XREALLOC(m->rows, sp_matrix_list_head_t *, m->rowc);
 120         m->last_row_el = XREALLOC(m->last_row_el, sp_matrix_list_head_t *, m->rowc);
 121
 122         for (p = start; p < m->rowc; ++p) {
 123                 m->rows[p] = XMALLOC(sp_matrix_list_head_t);
 124                 SP_MATRIX_INIT_LIST_HEAD(m->rows[p]);
 125                 m->last_row_el[p] = m->rows[p];
 126         }
 127 }
 128
 129 /**
 130  * Allocates space for @p count entries in the cols array and
 131  * initializes all entries from @p start to the end.
 132  */
 133 static inline void m_alloc_col(sp_matrix_t *m, int start, int count)
 134 {
 135         int p;
 136
 137         m->colc        = count;
 138         m->cols        = XREALLOC(m->cols, sp_matrix_list_head_t*, m->colc);
 139         m->last_col_el = XREALLOC(m->last_col_el, sp_matrix_list_head_t*, m->colc);
 140
 141         for (p = start; p < m->colc; ++p) {
 142                 m->cols[p] = XMALLOC(sp_matrix_list_head_t);
 143                 SP_MATRIX_INIT_LIST_HEAD(m->cols[p]);
 144                 m->last_col_el[p] = m->cols[p];
 145         }
 146 }
 147
 148 /**
 149  * Searches in row @p row for the matrix element m[row, col], starting at element @p start.
 150  * @return If the element exists:
 151  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 152  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 153  *         @p prev_prev always points to the previous element of @p prev
 154  */
 155 static inline matrix_elem_t *m_search_in_row_from(const sp_matrix_t *m,
 156                         int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 157 {
 158         sp_matrix_list_head_t *row_start;
 159         matrix_elem_t         *res = NULL;
 160
 161         row_start = m->rows[row];
 162         *prev     = start;
 163
 164         while ((*prev)->next != NULL && list_entry_by_row((*prev)->next)->col <= col) {
 165                 (*prev_prev) = (*prev);
 166                 *prev        = (*prev)->next;
 167         }
 168
 169         if (*prev != row_start) {
 170                 matrix_elem_t *me = list_entry_by_row(*prev);
 171
 172                 if (me->row == row && me->col == col)
 173                         res = me;
 174         }
 175
 176         if (res) {
 177                 m->last_row_el[row] = *prev;
 178         }
 179
 180         return res;
 181 }
 182
 183 /**
 184  * Searches in row @p row for the matrix element m[row, col].
 185  * @return If the element exists:
 186  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 187  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 188  *         @p prev_prev always points to the previous element of @p prev
 189  */
 190 static inline matrix_elem_t *m_search_in_row(const sp_matrix_t *m,
 191                         int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 192 {
 193         sp_matrix_list_head_t *start = m->rows[row];
 194
 195         *prev_prev = NULL;
 196
 197         if (m->last_row_el[row] != start) {
 198                 matrix_elem_t *el = list_entry_by_row(m->last_row_el[row]);
 199                 if (el->col < col) {
 200                         *prev_prev = start = m->last_row_el[row];
 201                 }
 202         }
 203
 204         return m_search_in_row_from(m, row, col, start, prev, prev_prev);
 205 }
 206
 207 /**
 208  * Searches in col @p col for the matrix element m[row, col], starting at @p start.
 209  * @return If the element exists:
 210  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 211  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 212  *         @p prev_prev always points to the previous element of @p prev
 213  */
 214 static inline matrix_elem_t *m_search_in_col_from(const sp_matrix_t *m,
 215                         int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 216 {
 217         sp_matrix_list_head_t *col_start;
 218         matrix_elem_t         *res = NULL;
 219
 220         col_start = m->cols[col];
 221         *prev     = start;
 222
 223         while ((*prev)->next != NULL && list_entry_by_col((*prev)->next)->row <= row) {
 224                 *prev_prev = (*prev);
 225                 *prev      = (*prev)->next;
 226         }
 227
 228         if (*prev != col_start) {
 229                 matrix_elem_t *me = list_entry_by_col(*prev);
 230
 231                 if (me->row == row && me->col == col)
 232                         res = me;
 233         }
 234
 235         if (res) {
 236                 m->last_col_el[col] = *prev;
 237         }
 238
 239         return res;
 240 }
 241
 242 /**
 243  * Searches in col @p col for the matrix element m[row, col].
 244  * @return If the element exists:
 245  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 246  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 247  *         @p prev_prev always points to the previous element of @p prev
 248  */
 249 static inline matrix_elem_t *m_search_in_col(const sp_matrix_t *m,
 250                         int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 251 {
 252         sp_matrix_list_head_t *start = m->cols[col];
 253
 254         *prev_prev = NULL;
 255
 256         if (m->last_col_el[col] != start) {
 257                 matrix_elem_t *el = list_entry_by_col(m->last_col_el[col]);
 258                 if (el->row < row) {
 259                         *prev_prev = start = m->last_col_el[col];
 260                 }
 261         }
 262
 263         return m_search_in_col_from(m, row, col, start, prev, prev_prev);
 264 }
 265
 266 sp_matrix_t *new_matrix(int row_init, int col_init)
 267 {
 268         sp_matrix_t *res = XMALLOCZ(sp_matrix_t);
 269         res->maxrow = -1;
 270         res->maxcol = -1;
 271         m_alloc_row(res, 0, MAX(0, row_init));
 272         m_alloc_col(res, 0, MAX(0, col_init));
 273         return res;
 274 }
 275
 276 void del_matrix(sp_matrix_t *m)
 277 {
 278         int i;
 279
 280         for (i = 0; i < m->rowc; ++i) {
 281                 if (! is_empty_row(i)) {
 282                         entry_t *e;
 283                         sp_matrix_list_head_t *n;
 284
 285                         n = m->rows[i]->next;
 286                         do {
 287                                 /* get current matrix element */
 288                                 e = _container_of(n, entry_t, row_chain);
 289                                 n = n->next;
 290                                 xfree(e);
 291                         } while (n != NULL);
 292
 293                 }
 294                 xfree(m->rows[i]);
 295         }
 296         for (i = 0; i < m->colc; ++i)
 297                 xfree(m->cols[i]);
 298         xfree(m->last_col_el);
 299         xfree(m->last_row_el);
 300         xfree(m->rows);
 301         xfree(m->cols);
 302         xfree(m);
 303 }
 304
 305 void matrix_set(sp_matrix_t *m, int row, int col, double val)
 306 {
 307         matrix_elem_t *me = NULL;
 308         entry_t       *entr;
 309         sp_matrix_list_head_t *leftof, *above;
 310         sp_matrix_list_head_t *prev_leftof, *prev_above;
 311
 312         /* if necessary enlarge the matrix */
 313         if (row > m->maxrow) {
 314                 m->maxrow = row;
 315                 if (row >= m->rowc)
 316                         m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
 317         }
 318         if (col > m->maxcol) {
 319                 m->maxcol = col;
 320                 if (col >= m->colc)
 321                         m_alloc_col(m, m->colc, m_new_size(m->colc, col));
 322         }
 323
 324         /* search for existing entry */
 325         if (m->maxrow < m->maxcol)
 326                 me = m_search_in_col(m, row, col, &above, &prev_above);
 327         else
 328                 me = m_search_in_row(m, row, col, &leftof, &prev_leftof);
 329
 330         /* if it exists, set the value and return */
 331         if (me) {
 332                 if (val != 0) {
 333                         me->val = (float)val;
 334                 } else {
 335                         entr = _container_of(me, entry_t, e);
 336
 337                         /* remove row_chain entry */
 338                         if (prev_leftof)
 339                                 prev_leftof->next = entr->row_chain.next;
 340                         else
 341                                 m->rows[row]->next = entr->row_chain.next;
 342
 343                         /* remove col_chain entry */
 344                         if (prev_above)
 345                                 prev_above->next = entr->col_chain.next;
 346                         else
 347                                 m->cols[col]->next = entr->col_chain.next;
 348
 349                         entr->row_chain.next = NULL;
 350                         entr->col_chain.next = NULL;
 351
 352                         /* set the last pointer to the "previous" element */
 353                         if (m->last_col_el[col] == &entr->col_chain ||
 354                                 m->last_row_el[row] == &entr->row_chain)
 355                         {
 356                                 m->last_col_el[col] = prev_above  ? prev_above  : m->cols[col];
 357                                 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
 358                         }
 359
 360                         free(entr);
 361                         m->entries--;
 362                 }
 363                 return;
 364         }
 365
 366         /* if it does not exist and 0 should be set just quit */
 367         if (val == 0)
 368                 return;
 369
 370         /* if it does not exist and val != 0 search the other direction */
 371         if (m->maxrow >= m->maxcol)
 372                 m_search_in_col(m, row, col, &above, &prev_above);
 373         else
 374                 m_search_in_row(m, row, col, &leftof, &prev_leftof);
 375         /* now leftof and above are the entry_t's prior the new one in each direction */
 376
 377         /* insert new entry */
 378         entr        = XMALLOC(entry_t);
 379         entr->e.row = row;
 380         entr->e.col = col;
 381         entr->e.val = (float)val;
 382
 383         /* add row_chain entry */
 384         entr->row_chain.next = leftof->next;
 385         leftof->next = &entr->row_chain;
 386
 387         /* add col_chain entry */
 388         entr->col_chain.next = above->next;
 389         above->next = &entr->col_chain;
 390
 391         m->last_col_el[col] = &entr->col_chain;
 392         m->last_row_el[row] = &entr->row_chain;
 393
 394         m->entries++;
 395 }
 396
 397 void matrix_set_row_bulk(sp_matrix_t *m, int row, int *cols, int num_cols, double val)
 398 {
 399         matrix_elem_t *me = NULL;
 400         entry_t       *entr;
 401         int           i;
 402         sp_matrix_list_head_t *leftof, *above;
 403         sp_matrix_list_head_t *prev_leftof, *prev_above;
 404
 405         /* if necessary enlarge the matrix */
 406         if (row > m->maxrow) {
 407                 m->maxrow = row;
 408                 if (row >= m->rowc)
 409                         m_alloc_row(m, m->rowc, m_new_size(m->rowc, row));
 410         }
 411         if (cols[num_cols - 1] > m->maxcol) {
 412                 m->maxcol = cols[num_cols - 1];
 413                 if (cols[num_cols - 1] >= m->colc)
 414                         m_alloc_col(m, m->colc, m_new_size(m->colc, cols[num_cols - 1]));
 415         }
 416
 417     /* set start values */
 418         prev_above  = NULL;
 419         prev_leftof = NULL;
 420
 421         for (i = 0; i < num_cols; ++i) {
 422                 /* search for existing entry */
 423                 me = m_search_in_row(m, row, cols[i], &leftof, &prev_leftof);
 424
 425                 /* if it exists, set the value and return */
 426                 if (me) {
 427                         if (val != 0) {
 428                                 me->val = (float)val;
 429                         } else {
 430                                 entr = _container_of(me, entry_t, e);
 431
 432                                 /* remove row_chain entry */
 433                                 if (prev_leftof)
 434                                         prev_leftof->next = entr->row_chain.next;
 435                                 else
 436                                         m->rows[row]->next = entr->row_chain.next;
 437
 438                                 /* remove col_chain entry */
 439                                 if (prev_above)
 440                                         prev_above->next = entr->col_chain.next;
 441                                 else
 442                                         m->cols[cols[i]]->next = entr->col_chain.next;
 443
 444                                 entr->row_chain.next = NULL;
 445                                 entr->col_chain.next = NULL;
 446
 447                                 /* set the last pointer to the "previous" element */
 448                                 if (m->last_col_el[cols[i]] == &entr->col_chain ||
 449                                         m->last_row_el[row]     == &entr->row_chain)
 450                                 {
 451                                         m->last_col_el[cols[i]] = prev_above  ? prev_above  : m->cols[cols[i]];
 452                                         m->last_row_el[row]     = prev_leftof ? prev_leftof : m->rows[row];
 453                                 }
 454
 455                                 free(entr);
 456                                 m->entries--;
 457                         }
 458
 459                         continue;
 460                 }
 461
 462                 /* if it does not exist and 0 should be set just quit */
 463                 if (val == 0)
 464                         continue;
 465
 466                 /* we have to search the col list as well, to get the above pointer */
 467                 m_search_in_col(m, row, cols[i], &above, &prev_above);
 468
 469                 /* now leftof and above are the entry_t's prior the new one in each direction */
 470
 471                 /* insert new entry */
 472                 entr        = XMALLOC(entry_t);
 473                 entr->e.row = row;
 474                 entr->e.col = cols[i];
 475                 entr->e.val = (float)val;
 476
 477                 m->last_col_el[cols[i]] = &entr->col_chain;
 478                 m->last_row_el[row]     = &entr->row_chain;
 479
 480                 /* add row_chain entry */
 481                 entr->row_chain.next = leftof->next;
 482                 leftof->next = &entr->row_chain;
 483
 484                 /* add col_chain entry */
 485                 entr->col_chain.next = above->next;
 486                 above->next = &entr->col_chain;
 487
 488                 m->entries++;
 489         }
 490 }
 491
 492 double matrix_get(const sp_matrix_t *m, int row, int col)
 493 {
 494         sp_matrix_list_head_t *dummy, *dummy2;
 495         matrix_elem_t *me;
 496
 497         if (is_empty_row(row) || is_empty_col(col))
 498                 return 0.0;
 499
 500         if (m->maxrow < m->maxcol)
 501                 me = m_search_in_col(m, row, col, &dummy, &dummy2);
 502         else
 503                 me = m_search_in_row(m, row, col, &dummy, &dummy2);
 504
 505         if (me)
 506                 assert(me->col == col && me->row == row);
 507
 508         return me ? me->val : 0.0;
 509 }
 510
 511 int matrix_get_entries(const sp_matrix_t *m)
 512 {
 513         return m->entries;
 514 }
 515
 516 int matrix_get_rowcount(const sp_matrix_t *m)
 517 {
 518         return m->maxrow + 1;
 519 }
 520
 521 int matrix_get_colcount(const sp_matrix_t *m)
 522 {
 523         return m->maxcol + 1;
 524 }
 525
 526 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row)
 527 {
 528         if (is_empty_row(row))
 529                 return NULL;
 530
 531         m->dir   = right;
 532         m->first = m->rows[row];
 533         m->last  = m->first->next;
 534         m->next  = m->last ? m->last->next : NULL;
 535
 536         assert (list_entry_by_row(m->last)->row == row);
 537
 538         return list_entry_by_row(m->last);
 539 }
 540
 541 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col)
 542 {
 543         if (is_empty_col(col))
 544                 return NULL;
 545
 546         m->dir   = down;
 547         m->first = m->cols[col];
 548         m->last  = m->first->next;
 549         m->next  = m->last ? m->last->next : NULL;
 550
 551         assert (list_entry_by_col(m->last)->col == col);
 552
 553         return list_entry_by_col(m->last);
 554 }
 555
 556 static inline const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow)
 557 {
 558         const matrix_elem_t *res;
 559         int i;
 560
 561         for (i = startrow; i <= m->maxrow; ++i) {
 562                 res = matrix_row_first(m, i);
 563                 if (res) {
 564                         m->iter_row = i;
 565                         m->dir      = all;
 566                         return res;
 567                 }
 568         }
 569
 570         return NULL;
 571 }
 572
 573 const matrix_elem_t *matrix_first(sp_matrix_t *m)
 574 {
 575         return matrix_first_from(m, 0);
 576 }
 577
 578 const matrix_elem_t *matrix_next(sp_matrix_t *m)
 579 {
 580         assert(m->first && "Start iteration with matrix_???_first, before calling me!");
 581
 582         if (m->next == NULL) {
 583                 if (m->dir == all)
 584                         return matrix_first_from(m, ++m->iter_row);
 585                 else
 586                         return NULL;
 587         }
 588
 589         m->last = m->next;
 590         m->next = m->next->next;
 591
 592         if (m->dir == down)
 593                 return list_entry_by_col(m->last);
 594         else /* right or all */
 595                 return list_entry_by_row(m->last);
 596 }
 597
 598 static int cmp_count(const void *e1, const void *e2)
 599 {
 600         return (int *)e2 - (int *)e1;
 601 }
 602
 603 static inline void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow)
 604 {
 605         const matrix_elem_t *e;
 606         bitset_set(fullrow, row);
 607         matrix_foreach_in_col(m, row, e) {
 608                 if (! bitset_is_set(fullrow, e->row)) {
 609                         assert(0.0 == matrix_get(m, e->col, e->row));
 610                         matrix_set(m, e->col, e->row, e->val);
 611                         matrix_set(m, e->row, e->col, 0.0);
 612                 }
 613         }
 614 }
 615
 616 void matrix_optimize(sp_matrix_t *m)
 617 {
 618         int i, size, redo;
 619         int *c;
 620         const matrix_elem_t *e;
 621         bitset_t *fullrow;
 622
 623         size = MAX(m->maxcol, m->maxrow)+1;
 624
 625         /* kill all double-entries (Mij and Mji are set) */
 626         matrix_foreach(m, e) {
 627                 double t_val;
 628
 629                 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will always have the same color as root");
 630                 t_val = matrix_get(m, e->col, e->row);
 631                 if (fabs(t_val) > 1e-10) {
 632                         matrix_set(m, e->col, e->row, 0);
 633                         matrix_set(m, e->row, e->col, e->val + t_val);
 634                 }
 635         }
 636
 637         c       = alloca(size * sizeof(*c));
 638         redo    = 1;
 639         fullrow = bitset_alloca(size);
 640
 641         /* kill 'all' rows containing only 1 entry */
 642         while (redo) {
 643                 redo = 0;
 644                 /* count elements in rows */
 645                 memset(c, 0, size * sizeof(*c));
 646
 647                 matrix_foreach(m, e)
 648                         c[e->row]++;
 649
 650                 for (i = 0; i<size; ++i)
 651                         if (c[i] == 1 && ! bitset_is_set(fullrow, i)) {
 652                                 redo = 1;
 653                                 /* if the other row isn't empty move the e in there, else fill e's row */
 654                                 if (e = matrix_row_first(m, i), e) {
 655                                         if (c[e->col] > 0)
 656                                                 matrix_fill_row(m, e->col, fullrow);
 657                                         else
 658                                                 matrix_fill_row(m, e->row, fullrow);
 659                                 }
 660                         }
 661         }
 662
 663
 664         memset(c, 0, size * sizeof(*c));
 665         matrix_foreach(m, e)
 666                 c[e->row]++;
 667
 668         qsort(c, size, sizeof(*c), cmp_count);
 669
 670         for (i = 0; i < size; ++i) {
 671                 if (! bitset_is_set(fullrow, i))
 672                         matrix_fill_row(m, i, fullrow);
 673         }
 674 }
 675
 676 void matrix_dump(sp_matrix_t *m, FILE *out, int factor)
 677 {
 678         int i, o, last_idx;
 679         const matrix_elem_t *e;
 680
 681         for (i = 0; i <= m->maxrow; ++i) {
 682                 last_idx = -1;
 683                 matrix_foreach_in_row(m, i, e) {
 684                         for (o = last_idx + 1; o < e->col; ++o)
 685                                 fprintf(out, " %4.1f" , 0.0);
 686
 687                         fprintf(out, " %4.1f", factor * e->val);
 688                         last_idx = e->col;
 689                 }
 690
 691                 for (o = last_idx + 1; o <= m->maxcol; ++o)
 692                         fprintf(out, " %4.1f" , 0.0);
 693
 694                 fprintf(out, "\n");
 695         }
 696 }
 697
 698 void matrix_self_test(int d)
 699 {
 700         int i, o;
 701         const matrix_elem_t *e;
 702         sp_matrix_t *m = new_matrix(10, 10);
 703
 704         for (i = 0; i < d; ++i)
 705                 for (o = 0; o < d; ++o)
 706                         matrix_set(m, i, o, i*o);
 707
 708         for (i = 0; i < d; ++i)
 709                 for (o = 0; o<d; ++o)
 710                         assert(matrix_get(m, i, o) == i*o);
 711
 712         i = 1;
 713         matrix_foreach_in_row(m,1,e) {
 714                 assert(e->val == i);
 715                 i++;
 716         }
 717         assert(!matrix_next(m)); /*iter must finish */
 718
 719         i = d-1;
 720         matrix_foreach_in_col(m,d-1,e) {
 721                 assert(e->val == i);
 722                 i += d-1;
 723         }
 724         assert(!matrix_next(m));
 725         del_matrix(m);
 726         m = new_matrix(16,16);
 727         matrix_set(m, 1,1,9);
 728         matrix_set(m, 1,2,8);
 729         matrix_set(m, 1,3,7);
 730
 731         matrix_set(m, 1,3,6);
 732         matrix_set(m, 1,2,5);
 733         matrix_set(m, 1,1,4);
 734         i = 1;
 735         matrix_foreach_in_row(m, 1, e) {
 736                 assert(e->row == 1 && e->col == i && e->val == i+3);
 737                 i++;
 738         }
 739         assert(i == 4);
 740         del_matrix(m);
 741
 742         m = new_matrix(5,5);
 743         matrix_set(m, 1,1,1);
 744         matrix_set(m, 2,2,2);
 745         matrix_set(m, 3,3,3);
 746         matrix_set(m, 3,5,4);
 747         matrix_set(m, 4,4,5);
 748         matrix_set(m, 5,5,6);
 749         for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
 750                 assert(e->val == i);
 751         assert(i == 7);
 752         matrix_set(m, 1,1,0);
 753         assert(5 == matrix_get_entries(m));
 754         del_matrix(m);
 755 }