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

int main ( );
void jacobi ( float *b, float *x1, float *x2, int n );

int main ( ) 
/*
  Purpose:

    C program for the Jacobi iterative solution of A*x=b.

  Modified:

    11 March 2017
*/
{
  float *b;
  clock_t clock1;
  clock_t clock2;
  float e;
  float e_norm;
  float *x1;
  float *x2;
  int i;
  int it;
  int n = 10000; 
  double value;

  printf ( "\n" );
  printf ( "JACOBI:\n" );
  printf ( "  C version.\n" );
  printf ( "  Apply Jacobi's iterative method to solve A*x=b\n" );
  printf ( "  System size is N = %d\n", n );
/*
  Allocate memory.
*/
  b = ( float * ) malloc ( n * sizeof ( float ) );
  x1 = ( float * ) malloc ( n * sizeof ( float ) );
  x2 = ( float * ) malloc ( n * sizeof ( float ) );
/*
  Assign B and X.
*/
  for ( i = 0; i < n; i++ )
  {
    x1[i] = ( float ) rand ( ) / ( float ) RAND_MAX;
  }

  b[0] = 2.0 * x1[0] - x1[1];
  for ( i = 1; i < n - 1; i++ )
  {
    b[i] = - x1[i-1] + 2.0 * x1[i] - x1[i+1];
  }
  b[n-1] = - x1[n-2] + 2.0 * x1[n-1];

  for ( i = 0; i < n; i++ )
  {
    x1[i] = 0.0;
  }
/*
  Call the function that does the work.
*/
  clock1 = clock ( );
  for ( it = 0; it < 10 * n; it++ )
  {
    jacobi ( b, x1, x2, n ); 
    jacobi ( b, x2, x1, n ); 
  }
  clock2 = clock ( );
  value = ( double ) ( clock2 - clock1 ) / ( double ) CLOCKS_PER_SEC;
/*
  Print a few values.
*/
  printf ( "\n" );
  printf ( "    I    B[i]      X[i]     Resid[I]\n" );
  printf ( "\n" );

  for ( i = 0; i < n && i < 10; i++ )
  {
    e = 2.0 * x1[i] - b[i];
    if ( 0 < i )
    {
      e = e - x1[i-1];
    }
    if ( i < n - 1 )
    {
      e = e - x1[i+1];
    }
    printf ( "  %3d  %8g  %8g  %8g\n", i, b[i], x1[i], e );
  }
/*
  Compute the error.
*/
  e_norm = 0.0;
  for ( i = 0; i < n; i++ )
  {
    e = 2.0 * x1[i] - b[i];
    if ( 0 < i )
    {
      e = e - x1[i-1];
    }
    if ( i < n - 1 )
    {
      e = e - x1[i+1];
    }
    e_norm = e_norm + e * e;
  }
  e_norm = sqrt ( e_norm / ( float ) ( n ) );
  printf ( "\n" );
  printf ( "  RMS error = %g\n", e_norm );
  printf ( "\n" );
  printf ( "  Computation required %g seconds\n", value );
/*
  Free memory.
*/
  free ( b );
  free ( x1 );
  free ( x2 );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "JACOBI:\n" );
  printf ( "  Normal end of execution.\n" );

  return 0;
}

void jacobi ( float *b, float *x1, float *x2, int n ) 
{
/*
  Purpose:

    JACOBI carries out one thread of one step of Jacobi iteration.

  Discussion:

    Specifically, for 0 < i < N - 1
      - x(i-1) + 2 x(i) - x(i+1) = b(i) 
    is rewritten
      x2(i) = ( b(i) + x(i-1) + x(i+1) ) / 2
    for i = 0, we have
      x2(i) = ( b(i) + x(i+1) ) / 2
    and for i = n - 1 we have
      x2(i) = ( b(i) + x(i-1) ) / 2

  Modified:

    11 March 2017
*/
  int i;
  float xim1;
  float xip1;

  for ( i = 0; i < n; i++ )
  {
    if ( 0 < i )
    {
      xim1 = x1[i-1];
    }
    else
    {
      xim1 = 0.0;
    }
    if ( i < n - 1 )
    {
      xip1 = x1[i+1];
    }
    else
    {
      xip1 = 0.0;
    }
    x2[i] = 0.5 * ( b[i] + xim1 + xip1 );

  }
  return;
}