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

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