// sor.cpp : Defines the entry point for the console application.
// SOR - Successive Over-Relaxation method for solving linear equations
//

#include "stdafx.h"
#define   N   4   //  For this example we pick a relatively small value

void solve(double a[N][N], double b[N], double x[N], int n)
{   //  This function performs one iteration of Gauss-Seidel method
	int i, j, k;

	for (i=0; i<n; i++) {
		x[i] = b[i];
		for (k=0; k<n; k++) {
			if (k!=i) x[i] = x[i] - a[i][k] * x[k];
		}
		x[i] = x[i] / a[i][i];
	}
	return;
}
int _tmain(int argc, _TCHAR* argv[])
{
	double alpha = 1.2;                //  The RELAXATION PARAMETER (usually between 1 and 2)
	double a[N][N], b[N], x[N], xgs[N];//  Allocate storage for matrix a, vectors b, x, and xgs
	int i, j, k, n, kk;
            // The matrix a is stored in an array of size N x N
	n = N;  // n must be less than or equal to N
            // n x n  is the size of the coefficient matrix for the linear system
	        // In other words, this program could be used to solve any system of
	        // size from  1 x 1  up to  N x N
	a[0][0] = -2.0;    a[0][1] = 1.0;     a[0][2] = 0.0;     a[0][3] = 0.0;    
	a[1][0] = 1.0;     a[1][1] = -2.0;    a[1][2] = 1.0;     a[1][3] = 0.0; 
	a[2][0] = 0.0;     a[2][1] = 1.0;     a[2][2] = -2.0;    a[2][3] = 1.0; 
	a[3][0] = 0.0;     a[3][1] = 0.0;     a[3][2] = 1.0;     a[3][3] = -2.0; 
	b[0] = 1.0;
	b[1] = 2.0;	
	b[2] = 3.0;	
	b[3] = 4.0;	
	xgs[0] = x[0] = 1.0;           //  Guess the solution vectors
	xgs[1] = x[1] = 1.0;
	xgs[2] = x[2] = 1.0;
	xgs[3] = x[3] = 1.0;
	printf(" The original  a  matrix and  b  vector are:\n");
	for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
			printf(" %10.5lf", a[i][j]);
		}
		printf("   %10.5lf\n", b[i]);
	}
	printf("------------------------------------------\n");
	for (kk=0; kk<40; kk++) {
        solve(a, b, xgs, n);
		for (i=0; i<n; i++) {
			x[i] = x[i] + alpha * (xgs[i]-x[i]);
			xgs[i] = x[i];
		}
	    printf("------------------------------------------\n");
	    printf(" After %d iterations, the solution  x  is:\n", kk+1);
	    for (i=0; i<n; i++) {
		    printf(" %10.5lf", x[i]);
	    }
	    printf("\n");
	}
	return 0;
}