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         return sizeof(entry_t);
  91 }
  92
  93 /**
  94  * Returns the new size for an array of size old_size,
  95  *  which must at least store an entry at position min.
  96  */
  97 static inline int _m_new_size(int old_size, int min) {
  98         unsigned bits = 0;
  99         assert(min >= old_size);
 100         while (min > 0) {
 101                 min >>= 1;
 102                 bits++;
 103         }
 104         assert(bits < sizeof(min) * 8 - 1);
 105         return 1 << bits;
 106 }
 107
 108 /**
 109  * Allocates space for @p count entries in the rows array and
 110  * initializes all entries from @p start to the end.
 111  */
 112 static inline void _m_alloc_row(sp_matrix_t *m, int start, int count) {
 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         int p;
 132
 133         m->colc        = count;
 134         m->cols        = XREALLOC(m->cols, sp_matrix_list_head_t*, m->colc);
 135         m->last_col_el = XREALLOC(m->last_col_el, sp_matrix_list_head_t*, m->colc);
 136
 137         for (p = start; p < m->colc; ++p) {
 138                 m->cols[p] = XMALLOC(sp_matrix_list_head_t);
 139                 SP_MATRIX_INIT_LIST_HEAD(m->cols[p]);
 140                 m->last_col_el[p] = m->cols[p];
 141         }
 142 }
 143
 144 /**
 145  * Searches in row @p row for the matrix element m[row, col], starting at element @p start.
 146  * @return If the element exists:
 147  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 148  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 149  *         @p prev_prev always points to the previous element of @p prev
 150  */
 151 static inline matrix_elem_t *_m_search_in_row_from(const sp_matrix_t *m,
 152                         int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 153 {
 154         sp_matrix_list_head_t *row_start;
 155         matrix_elem_t         *res = NULL;
 156
 157         row_start = m->rows[row];
 158         *prev     = start;
 159
 160         while ((*prev)->next != NULL && list_entry_by_row((*prev)->next)->col <= col) {
 161                 (*prev_prev) = (*prev);
 162                 *prev        = (*prev)->next;
 163         }
 164
 165         if (*prev != row_start) {
 166                 matrix_elem_t *me = list_entry_by_row(*prev);
 167
 168                 if (me->row == row && me->col == col)
 169                         res = me;
 170         }
 171
 172         if (res) {
 173                 m->last_row_el[row] = *prev;
 174         }
 175
 176         return res;
 177 }
 178
 179 /**
 180  * Searches in row @p row for the matrix element m[row, col].
 181  * @return If the element exists:
 182  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 183  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 184  *         @p prev_prev always points to the previous element of @p prev
 185  */
 186 static inline matrix_elem_t *_m_search_in_row(const sp_matrix_t *m,
 187                         int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 188 {
 189         sp_matrix_list_head_t *start = m->rows[row];
 190
 191         *prev_prev = NULL;
 192
 193         if (m->last_row_el[row] != start) {
 194                 matrix_elem_t *el = list_entry_by_row(m->last_row_el[row]);
 195                 if (el->col < col) {
 196                         *prev_prev = start = m->last_row_el[row];
 197                 }
 198         }
 199
 200         return _m_search_in_row_from(m, row, col, start, prev, prev_prev);
 201 }
 202
 203 /**
 204  * Searches in col @p col for the matrix element m[row, col], starting at @p start.
 205  * @return If the element exists:
 206  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 207  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 208  *         @p prev_prev always points to the previous element of @p prev
 209  */
 210 static inline matrix_elem_t *_m_search_in_col_from(const sp_matrix_t *m,
 211                         int row, int col, sp_matrix_list_head_t *start, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 212 {
 213         sp_matrix_list_head_t *col_start;
 214         matrix_elem_t         *res = NULL;
 215
 216         col_start = m->cols[col];
 217         *prev     = start;
 218
 219         while ((*prev)->next != NULL && list_entry_by_col((*prev)->next)->row <= row) {
 220                 *prev_prev = (*prev);
 221                 *prev      = (*prev)->next;
 222         }
 223
 224         if (*prev != col_start) {
 225                 matrix_elem_t *me = list_entry_by_col(*prev);
 226
 227                 if (me->row == row && me->col == col)
 228                         res = me;
 229         }
 230
 231         if (res) {
 232                 m->last_col_el[col] = *prev;
 233         }
 234
 235         return res;
 236 }
 237
 238 /**
 239  * Searches in col @p col for the matrix element m[row, col].
 240  * @return If the element exists:
 241  *            Element m[row, col] and @p prev points to the sp_matrix_list_head in the entry_t holding the element.
 242  *         Else: NULL and @p prev points to the sp_matrix_list_head after which the element would be inserted.
 243  *         @p prev_prev always points to the previous element of @p prev
 244  */
 245 static inline matrix_elem_t *_m_search_in_col(const sp_matrix_t *m,
 246                         int row, int col, sp_matrix_list_head_t **prev, sp_matrix_list_head_t **prev_prev)
 247 {
 248         sp_matrix_list_head_t *start = m->cols[col];
 249
 250         *prev_prev = NULL;
 251
 252         if (m->last_col_el[col] != start) {
 253                 matrix_elem_t *el = list_entry_by_col(m->last_col_el[col]);
 254                 if (el->row < row) {
 255                         *prev_prev = start = m->last_col_el[col];
 256                 }
 257         }
 258
 259         return _m_search_in_col_from(m, row, col, start, prev, prev_prev);
 260 }
 261
 262 sp_matrix_t *new_matrix(int row_init, int col_init) {
 263         sp_matrix_t *res = XMALLOCZ(sp_matrix_t);
 264         res->maxrow = -1;
 265         res->maxcol = -1;
 266         _m_alloc_row(res, 0, MAX(0, row_init));
 267         _m_alloc_col(res, 0, MAX(0, col_init));
 268         return res;
 269 }
 270
 271 void del_matrix(sp_matrix_t *m) {
 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         matrix_elem_t *me = NULL;
 301         entry_t       *entr;
 302         sp_matrix_list_head_t *leftof, *above;
 303         sp_matrix_list_head_t *prev_leftof, *prev_above;
 304
 305         /* if necessary enlarge the matrix */
 306         if (row > m->maxrow) {
 307                 m->maxrow = row;
 308                 if (row >= m->rowc)
 309                         _m_alloc_row(m, m->rowc, _m_new_size(m->rowc, row));
 310         }
 311         if (col > m->maxcol) {
 312                 m->maxcol = col;
 313                 if (col >= m->colc)
 314                         _m_alloc_col(m, m->colc, _m_new_size(m->colc, col));
 315         }
 316
 317         /* search for existing entry */
 318         if (m->maxrow < m->maxcol)
 319                 me = _m_search_in_col(m, row, col, &above, &prev_above);
 320         else
 321                 me = _m_search_in_row(m, row, col, &leftof, &prev_leftof);
 322
 323         /* if it exists, set the value and return */
 324         if (me) {
 325                 if (val != 0) {
 326                         me->val = (float)val;
 327                 } else {
 328                         entr = _container_of(me, entry_t, e);
 329
 330                         /* remove row_chain entry */
 331                         if (prev_leftof)
 332                                 prev_leftof->next = entr->row_chain.next;
 333                         else
 334                                 m->rows[row]->next = entr->row_chain.next;
 335
 336                         /* remove col_chain entry */
 337                         if (prev_above)
 338                                 prev_above->next = entr->col_chain.next;
 339                         else
 340                                 m->cols[col]->next = entr->col_chain.next;
 341
 342                         entr->row_chain.next = NULL;
 343                         entr->col_chain.next = NULL;
 344
 345                         /* set the last pointer to the "previous" element */
 346                         if (m->last_col_el[col] == &entr->col_chain ||
 347                                 m->last_row_el[row] == &entr->row_chain)
 348                         {
 349                                 m->last_col_el[col] = prev_above  ? prev_above  : m->cols[col];
 350                                 m->last_row_el[row] = prev_leftof ? prev_leftof : m->rows[row];
 351                         }
 352
 353                         free(entr);
 354                         m->entries--;
 355                 }
 356                 return;
 357         }
 358
 359         /* if it does not exist and 0 should be set just quit */
 360         if (val == 0)
 361                 return;
 362
 363         /* if it does not exist and val != 0 search the other direction */
 364         if (m->maxrow >= m->maxcol)
 365                 _m_search_in_col(m, row, col, &above, &prev_above);
 366         else
 367                 _m_search_in_row(m, row, col, &leftof, &prev_leftof);
 368         /* now leftof and above are the entry_t's prior the new one in each direction */
 369
 370         /* insert new entry */
 371         entr        = XMALLOC(entry_t);
 372         entr->e.row = row;
 373         entr->e.col = col;
 374         entr->e.val = (float)val;
 375
 376         /* add row_chain entry */
 377         entr->row_chain.next = leftof->next;
 378         leftof->next = &entr->row_chain;
 379
 380         /* add col_chain entry */
 381         entr->col_chain.next = above->next;
 382         above->next = &entr->col_chain;
 383
 384         m->last_col_el[col] = &entr->col_chain;
 385         m->last_row_el[row] = &entr->row_chain;
 386
 387         m->entries++;
 388 }
 389
 390 void matrix_set_row_bulk(sp_matrix_t *m, int row, int *cols, int num_cols, double val) {
 391         matrix_elem_t *me = NULL;
 392         entry_t       *entr;
 393         int           i;
 394         sp_matrix_list_head_t *leftof, *above;
 395         sp_matrix_list_head_t *prev_leftof, *prev_above;
 396
 397         /* if necessary enlarge the matrix */
 398         if (row > m->maxrow) {
 399                 m->maxrow = row;
 400                 if (row >= m->rowc)
 401                         _m_alloc_row(m, m->rowc, _m_new_size(m->rowc, row));
 402         }
 403         if (cols[num_cols - 1] > m->maxcol) {
 404                 m->maxcol = cols[num_cols - 1];
 405                 if (cols[num_cols - 1] >= m->colc)
 406                         _m_alloc_col(m, m->colc, _m_new_size(m->colc, cols[num_cols - 1]));
 407         }
 408
 409     /* set start values */
 410         prev_above  = NULL;
 411         prev_leftof = NULL;
 412
 413         for (i = 0; i < num_cols; ++i) {
 414                 /* search for existing entry */
 415                 me = _m_search_in_row(m, row, cols[i], &leftof, &prev_leftof);
 416
 417                 /* if it exists, set the value and return */
 418                 if (me) {
 419                         if (val != 0) {
 420                                 me->val = (float)val;
 421                         } else {
 422                                 entr = _container_of(me, entry_t, e);
 423
 424                                 /* remove row_chain entry */
 425                                 if (prev_leftof)
 426                                         prev_leftof->next = entr->row_chain.next;
 427                                 else
 428                                         m->rows[row]->next = entr->row_chain.next;
 429
 430                                 /* remove col_chain entry */
 431                                 if (prev_above)
 432                                         prev_above->next = entr->col_chain.next;
 433                                 else
 434                                         m->cols[cols[i]]->next = entr->col_chain.next;
 435
 436                                 entr->row_chain.next = NULL;
 437                                 entr->col_chain.next = NULL;
 438
 439                                 /* set the last pointer to the "previous" element */
 440                                 if (m->last_col_el[cols[i]] == &entr->col_chain ||
 441                                         m->last_row_el[row]     == &entr->row_chain)
 442                                 {
 443                                         m->last_col_el[cols[i]] = prev_above  ? prev_above  : m->cols[cols[i]];
 444                                         m->last_row_el[row]     = prev_leftof ? prev_leftof : m->rows[row];
 445                                 }
 446
 447                                 free(entr);
 448                                 m->entries--;
 449                         }
 450
 451                         continue;
 452                 }
 453
 454                 /* if it does not exist and 0 should be set just quit */
 455                 if (val == 0)
 456                         continue;
 457
 458                 /* we have to search the col list as well, to get the above pointer */
 459                 _m_search_in_col(m, row, cols[i], &above, &prev_above);
 460
 461                 /* now leftof and above are the entry_t's prior the new one in each direction */
 462
 463                 /* insert new entry */
 464                 entr        = XMALLOC(entry_t);
 465                 entr->e.row = row;
 466                 entr->e.col = cols[i];
 467                 entr->e.val = (float)val;
 468
 469                 m->last_col_el[cols[i]] = &entr->col_chain;
 470                 m->last_row_el[row]     = &entr->row_chain;
 471
 472                 /* add row_chain entry */
 473                 entr->row_chain.next = leftof->next;
 474                 leftof->next = &entr->row_chain;
 475
 476                 /* add col_chain entry */
 477                 entr->col_chain.next = above->next;
 478                 above->next = &entr->col_chain;
 479
 480                 m->entries++;
 481         }
 482 }
 483
 484 double matrix_get(const sp_matrix_t *m, int row, int col) {
 485         sp_matrix_list_head_t *dummy, *dummy2;
 486         matrix_elem_t *me;
 487
 488         if (is_empty_row(row) || is_empty_col(col))
 489                 return 0.0;
 490
 491         if (m->maxrow < m->maxcol)
 492                 me = _m_search_in_col(m, row, col, &dummy, &dummy2);
 493         else
 494                 me = _m_search_in_row(m, row, col, &dummy, &dummy2);
 495
 496         if (me)
 497                 assert(me->col == col && me->row == row);
 498
 499         return me ? me->val : 0.0;
 500 }
 501
 502 int matrix_get_entries(const sp_matrix_t *m) {
 503         return m->entries;
 504 }
 505
 506 int matrix_get_rowcount(const sp_matrix_t *m) {
 507         return m->maxrow + 1;
 508 }
 509
 510 int matrix_get_colcount(const sp_matrix_t *m) {
 511         return m->maxcol + 1;
 512 }
 513
 514 const matrix_elem_t *matrix_row_first(sp_matrix_t *m, int row) {
 515         if (is_empty_row(row))
 516                 return NULL;
 517
 518         m->dir   = right;
 519         m->first = m->rows[row];
 520         m->last  = m->first->next;
 521         m->next  = m->last ? m->last->next : NULL;
 522
 523         assert (list_entry_by_row(m->last)->row == row);
 524
 525         return list_entry_by_row(m->last);
 526 }
 527
 528 const matrix_elem_t *matrix_col_first(sp_matrix_t *m, int col) {
 529         if (is_empty_col(col))
 530                 return NULL;
 531
 532         m->dir   = down;
 533         m->first = m->cols[col];
 534         m->last  = m->first->next;
 535         m->next  = m->last ? m->last->next : NULL;
 536
 537         assert (list_entry_by_col(m->last)->col == col);
 538
 539         return list_entry_by_col(m->last);
 540 }
 541
 542 static inline const matrix_elem_t *matrix_first_from(sp_matrix_t *m, int startrow) {
 543         const matrix_elem_t *res;
 544         int i;
 545
 546         for (i = startrow; i <= m->maxrow; ++i) {
 547                 res = matrix_row_first(m, i);
 548                 if (res) {
 549                         m->iter_row = i;
 550                         m->dir      = all;
 551                         return res;
 552                 }
 553         }
 554
 555         return NULL;
 556 }
 557
 558 const matrix_elem_t *matrix_first(sp_matrix_t *m) {
 559         return matrix_first_from(m, 0);
 560 }
 561
 562 const matrix_elem_t *matrix_next(sp_matrix_t *m) {
 563         assert(m->first && "Start iteration with matrix_???_first, before calling me!");
 564
 565         if (m->next == NULL) {
 566                 if (m->dir == all)
 567                         return matrix_first_from(m, ++m->iter_row);
 568                 else
 569                         return NULL;
 570         }
 571
 572         m->last = m->next;
 573         m->next = m->next->next;
 574
 575         if (m->dir == down)
 576                 return list_entry_by_col(m->last);
 577         else /* right or all */
 578                 return list_entry_by_row(m->last);
 579 }
 580
 581 static int cmp_count(const void *e1, const void *e2) {
 582         return (int *)e2 - (int *)e1;
 583 }
 584
 585 static inline void matrix_fill_row(sp_matrix_t *m, int row, bitset_t *fullrow) {
 586         const matrix_elem_t *e;
 587         bitset_set(fullrow, row);
 588         matrix_foreach_in_col(m, row, e) {
 589                 if (! bitset_is_set(fullrow, e->row)) {
 590                         assert(0.0 == matrix_get(m, e->col, e->row));
 591                         matrix_set(m, e->col, e->row, e->val);
 592                         matrix_set(m, e->row, e->col, 0.0);
 593                 }
 594         }
 595 }
 596
 597 void matrix_optimize(sp_matrix_t *m) {
 598         int i, size, redo;
 599         int *c;
 600         const matrix_elem_t *e;
 601         bitset_t *fullrow;
 602
 603         size = MAX(m->maxcol, m->maxrow)+1;
 604
 605         /* kill all double-entries (Mij and Mji are set) */
 606         matrix_foreach(m, e) {
 607                 double t_val;
 608
 609                 assert(e->row != e->col && "Root has itself as arg. Ok. But the arg (=root) will always have the same color as root");
 610                 t_val = matrix_get(m, e->col, e->row);
 611                 if (fabs(t_val) > 1e-10) {
 612                         matrix_set(m, e->col, e->row, 0);
 613                         matrix_set(m, e->row, e->col, e->val + t_val);
 614                 }
 615         }
 616
 617         c       = alloca(size * sizeof(*c));
 618         redo    = 1;
 619         fullrow = bitset_alloca(size);
 620
 621         /* kill 'all' rows containing only 1 entry */
 622         while (redo) {
 623                 redo = 0;
 624                 /* count elements in rows */
 625                 memset(c, 0, size * sizeof(*c));
 626
 627                 matrix_foreach(m, e)
 628                         c[e->row]++;
 629
 630                 for (i = 0; i<size; ++i)
 631                         if (c[i] == 1 && ! bitset_is_set(fullrow, i)) {
 632                                 redo = 1;
 633                                 /* if the other row isn't empty move the e in there, else fill e's row */
 634                                 if (e = matrix_row_first(m, i), e) {
 635                                         if (c[e->col] > 0)
 636                                                 matrix_fill_row(m, e->col, fullrow);
 637                                         else
 638                                                 matrix_fill_row(m, e->row, fullrow);
 639                                 }
 640                         }
 641         }
 642
 643
 644         memset(c, 0, size * sizeof(*c));
 645         matrix_foreach(m, e)
 646                 c[e->row]++;
 647
 648         qsort(c, size, sizeof(*c), cmp_count);
 649
 650         for (i = 0; i < size; ++i) {
 651                 if (! bitset_is_set(fullrow, i))
 652                         matrix_fill_row(m, i, fullrow);
 653         }
 654 }
 655
 656 void matrix_dump(sp_matrix_t *m, FILE *out, int factor) {
 657         int i, o, last_idx;
 658         const matrix_elem_t *e;
 659
 660         for (i = 0; i <= m->maxrow; ++i) {
 661                 last_idx = -1;
 662                 matrix_foreach_in_row(m, i, e) {
 663                         for (o = last_idx + 1; o < e->col; ++o)
 664                                 fprintf(out, " %4.1f" , 0.0);
 665
 666                         fprintf(out, " %4.1f", factor * e->val);
 667                         last_idx = e->col;
 668                 }
 669
 670                 for (o = last_idx + 1; o <= m->maxcol; ++o)
 671                         fprintf(out, " %4.1f" , 0.0);
 672
 673                 fprintf(out, "\n");
 674         }
 675 }
 676
 677 void matrix_self_test(int d) {
 678         int i, o;
 679         const matrix_elem_t *e;
 680         sp_matrix_t *m = new_matrix(10, 10);
 681
 682         for (i = 0; i < d; ++i)
 683                 for (o = 0; o < d; ++o)
 684                         matrix_set(m, i, o, i*o);
 685
 686         for (i = 0; i < d; ++i)
 687                 for (o = 0; o<d; ++o)
 688                         assert(matrix_get(m, i, o) == i*o);
 689
 690         i = 1;
 691         matrix_foreach_in_row(m,1,e) {
 692                 assert(e->val == i);
 693                 i++;
 694         }
 695         assert(!matrix_next(m)); /*iter must finish */
 696
 697         i = d-1;
 698         matrix_foreach_in_col(m,d-1,e) {
 699                 assert(e->val == i);
 700                 i += d-1;
 701         }
 702         assert(!matrix_next(m));
 703         del_matrix(m);
 704         m = new_matrix(16,16);
 705         matrix_set(m, 1,1,9);
 706         matrix_set(m, 1,2,8);
 707         matrix_set(m, 1,3,7);
 708
 709         matrix_set(m, 1,3,6);
 710         matrix_set(m, 1,2,5);
 711         matrix_set(m, 1,1,4);
 712         i = 1;
 713         matrix_foreach_in_row(m, 1, e) {
 714                 assert(e->row == 1 && e->col == i && e->val == i+3);
 715                 i++;
 716         }
 717         assert(i == 4);
 718         del_matrix(m);
 719
 720         m = new_matrix(5,5);
 721         matrix_set(m, 1,1,1);
 722         matrix_set(m, 2,2,2);
 723         matrix_set(m, 3,3,3);
 724         matrix_set(m, 3,5,4);
 725         matrix_set(m, 4,4,5);
 726         matrix_set(m, 5,5,6);
 727         for (i=1, e = matrix_first(m); e; ++i, e=matrix_next(m))
 728                 assert(e->val == i);
 729         assert(i == 7);
 730         matrix_set(m, 1,1,0);
 731         assert(5 == matrix_get_entries(m));
 732         del_matrix(m);
 733 }