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

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