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

# include "r83p.h"

int main ( );
void r83p_det_test ( );
void r83p_fa_test ( );
void r83p_indicator_test ( );
void r83p_ml_test ( );
void r83p_mtv_test ( );
void r83p_mv_test ( );
void r83p_print_test ( );
void r83p_print_some_test ( );
void r83p_random_test ( );
void r83p_sl_test ( );
void r83p_to_r8ge_test ( );
void r83p_zeros_test ( );

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

int main ( )

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

    MAIN is the main program for LINPLUS_TEST.

  Discussion:

    LINPLUS_TEST tests LINPLUS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 May 2016

  Author:

    John Burkardt
*/
{
  timestamp ( );
  printf ( "\n" );
  printf ( "LINPLUS_TEST\n" );
  printf ( "  C version\n" );
  printf ( "  Test the LINPLUS library.\n" );

  r83p_det_test ( );
  r83p_fa_test ( );
  r83p_indicator_test ( );
  r83p_ml_test ( );
  r83p_mtv_test ( );
  r83p_mv_test ( );
  r83p_print_test ( );
  r83p_print_some_test ( );
  r83p_random_test ( );
  r83p_sl_test ( );
  r83p_to_r8ge_test ( );
  r83p_zeros_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "LINPLUS_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r83p_det_test ( )

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

    R83P_DET_TEST tests R83P_DET.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 March 2013

  Author:

    John Burkardt
*/
{
# define N 12

  double *a;
  double *b;
  double det;
  int info;
  int pivot[N];
  int seed = 123456789;
  double work2[N-1];
  double work3[N-1];
  double work4;

  printf ( "\n" );
  printf ( "R83P_DET_TEST\n" );
  printf ( "  R83P_DET, determinant of a tridiagonal\n" );
  printf ( "  periodic matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r83p_random ( N, &seed );

  r83p_print ( N, a, "  The periodic tridiagonal matrix:" );
/*
  Copy the matrix into a general array.
*/
  b = r83p_to_r8ge ( N, a );
/*
  Factor the matrix.
*/
  info = r83p_fa ( N, a, work2, work3, &work4 );
/*
  Compute the determinant.
*/
  det = r83p_det ( N, a, work4 );

  printf ( "\n" );
  printf ( "  R83P_DET computes the determinant = %g\n", det );
/*
  Factor the general matrix.
*/
  info = r8ge_fa ( N, b, pivot );
/*
  Compute the determinant.
*/
  det = r8ge_det ( N, b, pivot );

  printf ( "  R8GE_DET computes the determinant = %g\n", det );

  free ( a );
  free ( b );

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

void r83p_fa_test ( )

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

    R83P_FA_TEST tests R83P_FA, R83P_SL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  int i;
  int info;
  int job;
  int seed = 123456789;
  double work2[N-1];
  double work3[N-1];
  double work4;
  double *x;

  printf ( "\n" );
  printf ( "R83P_FA_TEST\n" );
  printf ( "  R83P_FA factors a tridiagonal periodic system.\n" );
  printf ( "  R83P_SL solves a tridiagonal periodic system.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r83p_random ( N, &seed );
/*
  Factor the matrix.
*/
  info = r83p_fa ( N, a, work2, work3, &work4 );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R83P_FA_TEST - Fatal error!\n" );
    printf ( "  R83P_FA returns INFO = %d\n", info );
    return;
  }

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    b = r83p_ml ( N, a, x, job );
/*
  Solve the linear system.
*/
    free ( x );

    x = r83p_sl ( N, a, b, job, work2, work3, work4 );

    if ( job == 0 )
    {
      r8vec_print ( N, x, "  Solution to A*x=b:" );
    }
    else
    {
      r8vec_print ( N, x, "  Solution to A'*x=b:" );
    }

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

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

void r83p_indicator_test ( )

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

    R83P_INDICATOR_TEST tests R83P_INDICATOR.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;

  printf ( "\n" );
  printf ( "R83P_INDICATOR_TEST\n" );
  printf ( "  R83P_INDICATOR sets up an R83P indicator matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_indicator ( n );

  r83p_print ( n, a, "  The R83P indicator matrix:" );

  free ( a );

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

void r83p_ml_test ( )

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

    R83P_ML_TEST tests R83P_ML.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    07 March 2013

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  double *b2;
  int info;
  int i;
  int job;
  int seed = 123456789;
  double work2[N-1];
  double work3[N-1];
  double work4;
  double *x;

  printf ( "\n" );
  printf ( "R83P_ML_TEST\n" );
  printf ( "  R83P_ML computes A*x or A'*X\n" );
  printf ( "  where A has been factored by R83P_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    a = r83p_random ( N, &seed );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    if ( job == 0 )
    {
      b = r83p_mv ( N, a, x );
    }
    else
    {
      b = r83p_mtv ( N, a, x );
    }
/*
  Factor the matrix.
*/
    info = r83p_fa ( N, a, work2, work3, &work4 );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R83P_ML_TEST - Fatal error!\n" );
      printf ( "  R83P_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 = r83p_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 ( x );
    free ( b );
    free ( b2 );
  }

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

void r83p_mtv_test ( )

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

    R83P_MTV_TEST tests R83P_MTV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int n = 5;
  double *x;

  printf ( "\n" );
  printf ( "R83P_MTV_TEST\n" );
  printf ( "  R83P_MTV computes A'*x=b for an R83P matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_indicator ( n );

  r83p_print ( n, a, "  The R83P indicator matrix:" );

  x = r8vec_indicator1_new ( n );

  r8vec_print ( n, x, "  Vector x:" );

  b = r83p_mtv ( n, a, x );

  r8vec_print ( n, b, "  b=A'*x:" );

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

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

void r83p_mv_test ( )

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

    R83P_MV_TEST tests R83P_MV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int n = 5;
  double *x;

  printf ( "\n" );
  printf ( "R83P_MV_TEST\n" );
  printf ( "  R83P_MV computes A*x=b for an R83P matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_indicator ( n );

  r83p_print ( n, a, "  The R83P indicator matrix:" );

  x = r8vec_indicator1_new ( n );

  r8vec_print ( n, x, "  Vector x:" );

  b = r83p_mv ( n, a, x );

  r8vec_print ( n, b, "  b=A*x:" );

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

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

void r83p_print_test ( )

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

    R83P_PRINT_TEST tests R83P_PRINT.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;

  printf ( "\n" );
  printf ( "R83P_PRINT_TEST\n" );
  printf ( "  R83P_PRINT prints an R83P matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_indicator ( n );

  r83p_print ( n, a, "  The R83P matrix:" );

  free ( a );

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

void r83p_print_some_test ( )

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

    R83P_PRINT_SOME_TEST tests R83P_PRINT_SOME.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

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

  printf ( "\n" );
  printf ( "R83P_PRINT_SOME_TEST\n" );
  printf ( "  R83P_PRINT_SOME prints some of an R83P matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_indicator ( n );

  r83p_print_some ( n, a, 1, 1, n, 2, "  Rows 1:N, Cols 1:2:" );

  free ( a );

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

void r83p_random_test ( )

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

    R83P_RANDOM_TEST tests R83P_RANDOM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int seed;

  printf ( "\n" );
  printf ( "R83P_RANDOM_TEST\n" );
  printf ( "  R83P_RANDOM sets up a random R83P matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  seed = 123456789;
  a = r83p_random ( n, &seed );

  r83p_print ( n, a, "  The random R83P matrix:" );

  free ( a );

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

void r83p_sl_test ( )

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

    R83P_SL_TEST tests R83P_SL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
# define N 10

  double *a;
  double *b;
  int i;
  int info;
  int job;
  int seed = 123456789;
  double work2[N-1];
  double work3[N-1];
  double work4;
  double *x;

  printf ( "\n" );
  printf ( "R83P_SL_TEST\n" );
  printf ( "  R83P_SL solves a tridiagonal periodic system\n" );
  printf ( "  after it has been factored by R83P_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r83p_random ( N, &seed );
/*
  Factor the matrix.
*/
  info = r83p_fa ( N, a, work2, work3, &work4 );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R83P_SL_TEST - Fatal error!\n" );
    printf ( "  R83P_FA returns INFO = %d\n", info );
    return;
  }

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    b = r83p_ml ( N, a, x, job );
/*
  Solve the linear system.
*/
    free ( x );

    x = r83p_sl ( N, a, b, job, work2, work3, work4 );

    if ( job == 0 )
    {
      r8vec_print ( N, x, "  Solution to A*x=b:" );
    }
    else
    {
      r8vec_print ( N, x, "  Solution to A'*x=b:" );
    }

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

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

void r83p_to_r8ge_test ( )

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

    R83P_TO_R8GE_TEST tests R83P_TO_R8GE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a_r83p;
  double *a_r8ge;
  int n = 5;
  int seed;

  printf ( "\n" );
  printf ( "R83P_TO_R8GE_TEST\n" );
  printf ( "  R83P_TO_R8GE convert a matrix from R83P to R8GE format.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  seed = 123456789;
  a_r83p = r83p_random ( n, &seed );
  r83p_print ( n, a_r83p, "  The R83P matrix:" );

  a_r8ge = r83p_to_r8ge ( n, a_r83p );
  r8ge_print ( n, n, a_r8ge, "  The R8GE matrix:" );

  free ( a_r83p );
  free ( a_r8ge );

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

void r83p_zeros_test ( )

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

    R83P_ZEROS_TEST tests R83P_ZEROS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    15 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;

  printf ( "\n" );
  printf ( "R83P_ZEROS_TEST\n" );
  printf ( "  R83P_ZEROS sets up a R83P zero matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );

  a = r83p_zeros ( n );

  r83p_print ( n, a, "  The R83P zero matrix:" );

  free ( a );

  return;
}