#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/*******************************************************************/

void print_matrix(double **A, int n) {
    int i, j;
    for(i = 0; i < n ; i++) {
       for(j = 0; j < n; j++)
          printf("% 6.2e  ",A[i][j]);
       printf("\n");
    }
    printf("\n");
}

/*******************************************************************/

void print_vector(double *c, int n) {
   int i;
   for(i = 0; i < n ; i++)
      printf("% 6.2e  ", c[i]);
   printf("\n");
}

/*******************************************************************/

int main() {
   double **a , *d, *x, *x_old, *e, s;
   double om = 1.7;       /* omega */
   int    i, j, k, n;
   double eps=1.0e-10;    /* epsilon */
   double n_rhs, n_lhs;   /* norm right- and left hand side of inequality */ 

   n = 8;   /* dimension of the problem */

   a = (double **)malloc( n * sizeof(double*) );
   d = (double *) malloc( n * sizeof(double) );   
   x = (double *) malloc( n * sizeof(double) ); 
   e = (double *) malloc( n * sizeof(double) );
   x_old = (double *) malloc( n * sizeof(double) ); 

   for(i=0; i < n; i++)
      a[i] = (double *)malloc( n * sizeof(double) );

   /* initialization */
   /* srand( (unsigned)time( (long *)0 )); */
   for(i = 0; i < n ; i++) {
      for(j = 0; j < n ; j++)
         a[i][j] = rand() /2147483647.0;
      a[i][i] += 11;  
      d[i] = rand() /2147483647.0;   
      x[i] = 0.0;                  /* start vector */
   }

   /* print_matrix(a, n); */

   do {
      memcpy( x_old, x, n * sizeof(double) );
      n_lhs = 0.0 ; 
      n_rhs = 0.0 ; 
      for(i = 0; i < n; i++) 
	 n_rhs += x_old[i] * x_old[i];
      n_rhs = sqrt( n_rhs );
      for(i = 0; i < n; i++) {    
         s = d[i];
         for(j = 0; j < n; j++) 
            s += ( (i!=j) ? (-a[i][j]) : ((1-om)/om * a[i][j]) ) * x[j];
         x[i] = om / a[i][i] * s;
	 n_lhs += (x[i]-x_old[i]) * (x[i]-x_old[i]);
      }
      n_lhs = sqrt( n_lhs );
      printf("n_lhs = % 6.2e   n_rhs = % 6.2e\n", n_lhs, n_rhs);
   } while ( n_lhs > eps * n_rhs );

   for(i = 0; i < n; i++) {     /* calculation of the error */
      e[i] = 0.0; 
      for(j = 0; j < n; j++) 
         e[i] += a[i][j]*x[j];
      e[i] = d[i]-e[i];
   }
   printf("\n");
   print_vector(x,n);           /* printing of x vector */
   print_vector(e,n);           /* printing of error vector */
}

/*******************************************************************/