Commit 0553e2d5 authored by Matthias Braun's avatar Matthias Braun
Browse files

gaussseidel: Cleanup

parent d3b482aa
#ifndef MATRIX_H_
#define MATRIX_H_
#ifndef FIRM_ADT_GAUSSSEIDEL_H
#define FIRM_ADT_GAUSSSEIDEL_H
#include <stdio.h>
......@@ -12,36 +12,37 @@ typedef struct gs_matrix_t gs_matrix_t;
* row_init X col_init. Matrix cannot grow beyond these init values.
* All elements are initially (implicitly) set to 0.
*/
gs_matrix_t *gs_new_matrix(int n_init_rows, int n_init_cols);
FIRM_API gs_matrix_t *gs_new_matrix(unsigned n_init_rows, unsigned n_init_cols);
/**
* Free space used by matrix m
*/
void gs_delete_matrix(gs_matrix_t *m);
FIRM_API void gs_delete_matrix(gs_matrix_t *m);
/**
* Sets m[row, col] to val
*/
void gs_matrix_set(gs_matrix_t *m, unsigned row, unsigned col, double val);
FIRM_API void gs_matrix_set(gs_matrix_t *m, unsigned row, unsigned col,
double val);
/**
* Returns the value stored in m[row, col].
*/
double gs_matrix_get(const gs_matrix_t *m, unsigned row, unsigned col);
FIRM_API double gs_matrix_get(const gs_matrix_t *m, unsigned row, unsigned col);
/**
* Performs one step of the Gauss-Seidel algorithm
* @p m The iteration matrix
* @p x The iteration vector
*/
double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x);
FIRM_API double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x);
unsigned gs_matrix_get_n_entries(const gs_matrix_t *m);
FIRM_API unsigned gs_matrix_get_n_entries(const gs_matrix_t *m);
/**
* Dumps the matrix factor*m to the stream @p out.
*/
void gs_matrix_dump(const gs_matrix_t *m, FILE *out);
FIRM_API void gs_matrix_dump(const gs_matrix_t *m, FILE *out);
#include "../end.h"
......
#include "gaussseidel.h"
#include <assert.h>
#include <math.h>
#include <string.h>
#include "xmalloc.h"
#include "gaussseidel.h"
#include "util.h"
typedef struct col_val_t {
typedef struct {
double v;
unsigned col_idx;
} col_val_t;
typedef struct row_col_t {
typedef struct {
unsigned c_cols;
unsigned n_cols;
double diag;
......@@ -48,7 +49,7 @@ static void alloc_rows(gs_matrix_t *m, unsigned c_rows, unsigned c_cols,
}
}
gs_matrix_t *gs_new_matrix(int n_init_rows, int n_init_cols)
gs_matrix_t *gs_new_matrix(unsigned n_init_rows, unsigned n_init_cols)
{
gs_matrix_t *res = XMALLOCZ(gs_matrix_t);
alloc_rows(res, n_init_rows, n_init_cols, 0);
......@@ -91,7 +92,7 @@ void gs_matrix_set(gs_matrix_t *m, unsigned row, unsigned col, double val)
return;
}
// Search for correct column
/* Search for correct column */
col_val_t *cols = the_row->cols;
unsigned min = 0;
unsigned max = the_row->n_cols;
......@@ -104,10 +105,10 @@ void gs_matrix_set(gs_matrix_t *m, unsigned row, unsigned col, double val)
max = MIN(c, max-1);
else
break;
c = (max+min)/2;
c = (min+max)/2;
}
// Have we found the entry?
/* Have we found the entry? */
if (c < the_row->n_cols && the_row->cols[c].col_idx == col) {
the_row->cols[c].v = val;
if (val == 0.0)
......@@ -115,23 +116,24 @@ void gs_matrix_set(gs_matrix_t *m, unsigned row, unsigned col, double val)
return;
}
// We haven't found the entry, so we must create a new one.
// Is there enough space?
/* We haven't found the entry, so we must create a new one.
* Is there enough space? */
if (the_row->n_cols <= the_row->c_cols)
alloc_cols(the_row, the_row->c_cols + 16);
// Shift right-most entries to the right by one
/* Shift right-most entries to the right by one */
for (unsigned i = the_row->n_cols; i > c; --i)
the_row->cols[i] = the_row->cols[i-1];
// Finally insert the new entry
/* Finally insert the new entry */
the_row->n_cols++;
the_row->cols[c].col_idx = col;
the_row->cols[c].v = val;
// Check that the entries are sorted
/* Check that the entries are sorted */
assert(c==0 || the_row->cols[c-1].col_idx < the_row->cols[c].col_idx);
assert(c>=the_row->n_cols-1 || the_row->cols[c].col_idx < the_row->cols[c+1].col_idx);
assert(c>=the_row->n_cols-1
|| the_row->cols[c].col_idx < the_row->cols[c+1].col_idx);
}
double gs_matrix_get(const gs_matrix_t *m, unsigned row, unsigned col)
......@@ -141,7 +143,7 @@ double gs_matrix_get(const gs_matrix_t *m, unsigned row, unsigned col)
if (row == col)
return the_row->diag != 0.0 ? 1.0 / the_row->diag : 0.0;
// Search for correct column
/* Search for correct column */
unsigned c;
for (c = 0; c < the_row->n_cols && the_row->cols[c].col_idx < col; ++c) {
}
......@@ -177,7 +179,7 @@ double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x)
double old = x[r];
double nw = - sum * row->diag;
// nw = old - overdrive * (old + sum * row->diag);
/* nw = old - overdrive * (old + sum * row->diag); */
res += fabs(old - nw);
x[r] = nw;
}
......@@ -187,10 +189,10 @@ double gs_matrix_gauss_seidel(const gs_matrix_t *m, double *x)
void gs_matrix_dump(const gs_matrix_t *m, FILE *out)
{
unsigned size = m->c_rows;
double *elems = XMALLOCN(double, size);
unsigned size = m->c_rows;
double *elems = XMALLOCN(double, size);
// The rows which have some content
/* The rows which have some content */
for (unsigned r = 0; r < size; ++r) {
row_col_t *row = &m->rows[r];
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment