#include <stdio.h>
#include <math.h>
/* seidel.c : 1999 3 13 by Oguni */
#define N 5
main(){
  double A[N][N]={5.,-1.,0.,0.,0.,-1.,5.,-1.,0.,0.,0.,-1.,5.,-1.,0.,0.,0.,-1.,5.,-1.,0.,0.,0.,-1.,5.};
  double b[N]={4.,3.,3.,3.,4.}, x[N]={1.,0.,0.,0.,0.};
  double tol[2]={0.000001, 0.000001};
  int n=5, itr[2]={20,1}, i;
  seidel(A,b,x,tol,n,itr);
  printf("反復回数 %d 誤差%f\n", itr[1], tol[1]);
  for (i=0; i<n; i++){ printf("x(%d)=%f\n", i, x[i]); }
}
/* function seidel */
seidel(A,b,x,tol,n,itr)
double A[N][N], b[N], x[N], tol[2];
int n, itr[2];
{
  double L[N][N], U[N][N], D[N][N], C[N][N];
  double xnew[N], w[N], s1, s2, ep;
  int i,j,k;
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      L[i][j]=0.; U[i][j]=0.; D[i][j]=0.; }}
  for (i=1; i<n; i++){
    for (j=0; j<=i-1; j++){
      L[i][j]=A[i][j]; }
  }
  for (i=0; i<n-1; i++){
    for (j=i+1; j<n; j++){
      U[i][j]=A[i][j]; }
  }
  for (i=0; i<n; i++){ D[i][i]=A[i][i]; }
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      C[i][j]=L[i][j]+D[i][j]; }
  }
  k=0;
  while (k<=itr[0]){
    ++k; mvsub(U,x,w,n);
    for (i=0; i<n; i++){
      w[i]=b[i]-w[i]; }
    forward(C,w,xnew,n);
    s1=0.0; s2=0.0;
    for (i=0; i<n; i++){
      s1+=(xnew[i]-x[i])*(xnew[i]-x[i]); s2+=x[i]*x[i]; }
    ep=sqrt(fabs(s1/s2));
    if (ep<tol[0]){
     itr[1]=k; tol[1]=ep; return;}
    for (i=0; i<n; i++){
      x[i]=xnew[i]; }
  }
  itr[1]=k; tol[1]=ep;
  return;
}
/* function */
forward(L,b,x,n)
double L[N][N], b[N], x[N];
int n;
{
  double tol1=0.0000001, s;
  int i,j;
  x[0]=b[0]/L[0][0];
  for (i=1; i<n; i++){
    if (fabs(L[i][i])<tol1){
      printf("The pivot is too small at step #\n",i);}
    s=0.;
    for (j=0; j<=i-1; j++){
      s+=L[i][j]*x[j]; }
    x[i]=(b[i]-s)/L[i][i];
  }
  return;
}
/* function */
mvsub(A,b,c,n)
double A[N][N], b[N], c[N];
int n;
{
/* 計算部分 */
  int i, j;
  double s;
  for (i=0; i<n; i++){
    s=0.0;
    for (j=0; j<n; j++){
      s+=A[i][j]*b[j];
    }
    c[i]=s;
  }
}