#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define I 4

typedef double real;

real* solveLS_GS(real **a, real *b, int N) {
real *x0 = (real*)calloc(N, sizeof(real));
real *x = (real*)calloc(N, sizeof(real));
int i, j, iter;
real sum;

    for (iter = 0; iter < 32; iter++) { // adjust number of interations accordingly
        for (i = 0; i < N; i++) {
            sum = 0.;
            for (j = 0; j < N; j++) {
                if (j < i) sum += a[i][j] * x[j];
                else if (j > i) sum += a[i][j] * x0[j];
                else {} // j == i
            } // for_j
            x[i] = 1. / a[i][i] * (b[i] - sum);
        } // for_i
        printf("%d\t", iter);
        for (i = 0; i < N; i++) { x0[i] = x[i]; printf("%.4lf\t", x[i]); }
        printf("\n");
    } // for_iter

    free(x0);

    return x;
} // solveLS_GS

real* solveLS_J(real **a, real *b, int N) {
real *x0 = (real*)calloc(N, sizeof(real));
real *x = (real*)calloc(N, sizeof(real));
int i, j, iter;
real sum;

    for (iter = 0; iter < 60; iter++) { // adjust number of interations accordingly
        for (i = 0; i < N; i++) {
            sum = 0.;
            for (j = 0; j < N; j++) if (j != i) sum += a[i][j] * x0[j];
            x[i] = 1. / a[i][i] * (b[i] - sum);
        } // for_i
        printf("%d\t", iter);
        for (i = 0; i < N; i++) { x0[i] = x[i]; printf("%.4lf\t", x[i]); }
        printf("\n");
    } // for_iter

    free(x0);

    return x;
} // solveLS_J

int main() {
FILE *pF = NULL;
int i;

	real **a = (real**)calloc(I, sizeof(real*));
	for (int i = 0; i < I; i++) a[i] = (real*)calloc(I, sizeof(real));
	memcpy(a, (real *[]) {
						(real[]){-2.,  1.,  0.,  0. },
						(real[]){ 1., -2.,  1.,  0. },
						(real[]){ 0.,  1., -2.,  1. },
						(real[]){ 0.,  0.,  1., -2. }
						}, I * I * sizeof(real));

	real *b = (real*)calloc(I, sizeof(real));
	memcpy(b, (real[I]) { 5., 1., 0., 8. }, I * sizeof(real));

	real *x;
#ifndef GAUSS // switch between #ifdef and #ifndef
    x = solveLS_GS(a, b, I);
#else
    x = solveLS_J(a, b, I);
#endif // GAUSS

    for (i = 0; i < I; i++) free(a[i]);
	free(a);
    free(b);
	free(x);

    return 0;
} // main