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

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