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

int main ( );
void gauss_seidel_test01 ( );

# include "gauss_seidel.h"

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

int main ( )

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

    gauss_seidel_test() tests gauss_seidel().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 September 2022

  Author:

    John Burkardt
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "gauss_seidel_test():\n" );
  printf ( "  C version.\n" );
  printf ( "  Test gauss_seidel().\n" );

  gauss_seidel_test01 ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "gauss_seidel_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void gauss_seidel_test01 ( )

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

    gauss_seidel_test01() tests gauss_seidel1().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 September 2022

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int i;
  int it;
  int it_num;
  int n;
  double r;
  double t;
  double *x;
  double *x_exact;
  double *x_old;

  printf ( "\n" );
  printf ( "gauss_seidel_test01():\n" );

  it_num = 2000;
  n = 33;
/*
  Set the matrix A.
*/
  a = dif2 ( n, n );
/*
  Determine the right hand side vector B.
*/
  x_exact = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    t = ( double ) i / ( double ) ( n - 1 );
    x_exact[i] = exp ( t ) * ( t - 1 ) * t;
  }
  b = r8mat_mv_new ( n, n, a, x_exact );
/*
  Carry out the iteration.
*/
  x_old = ( double * ) malloc ( n * sizeof ( double ) );
  for ( it = 0; it <= it_num; it++ )
  {
    if ( it == 0 )
    {
      x = ( double * ) malloc ( n * sizeof ( double ) );
      for ( i = 0; i < n; i++ )
      {
        x[i] = 0.0;
      }
    }
    else
    {
      x = gauss_seidel1 ( n, a, b, x_old );
    }
    r = r8mat_residual_norm ( n, n, a, x, b );
    if ( ( it <= 10 ) | ( ( it % 20 ) == 0 ) )
    {
      printf ( "  %3d  %g\n", it, r );
    }
    for ( i = 0; i < n; i++ )
    {
      x_old[i] = x[i];
    }
  }
/*
  Free memory.
*/
  free ( x );
  free ( x_old );

  return;
}