5 #include "gaussseidel.h"
6 #include "firm_config.h"
9 #define MAX(x,y) ((x) > (y) ? (x) : (y))
10 #define MIN(x,y) ((x) < (y) ? (x) : (y))
13 * The number of newly allocated rows (realloc)
14 * when there is no more room. Must be >= 1.
16 #define ROW_INCREASE_FACTOR 1.2
19 * The number of newly allocated cols (realloc)
20 * when there is no more room. Must be >= 1.
22 #define COL_INCREASE 2
24 typedef struct _col_val_t {
29 typedef struct _row_col_t {
37 int initial_col_increase;
39 int n_zero_entries; ///< Upper bound on number of entries equal to 0.0
43 static INLINE void alloc_cols(row_col_t *row, int c_cols) {
44 assert(c_cols > row->c_cols);
46 row->cols = xrealloc(row->cols, c_cols * sizeof(*row->cols));
49 static INLINE void alloc_rows(gs_matrix_t *m, int c_rows, int c_cols, int begin_init) {
51 assert(c_rows > m->c_rows);
54 m->rows = xrealloc(m->rows, c_rows * sizeof(*m->rows));
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;
62 alloc_cols(&m->rows[i], c_cols);
66 gs_matrix_t *gs_new_matrix(int n_init_rows, int n_init_cols) {
67 gs_matrix_t *res = xmalloc(sizeof(*res));
68 memset(res, 0, sizeof(*res));
71 res->initial_col_increase = n_init_cols;
72 alloc_rows(res, n_init_rows, n_init_cols, 0);
76 void gs_delete_matrix(gs_matrix_t *m) {
78 for (i = 0; i < m->c_rows; ++i) {
79 if (m->rows[i].c_cols)
80 xfree(m->rows[i].cols);
87 unsigned gs_matrix_get_n_entries(const gs_matrix_t *m) {
89 unsigned n_entries = 0;
91 for (i = 0; i < m->c_rows; ++i) {
92 n_entries += m->rows[i].n_cols;
93 n_entries += (m->rows[i].diag != 0.0) ? 1 : 0;
96 return n_entries - m->n_zero_entries;
99 int gs_matrix_get_sizeof_allocated_memory(const gs_matrix_t *m) {
100 int i, n_col_val_ts = 0;
101 for (i = 0; i < m->c_rows; ++i)
102 n_col_val_ts += m->rows[i].c_cols;
104 return n_col_val_ts * sizeof(col_val_t) + m->c_rows * sizeof(row_col_t) + sizeof(gs_matrix_t);
107 void gs_matrix_assure_row_capacity(gs_matrix_t *m, int row, int min_capacity) {
108 row_col_t *the_row = &m->rows[row];
109 if (the_row->c_cols < min_capacity)
110 alloc_cols(the_row, min_capacity);
113 void gs_matrix_trim_row_capacities(gs_matrix_t *m) {
115 for (i = 0; i < m->c_rows; ++i) {
116 row_col_t *the_row = &m->rows[i];
117 if (the_row->c_cols) {
118 the_row->c_cols = the_row->n_cols;
120 the_row->cols = xrealloc(the_row->cols, the_row->c_cols * sizeof(*the_row->cols));
122 xfree(the_row->cols);
127 void gs_matrix_delete_zero_entries(gs_matrix_t *m) {
129 for (i = 0; i < m->c_rows; ++i) {
130 row_col_t *the_row = &m->rows[i];
133 for (read_pos = 0; read_pos < the_row->n_cols; ++read_pos)
134 if (the_row->cols[read_pos].v != 0.0 && read_pos != write_pos)
135 the_row->cols[write_pos++] = the_row->cols[read_pos];
137 the_row->n_cols = write_pos;
139 m->n_zero_entries = 0;
142 void gs_matrix_set(gs_matrix_t *m, int row, int col, double val) {
147 if (row >= m->c_rows) {
148 int new_c_rows = (int)(ROW_INCREASE_FACTOR * row);
149 alloc_rows(m, new_c_rows, m->initial_col_increase, m->c_rows);
152 the_row = &m->rows[row];
155 /* Note that we store the diagonal inverted to turn divisions to mults in
156 * matrix_gauss_seidel(). */
158 the_row->diag = 1.0 / val;
162 // Search for correct column
163 cols = the_row->cols;
165 max = the_row->n_cols;
168 int idx = cols[c].col_idx;
178 // Have we found the entry?
179 if (c < the_row->n_cols && the_row->cols[c].col_idx == col) {
180 the_row->cols[c].v = val;
186 // We haven't found the entry, so we must create a new one.
187 // Is there enough space?
188 if (the_row->c_cols == the_row->n_cols)
189 alloc_cols(the_row, the_row->c_cols + COL_INCREASE);
191 // Shift right-most entries to the right by one
192 for (i = the_row->n_cols; i > c; --i)
193 the_row->cols[i] = the_row->cols[i-1];
195 // Finally insert the new entry
197 the_row->cols[c].col_idx = col;
198 the_row->cols[c].v = val;
200 // Check that the entries are sorted
201 assert(c==0 || the_row->cols[c-1].col_idx < the_row->cols[c].col_idx);
202 assert(c>=the_row->n_cols-1 || the_row->cols[c].col_idx < the_row->cols[c+1].col_idx);
205 double gs_matrix_get(const gs_matrix_t *m, int row, int col) {
209 if (row >= m->c_rows)
212 the_row = &m->rows[row];
215 return the_row->diag != 0.0 ? 1.0 / the_row->diag : 0.0;
217 // Search for correct column
218 for (c = 0; c < the_row->n_cols && the_row->cols[c].col_idx < col; ++c);
220 if (c >= the_row->n_cols || the_row->cols[c].col_idx > col)
223 assert(the_row->cols[c].col_idx == col);
224 return the_row->cols[c].v;
227 /* NOTE: You can slice out miss_rate and weights.
228 * This does ONE step of gauss_seidel. Termination must be checked outside!
229 * This solves m*x=0. You must add stuff for m*x=b. See wikipedia german and english article. Should be simple.
230 * param a is the number of rows in the matrix that should be considered.
232 * Note that the diagonal element is stored separately in this matrix implementation.
234 double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x, int n) {
238 assert(n <= m->c_rows);
240 for (r = 0; r < n; ++r) {
241 row_col_t *row = &m->rows[r];
242 col_val_t *cols = row->cols;
247 for (c = 0; c < row->n_cols; ++c) {
248 int col_idx = cols[c].col_idx;
249 sum += cols[c].v * x[col_idx];
253 nw = - sum * row->diag;
254 // nw = old - overdrive * (old + sum * row->diag);
255 res += fabs(old - nw);
262 void gs_matrix_export(const gs_matrix_t *m, double *nw, int size)
264 int effective_rows = MIN(size, m->c_rows);
267 memset(nw, 0, size * size * sizeof(*nw));
268 for (r=0; r < effective_rows; ++r) {
269 row_col_t *row = &m->rows[r];
272 assert(row->diag != 0.0);
273 nw[base + r] = 1.0 / row->diag;
274 for (c = 0; c < row->n_cols; ++c) {
275 int col_idx = row->cols[c].col_idx;
276 nw[base + col_idx] = row->cols[c].v;
281 void gs_matrix_dump(const gs_matrix_t *m, int a, int b, FILE *out) {
282 int effective_rows = MIN(a, m->c_rows);
284 double *elems = xmalloc(b * sizeof(*elems));
286 // The rows which have some content
287 for (r=0; r < effective_rows; ++r) {
288 row_col_t *row = &m->rows[r];
290 memset(elems, 0, b * sizeof(*elems));
292 for (c = 0; c < row->n_cols; ++c) {
293 int col_idx = row->cols[c].col_idx;
294 elems[col_idx] = row->cols[c].v;
296 elems[r] = row->diag != 0.0 ? 1.0 / row->diag : 0.0;
298 for (i = 0; i < b; ++i)
300 fprintf(out, "%+4.4f ", elems[i]);
306 // Append 0-rows to fit height of matrix
307 for (r=effective_rows; r < a; ++r) {
308 for (c=0; c < b; ++c)
316 void gs_matrix_self_test(int d) {
318 gs_matrix_t *m = gs_new_matrix(10, 10);
322 gs_matrix_set(m, i, o, i*o);
326 assert(gs_matrix_get(m, i, o) == i*o);