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

# include "poisson_1d.h"

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

void gauss_seidel ( int n, double r[], double u[], double *dif_l1 )

/******************************************************************************/
/*
  Purpose:

    gauss_seidel() carries out one step of a Gauss-Seidel iteration.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    07 December 2011

  Author:

    John Burkardt

  Reference:

    William Hager,
    Applied Numerical Linear Algebra,
    Prentice-Hall, 1988,
    ISBN13: 978-0130412942,
    LC: QA184.H33.

  Input:

    int N, the number of unknowns.

    double R[N], the right hand side.

    double U[N], the estimated solution.

  Output:

    double U[N], the updated solution.

    double *DIF_L1, the L1 norm of the change in the solution estimates.
*/
{
  int i;
  double u_old;

  *dif_l1 = 0.0;

  for ( i = 1; i < n - 1; i++ )
  {
    u_old = u[i];
    u[i] = 0.5 * ( u[i-1] + u[i+1] + r[i] );
    *dif_l1 = *dif_l1 + fabs ( u[i] - u_old );
  }

  return;
}
/******************************************************************************/

void poisson_1d ( int n, double a, double b, double ua, double ub,
  double force ( double x ), double exact ( double x ), int *it_num, 
  double u[] )

/******************************************************************************/
/*
  Purpose:

    poisson_1d() solves a 1D PDE, using the Gauss-Seidel method.

  Discussion:

    This routine solves a 1D boundary value problem of the form

      - U''(X) = F(X) for A < X < B,

    with boundary conditions U(A) = UA, U(B) = UB.

    The Gauss-Seidel method is used to solve the corresponding linear system.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 July 2014

  Author:

    John Burkardt

  Reference:

    William Hager,
    Applied Numerical Linear Algebra,
    Prentice-Hall, 1988,
    ISBN13: 978-0130412942,
    LC: QA184.H33.

  Input:

    int N, the number of intervals.

    double A, B, the left and right endpoints.

    double UA, UB, the left and right boundary values.

    double FORCE ( double x ), the name of the function 
    which evaluates the right hand side.

    double EXACT ( double x ), the name of the function 
    which evaluates the exact solution.

  Output:

    int *IT_NUM, the number of iterations.

    double U[N+1], the computed solution.
*/
{
  double d1;
  double h;
  int i;
  double *r;
  double tol;
  double *x;
/*
  Initialization.
*/
  tol = 0.0001;
/*
  Set the nodes.
*/
  x = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );
  for ( i = 0; i <= n; i++ )
  {
    x[i] = ( ( n - i ) * a + i * b ) / n;
  }
/*
  Set the right hand side.
*/
  r = ( double * ) malloc ( ( n + 1 ) * sizeof ( double ) );

  r[0] = ua;
  h = ( b - a ) / ( double ) ( n );
  for ( i = 1; i < n; i++ )
  {
    r[i] = h * h * force ( x[i] );
  }
  r[n] = ub;
/*
  Initialize the solution estimate.
*/
  for ( i = 0; i <= n; i++ )
  {
    u[i] = 0.0;
  }
  *it_num = 0;
/*
  Gauss Seidel iteration to improve U.
*/
  for ( ; ; )
  {
    *it_num = *it_num + 1;

    gauss_seidel ( n + 1, r, u, &d1 );

    if ( d1 <= tol )
    {
      break;
    }
  }
/*
  Free memory.
*/
  free ( r );
  free ( x );

  return;
}