# include <cmath>
# include <cstdlib>
# include <cstring>
# include <iomanip>
# include <iostream>

using namespace std;

# include "poisson_1d.hpp"

//****************************************************************************80

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

//****************************************************************************80
//
//  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 estimated solution.
//
//    double &DIF_L1, the L1 norm of the difference between the
//    input and output 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;
}
//****************************************************************************80

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[] )

//****************************************************************************80
//
//  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. 
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 October 2024
//
//  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 endpoints.
//
//    double UA, UB, the boundary values at the endpoints.
//
//    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 = new double[n+1];
  for ( i = 0; i <= n; i++ )
  {
    x[i] = ( ( n - i ) * a + i * b ) / n;
  }
//
//  Set the right hand side.
//
  r = new double[n+1];

  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.
//
  u[0] = ua;
  for ( i = 1; i < n; i++ )
  {
    u[i] = 0.0;
  }
  u[n] = ub;
//
//  Gauss-Seidel iteration.
//
  it_num = 0;

  for ( ; ; )
  {
    it_num = it_num + 1;

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

    if ( d1 <= tol )
    {
      break;
    }
  }

  delete [] r;
  delete [] x;

  return;
}