/*
 * gauss-seidel.c --- Gauss-Seidel 法
 *
 *   コンパイル
 *     gcc -c matutil.c
 *     gcc -o gauss-seidel gauss-seidel.c matutil.o -lm
 */

#include <stdio.h>
#include <math.h>
#include "matutil.h"

double dist(int n, vector x, vector y)
{
  int i;
  double s;
  s = 0;
  for (i = 0; i < n; i++) s += (x[i] - y[i]) * (x[i] - y[i]);
  return sqrt(s);
}

main()
{
  int i, j, k, n;
  double s;
  matrix a;
  vector exact_x, x, b;

  n = 3;
  a = new_matrix(n, n);
  x = new_vector(n);
  exact_x = new_vector(n);
  b = new_vector(n);

  a[0][0] = 2; a[0][1] = 1; a[0][2] = 0;
  a[1][0] = 1; a[1][1] = 4; a[1][2] = 2;
  a[2][0] = 2; a[2][1] = 3; a[2][2] = 8;

  x[0] = 1; x[1] = 2; x[2] = 3;
  for (i = 0; i < n; i++)
    exact_x[i] = x[i];

  printf("A=\n");
  for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++)
      printf("%g ", a[i][j]);
    printf("\n");
  }

  printf("b=\n");
  for (i = 0; i < n; i++)
    b[i] = dotprod(n, a[i], x);
  for (i = 0; i < n; i++)
    printf("%g ", b[i]);
  printf("\n");

  for (i = 0; i < n; i++)
    x[i] = 0;

  printf("反復開始\n");
  printf("        x1                    x2               x3           誤差\n");

  for (i = 0; i < n; i++)
    printf("%18.15f ", x[i]);
  printf(" %e\n", dist(n, x, exact_x));

  for (k = 0; k < 30; k++) {
    for (i = 0; i < n; i++) {
      s = b[i];
      for (j = 0; j < n; j++)
        if (j != i)
          s -= a[i][j] * x[j];
      x[i] = s / a[i][i];
    }
    for (i = 0; i < n; i++)
      printf("%18.15f ", x[i]);
    printf(" %e\n", dist(n, x, exact_x));
  }
}