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

# include "r83_np.h"

int main ( );
void r83_np_det_test ( );
void r83_np_fa_test ( );
void r83_np_fs_test ( );
void r83_np_ml_test ( );
void r83_np_sl_test ( );
void timestamp ( );

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

int main ( )

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

    r83_np_test() tests r83_np().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    15 August 2022

  Author:

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

  r83_np_det_test ( );
  r83_np_fa_test ( );
  r83_np_fs_test ( );
  r83_np_ml_test ( );
  r83_np_sl_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "r83_np_test():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r83_np_det_test ( )

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

    r83_np_det_test() tests r83_np_det().

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    20 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double det;
  int info;
  int n = 10;

  printf ( "\n" );
  printf ( "r83_np_det_test():\n" );
  printf ( "  r83_np_det() computes the determinant of a tridiagonal matrix\n" );
  printf ( "  after it has been factored by R83_NP_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
/*
  Set the matrix.
*/
  a = r83_dif2 ( n, n );

  r83_print ( n, n, a, "  The R83 matrix:" );
/*
  Factor the matrix.
*/
  info = r83_np_fa ( n, a );

  r83_print ( n, n, a, "  The factored matrix:" );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R83_NP_DET_TEST - Warning!\n" );
    printf ( "  R83_NP_FA returns INFO = %d\n", info );
    return;
  }
/*
  Compute the determinant.
*/
  det = r83_np_det ( n, a );

  printf ( "\n" );
  printf ( "  R83_NP_DET computes determinant =  %g\n", det );
  printf ( "  Exact determinant = %d\n", n + 1 );

  free ( a );

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

void r83_np_fa_test ( )

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

    R83_NP_FA_TEST tests R83_NP_FA.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    05 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  int info;
  int job;
  int seed = 123456789;
  double *x;

  printf ( "\n" );
  printf ( "r83_np_fa_test():\n" );
  printf ( "  r83_np_fa() factors a tridiagonal matrix with no pivoting.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r83_random ( N, N, &seed );

  r83_print ( N, N, a, "  The tridiagonal matrix:" );
/*
  Set the desired solution.
*/
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  b = r83_mv ( N, N, a, x );
  free ( x );
/*
  Factor the matrix.
*/
  info = r83_np_fa ( N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R83_NP_FA_TEST - Fatal error!\n" );
    printf ( "  The test matrix is singular.\n" );
    return;
  }
/*
  Solve the linear system.
*/
  job = 0;
  x = r83_np_sl ( N, a, b, job );
 
  r8vec_print ( N, x, "  Solution to A*x=b:" );
/*
  Set the desired solution
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side, using the factored matrix.
*/
  job = 1;
  free ( b );
  b = r83_np_ml ( N, a, x, job );
/*
  Solve the linear system.
*/
  job = 1;
  free ( x );
  x = r83_np_sl ( N, a, b, job );
 
  r8vec_print ( N, x, "  Solution to A'*x=b:" );
 
  free ( a );
  free ( b );
  free ( x );

  return;
# undef N
}
/******************************************************************************/

void r83_np_fs_test ( )

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

    R83_NP_FS_TEST tests R83_NP_FS.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    04 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  int seed = 123456789;
  double *x;

  printf ( "\n" );
  printf ( "r83_np_fs_test():\n" );
  printf ( "  r83_np_fs() factors and solves a tridiagonal\n" );
  printf ( "  linear system.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix elements.
*/
  a = r83_random ( N, N, &seed );
/*
  Set the desired solution.
*/
  x = r8vec_indicator1_new ( N );
/*
  Compute b = A * x.
*/
  b = r83_mv ( N, N, a, x );
/*
  Wipe out the solution.
*/
  free ( x );
/*
  Solve the system.
*/
  x = r83_np_fs ( N, a, b );

  r8vec_print ( N, x, "  Solution to A*x=b:" );

  free ( a );
  free ( b );
  free ( x );

  return;
# undef N
}
/******************************************************************************/

void r83_np_ml_test ( )

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

    R83_NP_ML_TEST tests R83_NP_ML.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    04 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  double *b2;
  int info;
  int job;
  int seed = 123456789;
  double *x;

  printf ( "\n" );
  printf ( "r83_np_ml_test():\n" );
  printf ( "  r83_np_ml() computes A*x or A'*x\n" );
  printf ( "  where A has been factored by R83_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    a = r83_random ( N, N, &seed );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    if ( job == 0 )
    {
      b = r83_mv ( N, N, a, x );
    }
    else
    {
      b = r83_mtv ( N, N, a, x );
    }
/*
  Factor the matrix.
*/
    info = r83_np_fa ( N, a );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R83_NP_ML_TEST - Fatal error!\n" );
      printf ( "  R83_NP_FA declares the matrix is singular!\n" );
      printf ( "  The value of INFO is %d\n", info );
      return;
    }
/*
  Now multiply factored matrix times solution to get right hand side again.
*/
    b2 = r83_np_ml ( N, a, x, job );

    if ( job == 0 )
    {
      r8vec2_print_some ( N, b, b2, 10, "  A*x and PLU*x:" );
    }
    else
    {
      r8vec2_print_some ( N, b, b2, 10, "  A'*x and (PLU)'*x" );
    }

    free ( a );
    free ( b );
    free ( b2 );
    free ( x );
  }

  return;
# undef N
}
/******************************************************************************/

void r83_np_sl_test ( )

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

    R83_NP_SL_TEST tests R83_NP_Sl.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    05 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  int info;
  int job;
  int seed = 123456789;
  double *x;

  printf ( "\n" );
  printf ( "r83_np_sl_test():\n" );
  printf ( "  r83_np_sl() solves a linear system that was factored by R83_NP_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r83_random ( N, N, &seed );

  r83_print ( N, N, a, "  The tridiagonal matrix:" );
/*
  Set the desired solution.
*/
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  b = r83_mv ( N, N, a, x );
  free ( x );
/*
  Factor the matrix.
*/
  info = r83_np_fa ( N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R83_NP_SL_TEST - Fatal error!\n" );
    printf ( "  The test matrix is singular.\n" );
    return;
  }
/*
  Solve the linear system.
*/
  job = 0;
  x = r83_np_sl ( N, a, b, job );
 
  r8vec_print ( N, x, "  Solution to A*x=b:" );
/*
  Set the desired solution
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side, using the factored matrix.
*/
  job = 1;
  free ( b );
  b = r83_np_ml ( N, a, x, job );
/*
  Solve the linear system.
*/
  job = 1;
  free ( x );
  x = r83_np_sl ( N, a, b, job );
 
  r8vec_print ( N, x, "  Solution to A'*x=b:" );
 
  free ( a );
  free ( b );
  free ( x );

  return;
# undef N
}
/******************************************************************************/

void timestamp ( )

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

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    31 May 2001 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    24 September 2003

  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
}