#include <stdio.h>
#include <stdlib.h>
#include <time.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, *e, s;
   double om = 1.7;      /* omega */
   int    i, j, k, n;
   int    k_stop = 100; 

   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) );

   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); */

   for(k = 0; k <= k_stop  ; k++) {
      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;
      }
   }

   for(i = 0; i < n; i++) {
      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(e,n);  
}

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

