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