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

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