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

using namespace std;

int main ( );
void jacobi_test01 ( );

# include "jacobi.hpp"

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_test() tests jacobi().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    14 April 2020
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "jacobi_test():\n";
  cout << "  C++ version.\n";
  cout << "  Test jacobi().\n";

  jacobi_test01 ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "jacobi_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void jacobi_test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    jacobi_test01() tests jacobi1().
//
//  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;

  cout << "\n";
  cout << "jacobi_test01():\n";

  it_num = 2000;
  n = 33;
//
//  Set the matrix A.
//
  a = dif2 ( n, n );
//
//  Determine the right hand side vector B.
//
  x_exact = new double[n];
  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.
//
  for ( it = 0; it <= it_num; it++ )
  {
    if ( it == 0 )
    {
      x = new double[n];
      for ( i = 0; i < n; i++ )
      {
        x[i] = 0.0;
      }
      x_old = new double[n];
    }
    else
    {
      x = jacobi1 ( 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];
    }
    delete [] x;
  }

  r8vec_print ( n, x, "Solution" );

  return;
}