nsz Git - libfirm/blob - ir/adt/gaussseidel.c

   1 #ifdef HAVE_CONFIG_H
   2 # include "config.h"
   3 #endif
   4
   5 #include <assert.h>
   6 #include <math.h>
   7 #include <string.h>
   8 #include "xmalloc.h"
   9 #include "gaussseidel.h"
  10 #include "firm_config.h"
  11 #include "util.h"
  12
  13 #define MAX(x,y)   ((x) > (y) ? (x) : (y))
  14 #define MIN(x,y)   ((x) < (y) ? (x) : (y))
  15
  16 /**
  17  * The number of newly allocated rows (realloc)
  18  * when there is no more room. Must be >= 1.
  19  */
  20 #define ROW_INCREASE_FACTOR 1.2
  21
  22 /**
  23  * The number of newly allocated cols (realloc)
  24  * when there is no more room. Must be >= 1.
  25  */
  26 #define COL_INCREASE 2
  27
  28 typedef struct _col_val_t {
  29         double v;
  30         int col_idx;
  31 } col_val_t;
  32
  33 typedef struct _row_col_t {
  34         int c_cols;
  35         int n_cols;
  36         double diag;
  37         col_val_t *cols;
  38 } row_col_t;
  39
  40 struct _gs_matrix_t {
  41         int initial_col_increase;
  42         int c_rows;
  43         int n_zero_entries;           ///< Upper bound on number of entries equal to 0.0
  44         row_col_t *rows;
  45 };
  46
  47 static INLINE void alloc_cols(row_col_t *row, int c_cols) {
  48         assert(c_cols > row->c_cols);
  49         row->c_cols = c_cols;
  50         row->cols   = XREALLOC(row->cols, col_val_t, c_cols);
  51 }
  52
  53 static INLINE void alloc_rows(gs_matrix_t *m, int c_rows, int c_cols, int begin_init) {
  54         int i;
  55         assert(c_rows > m->c_rows);
  56
  57         m->c_rows = c_rows;
  58         m->rows   = XREALLOC(m->rows, row_col_t, c_rows);
  59
  60         for (i = begin_init; i < c_rows; ++i) {
  61                 m->rows[i].c_cols = 0;
  62                 m->rows[i].n_cols = 0;
  63                 m->rows[i].diag   = 0.0;
  64                 m->rows[i].cols   = NULL;
  65                 if (c_cols > 0)
  66                         alloc_cols(&m->rows[i], c_cols);
  67         }
  68 }
  69
  70 gs_matrix_t *gs_new_matrix(int n_init_rows, int n_init_cols) {
  71         gs_matrix_t *res = XMALLOCZ(gs_matrix_t);
  72         if (n_init_rows < 16)
  73                 n_init_rows = 16;
  74         res->initial_col_increase = n_init_cols;
  75         alloc_rows(res, n_init_rows, n_init_cols, 0);
  76         return res;
  77 }
  78
  79 void gs_delete_matrix(gs_matrix_t *m) {
  80         int i;
  81         for (i = 0; i < m->c_rows; ++i) {
  82                 if (m->rows[i].c_cols)
  83                         xfree(m->rows[i].cols);
  84         }
  85         if (m->c_rows)
  86                 xfree(m->rows);
  87         xfree(m);
  88 }
  89
  90 unsigned gs_matrix_get_n_entries(const gs_matrix_t *m) {
  91         int i;
  92         unsigned n_entries = 0;
  93
  94         for (i = 0; i < m->c_rows; ++i) {
  95                 n_entries += m->rows[i].n_cols;
  96                 n_entries += (m->rows[i].diag != 0.0) ? 1 : 0;
  97         }
  98
  99         return n_entries - m->n_zero_entries;
 100 }
 101
 102 int gs_matrix_get_sizeof_allocated_memory(const gs_matrix_t *m) {
 103         int i, n_col_val_ts = 0;
 104         for (i = 0; i < m->c_rows; ++i)
 105                 n_col_val_ts += m->rows[i].c_cols;
 106
 107         return n_col_val_ts * sizeof(col_val_t) + m->c_rows * sizeof(row_col_t) + sizeof(gs_matrix_t);
 108 }
 109
 110 void gs_matrix_assure_row_capacity(gs_matrix_t *m, int row, int min_capacity) {
 111         row_col_t *the_row = &m->rows[row];
 112         if (the_row->c_cols < min_capacity)
 113                 alloc_cols(the_row, min_capacity);
 114 }
 115
 116 void gs_matrix_trim_row_capacities(gs_matrix_t *m) {
 117         int i;
 118         for (i = 0; i < m->c_rows; ++i) {
 119                 row_col_t *the_row = &m->rows[i];
 120                 if (the_row->c_cols) {
 121                         the_row->c_cols    = the_row->n_cols;
 122                         if (the_row->c_cols)
 123                                 the_row->cols = XREALLOC(the_row->cols, col_val_t, the_row->c_cols);
 124                         else
 125                                 xfree(the_row->cols);
 126                 }
 127         }
 128 }
 129
 130 void gs_matrix_delete_zero_entries(gs_matrix_t *m) {
 131         int i, read_pos;
 132         for (i = 0; i < m->c_rows; ++i) {
 133                 row_col_t *the_row = &m->rows[i];
 134                 int write_pos = 0;
 135
 136                 for (read_pos = 0; read_pos < the_row->n_cols; ++read_pos)
 137                         if (the_row->cols[read_pos].v != 0.0 && read_pos != write_pos)
 138                                 the_row->cols[write_pos++] = the_row->cols[read_pos];
 139
 140                 the_row->n_cols = write_pos;
 141         }
 142         m->n_zero_entries = 0;
 143 }
 144
 145 void gs_matrix_set(gs_matrix_t *m, int row, int col, double val) {
 146         row_col_t *the_row;
 147         col_val_t *cols;
 148         int min, max, c, i;
 149
 150         if (row >= m->c_rows) {
 151                 int new_c_rows = (int)(ROW_INCREASE_FACTOR * row);
 152                 alloc_rows(m, new_c_rows, m->initial_col_increase, m->c_rows);
 153         }
 154
 155         the_row = &m->rows[row];
 156
 157         if (row == col) {
 158                 /* Note that we store the diagonal inverted to turn divisions to mults in
 159                  * matrix_gauss_seidel(). */
 160                 assert(val != 0.0);
 161                 the_row->diag = 1.0 / val;
 162                 return;
 163         }
 164
 165         // Search for correct column
 166         cols = the_row->cols;
 167         min  = 0;
 168         max  = the_row->n_cols;
 169         c    = (max+min)/2;
 170         while (min < max) {
 171                 int idx = cols[c].col_idx;
 172                 if (idx < col)
 173                         min = MAX(c, min+1);
 174                 else if (idx > col)
 175                         max = MIN(c, max-1);
 176                 else
 177                         break;
 178                 c = (max+min)/2;
 179         }
 180
 181         // Have we found the entry?
 182         if (c < the_row->n_cols && the_row->cols[c].col_idx == col) {
 183                 the_row->cols[c].v = val;
 184                 if (val == 0.0)
 185                         m->n_zero_entries++;
 186                 return;
 187         }
 188
 189         // We haven't found the entry, so we must create a new one.
 190         // Is there enough space?
 191         if (the_row->c_cols == the_row->n_cols)
 192                 alloc_cols(the_row, the_row->c_cols + COL_INCREASE);
 193
 194         // Shift right-most entries to the right by one
 195         for (i = the_row->n_cols; i > c; --i)
 196                 the_row->cols[i] = the_row->cols[i-1];
 197
 198         // Finally insert the new entry
 199         the_row->n_cols++;
 200         the_row->cols[c].col_idx = col;
 201         the_row->cols[c].v = val;
 202
 203         // Check that the entries are sorted
 204         assert(c==0 || the_row->cols[c-1].col_idx < the_row->cols[c].col_idx);
 205         assert(c>=the_row->n_cols-1 || the_row->cols[c].col_idx < the_row->cols[c+1].col_idx);
 206 }
 207
 208 double gs_matrix_get(const gs_matrix_t *m, int row, int col) {
 209         row_col_t *the_row;
 210         int c;
 211
 212         if (row >= m->c_rows)
 213                 return 0.0;
 214
 215         the_row = &m->rows[row];
 216
 217         if (row == col)
 218                 return the_row->diag != 0.0 ? 1.0 / the_row->diag : 0.0;
 219
 220         // Search for correct column
 221         for (c = 0; c < the_row->n_cols && the_row->cols[c].col_idx < col; ++c);
 222
 223         if (c >= the_row->n_cols || the_row->cols[c].col_idx > col)
 224                 return 0.0;
 225
 226         assert(the_row->cols[c].col_idx == col);
 227         return the_row->cols[c].v;
 228 }
 229
 230 /* NOTE: You can slice out miss_rate and weights.
 231  * This does ONE step of gauss_seidel. Termination must be checked outside!
 232  * This solves m*x=0. You must add stuff for m*x=b. See wikipedia german and english article. Should be simple.
 233  * param a is the number of rows in the matrix that should be considered.
 234  *
 235  * Note that the diagonal element is stored separately in this matrix implementation.
 236  * */
 237 double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x, int n) {
 238         double res = 0.0;
 239         int r;
 240
 241         assert(n <= m->c_rows);
 242
 243         for (r = 0; r < n; ++r) {
 244                 row_col_t *row  = &m->rows[r];
 245                 col_val_t *cols = row->cols;
 246                 double sum, old, nw;
 247                 int c;
 248
 249                 sum = 0.0;
 250                 for (c = 0; c < row->n_cols; ++c) {
 251                         int col_idx = cols[c].col_idx;
 252                         sum += cols[c].v * x[col_idx];
 253                 }
 254
 255                 old  = x[r];
 256                 nw   = - sum * row->diag;
 257                 // nw   = old - overdrive * (old + sum * row->diag);
 258                 res += fabs(old - nw);
 259                 x[r] = nw;
 260         }
 261
 262         return res;
 263 }
 264
 265 void gs_matrix_export(const gs_matrix_t *m, double *nw, int size)
 266 {
 267         int effective_rows = MIN(size, m->c_rows);
 268         int c, r;
 269
 270         memset(nw, 0, size * size * sizeof(*nw));
 271         for (r=0; r < effective_rows; ++r) {
 272                 row_col_t *row = &m->rows[r];
 273                 int base       = r * size;
 274
 275                 assert(row->diag != 0.0);
 276                 nw[base + r] = 1.0 / row->diag;
 277                 for (c = 0; c < row->n_cols; ++c) {
 278                         int col_idx = row->cols[c].col_idx;
 279                         nw[base + col_idx] = row->cols[c].v;
 280                 }
 281         }
 282 }
 283
 284 void gs_matrix_dump(const gs_matrix_t *m, int a, int b, FILE *out) {
 285         int effective_rows = MIN(a, m->c_rows);
 286         int r, c, i;
 287         double *elems = XMALLOCN(double, b);
 288
 289         // The rows which have some content
 290         for (r=0; r < effective_rows; ++r) {
 291                 row_col_t *row = &m->rows[r];
 292
 293                 memset(elems, 0, b * sizeof(*elems));
 294
 295                 for (c = 0; c < row->n_cols; ++c) {
 296                         int col_idx = row->cols[c].col_idx;
 297                         elems[col_idx] = row->cols[c].v;
 298                 }
 299                 elems[r] = row->diag != 0.0 ? 1.0 / row->diag : 0.0;
 300
 301                 for (i = 0; i < b; ++i)
 302                         if (elems[i] != 0.0)
 303                                 fprintf(out, "%+4.4f ", elems[i]);
 304                         else
 305                                 fprintf(out, "        ");
 306                 fprintf(out, "\n");
 307         }
 308
 309         // Append 0-rows to fit height of matrix
 310         for (r=effective_rows; r < a; ++r) {
 311                 for (c=0; c < b; ++c)
 312                                 fprintf(out, "        ");
 313                 fprintf(out, "\n");
 314         }
 315
 316         xfree(elems);
 317 }
 318
 319 void gs_matrix_self_test(int d) {
 320         int i, o;
 321         gs_matrix_t *m = gs_new_matrix(10, 10);
 322
 323         for (i=0; i<d; ++i)
 324                 for (o=0; o<d; ++o)
 325                         gs_matrix_set(m, i, o, i*o);
 326
 327         for (i=0; i<d; ++i)
 328                 for (o=0; o<d; ++o)
 329                         assert(gs_matrix_get(m, i, o) == i*o);
 330         gs_delete_matrix(m);
 331 }