SOR iteration test

sor_iteration.c

/*
 * sor_iteration.c
 *
 *  Created on: Jul 7, 2014
 *      Author: du
 */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stddef.h>

#define real float

static real distance(size_t N, real *x, real *y) {
	int i;
	real s = 0.0;
	real d;
	for (i = 0; i < N; i++) {
		d = x[i] - y[i];
		s += d * d;
	}
	return sqrt(s / N);
}

static real sor_update(size_t N, size_t i, real *A, real *b, real *x, real w) {
	int j;
	real s = 0.0;
	for (j = 0; j < N; j++) {
		if (i == j) {
			continue;
		}
		s += A[i * N + j] * x[j];
	}
	return (1.0 - w) * x[i] + w * (b[i] - s) / A[i * N + i];
}

/* solve Ax = b by SOR iteration method
 * A[N][N], x[N], b[N]
 */
int sor_iteration(size_t N, real *A, real *b, real *initial_x, real *x, real w,
		int iter_max, real eps) {
	int iter;
	int i, j, k;
	int converged = 0;

	real old_x[N];
	real d = 1e4;

	for (i = 0; i < N; i++) {
		x[i] = initial_x[i];
	}

	for (iter = 0; iter < iter_max; iter++) {
		memcpy(old_x, x, sizeof(old_x));
		for (i = 0; i < N; i++) {
			x[i] = sor_update(N, i, A, b, x, w);
		}
		printf("iter = %d\n", iter);
		for (j = 0; j < N; j++) {
			printf("%f ", x[j]);
		}
		printf("distance = %f\n", distance(N, x, old_x));

		if (distance(N, x, old_x) < eps) {
			converged = 1;
			break;
		}
	}

	return converged;
}

#define n 2
int main(int argc, char **argv) {
	real A[n * n] = { 2, 1, 3, 4 };
	real b[n] = { 10, 7 };
	real x[n] = { 0, 0 };

//	A[0][0] = 3.0; A[0][1] = 2.0;
//	A[1][0] = 7.0; A[1][1] = 5.0;

	int conv = sor_iteration(n, A, b, b, x, 1.0, 100, 1e-5);
	int i;
	printf("conv = %d\n", conv);
	for (i = 0; i < n; i++) {
		printf("x[%d] = %f\n", i, x[i]);
	}
}