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

# include "sor.h"

int main ( );
void sor_test01 ( double w );
void timestamp ( );

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

int main ( )

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

    sor_test() tests sor().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 August 2024

  Author:

    John Burkardt
*/
{
  double w;

  timestamp ( );
  printf ( "\n" );
  printf ( "sor_test():\n" );
  printf ( "  C version.\n" );
  printf ( "  Test sor().\n" );

  w = 1.5;
  sor_test01 ( w );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "sor_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void sor_test01 ( double w )

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

    sor_test01() tests sor1().

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    14 August 2024

  Author:

    John Burkardt

  Input:

    double W, the relaxation factor.
    0 < W < 2 is required.
*/
{
  double *a;
  double *b;
  FILE *command;
  FILE *data;
  int i;
  int it;
  int it_num;
  int j;
  double *m_plot;
  double *ml_plot;
  int n;
  double *r_plot;
  double *rl_plot;
  double *s_plot;
  double t;
  double *x;
  double *x_exact;
  double *x_new;
  double *x_plot;

  printf ( "\n" );
  printf ( "sor1_test01():\n" );
  printf ( "  Relaxation parameter W = %g\n", w );

  it_num = 2000;
  n = 33;
/*
  Set the matrix A.
*/
  a = dif2_matrix ( 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 );
/*
  Set the initial estimate for the solution.
*/
  x = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    x[i] = 0.0;
  }
/*
  Allocate plot arrays.
*/
  m_plot = ( double * ) malloc ( (it_num+1) * sizeof ( double ) );
  r_plot = ( double * ) malloc ( (it_num+1) * sizeof ( double ) );
  s_plot = ( double * ) malloc ( (it_num+1) * sizeof ( double ) );
  x_plot = ( double * ) malloc ( (n*(it_num+1)) * sizeof ( double ) );
/*
  Initialize plot arrays.
*/
  r_plot[0] = r8mat_residual_norm ( n, n, a, x, b );
  m_plot[0] = 1.0;
  for ( i = 0; i < n; i++ )
  {
    x_plot[i+0*n] = x[i];
  }
  for ( j = 0; j <= it_num; j++ )
  {
    s_plot[j] = ( double ) j;
  }
/*
  Carry out the iteration.
*/
  for ( it = 1; it <= it_num; it++ )
  {
    x_new = sor1 ( n, a, b, x, w );

    r_plot[it] = r8mat_residual_norm ( n, n, a, x_new, b );
/*
  Compute the average point motion.
*/
    m_plot[it] = r8vec_diff_norm_squared ( n, x, x_new ) / ( double ) n;
/*
  Update the solution
*/
    r8vec_copy ( n, x_new, x );

    for ( i = 0; i < n; i++ )
    {
      x_plot[i+0*n] = x[i];
    }

    free ( x_new );
  }
    r8vec_print ( n, x, "Solution" );
/*
  Plot the residual.
*/
  rl_plot = ( double * ) malloc ( (it_num+1) * sizeof ( double ) );
  for ( j = 0; j <= it_num; j++ )
  {
    rl_plot[j] = log ( r_plot[j] );
  }
/*
  Create the data file.
*/
  data = fopen ( "residual_data.txt", "wt" );

  for ( j = 0; j <= it_num; j++ )
  {
    fprintf ( data, "%d  %g\n", j, rl_plot[j] );
  }
  fclose ( data );

  printf ( " \n" );
  printf ( "  Data stored in '%s'\n", "residual_data.txt" );
/*
  Create the command file.
*/
  command = fopen ( "residual_commands.txt", "wt" );

  fprintf ( command, "# residual_commands.txt\n" );
  fprintf ( command, "#\n" );
  fprintf ( command, "# Usage:\n" );
  fprintf ( command, "#  gnuplot < residual_commands.txt\n" );
  fprintf ( command, "#\n" );
  fprintf ( command, "set term png\n" );
  fprintf ( command, "set output 'residual.png'\n" );
  fprintf ( command, "set style data lines\n" );
  fprintf ( command, "set xlabel 'Iteration'\n" );
  fprintf ( command, "set ylabel 'Residual'\n" );
  fprintf ( command, "set title 'Log(Residual) over Iterations'\n" );
  fprintf ( command, "set grid\n" );
  fprintf ( command, "plot 'residual_data.txt' using 1:2 lw 2\n" );
  fprintf ( command, "quit\n" );

  fclose ( command );

  printf ( "  Plot commands stored in '%s'\n", "residual_commands.txt" );
/*
  Plot the average point motion.
*/
  ml_plot = ( double * ) malloc ( (it_num+1) * sizeof ( double ) );
  for ( j = 0; j <= it_num; j++ )
  {
    ml_plot[j] = log ( m_plot[j] );
  }
/*
  Create the data file.
*/
  data = fopen ( "motion_data.txt", "wt" );

  for ( j = 0; j <= it_num; j++ )
  {
    fprintf ( data, "%d  %g\n", j, ml_plot[j] );
  }
  fclose ( data );

  printf ( " \n" );
  printf ( "  Data stored in '%s'\n", "motion_data.txt" );
/*
  Create the command file.
*/
  command = fopen ( "motion_commands.txt", "wt" );

  fprintf ( command, "# motion_commands.txt\n" );
  fprintf ( command, "#\n" );
  fprintf ( command, "# Usage:\n" );
  fprintf ( command, "#  gnuplot < motion_commands.txt\n" );
  fprintf ( command, "#\n" );
  fprintf ( command, "set term png\n" );
  fprintf ( command, "set output 'motion.png'\n" );
  fprintf ( command, "set style data lines\n" );
  fprintf ( command, "set xlabel 'Iteration'\n" );
  fprintf ( command, "set ylabel 'Motion'\n" );
  fprintf ( command, "set title 'Log(Motion) over Iterations'\n" );
  fprintf ( command, "set grid\n" );
  fprintf ( command, "plot 'motion_data.txt' using 1:2 lw 2\n" );
  fprintf ( command, "quit\n" );

  fclose ( command );

  printf ( "  Plot commands stored in '%s'\n", "motion_commands.txt" );
/*
  Free memory.
*/
  free ( a );
  free ( b );
  free ( m_plot );
  free ( ml_plot );
  free ( r_plot );
  free ( rl_plot );
  free ( s_plot );
  free ( x );
  free ( x_exact );
  free ( x_plot );

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

void timestamp ( )

/******************************************************************************/
/*
  Purpose:
 
    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 May 2021

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}