#include <stdio.h>
#include <math.h>
/* sor.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 bb[N]={4.,3.,3.,3.,4.}, xx[N]={1.,0.,0.,0.,0.};
  double ome[N]={0.8,0.9,1.,1.1,1.2}, tol[2]={0.000001, 0.0};
  double x[N], b[N], omega;
  int n=5, i, k, itr[2]={20,1};
  for (k=0; k<n; k++){
    for (i=0; i<n; i++){
      x[i]=xx[i]; b[i]=bb[i]; }
    itr[1]=1; tol[1]=0.; omega=ome[k];
    sor(A,b,x,tol,omega,n,itr);
    printf("%f %d %f\n", ome[k], itr[1], tol[1]);
    for (i=0; i<n; i++){
      printf("%f\n", x[i]); }
  }
}
/* function */
sor(A,b,x,tol,omega,n,itr)
double A[N][N], b[N], x[N], tol[2], omega;
int n, itr[2];
{
  double L[N][N],U[N][N],D[N][N],C[N][N],xnew[N],w[N],y[N];
  double ep,ep1,ep2;
  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,y,n);
    for (i=0; i<n; i++){
      w[i]=b[i]-y[i]; y[i]=0.; }
    forward(C,w,y,n);
    for (i=0; i<n; i++){
      xnew[i]=(1.0-omega)*x[i]+omega*y[i]; }
    ep1=0.; ep2=0.;
    for (i=0; i<n; i++){
      ep1+=(xnew[i]-x[i])*(xnew[i]-x[i]);
      ep2+=x[i]*x[i]; }
    ep=sqrt(fabs(ep1/ep2));
    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;
  }
}