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