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

# include "r8cb.h"

int main ( );

void r8cb_det_test ( );
void r8cb_dif2_test ( );
void r8cb_indicator_test ( );
void r8cb_ml_test ( );
void r8cb_mtv_test ( );
void r8cb_mv_test ( );
void r8cb_np_fa_test ( );
void r8cb_np_sl_test ( );
void r8cb_print_test ( );
void r8cb_print_some_test ( );
void r8cb_random_test ( );
void r8cb_to_r8ge_test ( );
void r8cb_to_r8vec_test ( );
void r8cb_zeros_test ( );
void r8vec_to_r8cb_test ( );

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

int main ( )

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

    MAIN is the main program for R8CB_TEST.

  Discussion:

    R8CB_TEST tests R8CB.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

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

  r8cb_det_test ( );
  r8cb_dif2_test ( );
  r8cb_indicator_test ( );
  r8cb_ml_test ( );
  r8cb_mtv_test ( );
  r8cb_mv_test ( );
  r8cb_np_fa_test ( );
  r8cb_np_sl_test ( );
  r8cb_print_test ( );
  r8cb_print_some_test ( );
  r8cb_random_test ( );
  r8cb_to_r8ge_test ( );
  r8cb_to_r8vec_test ( );
  r8cb_zeros_test ( );
  r8vec_to_r8cb_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "R8CB_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r8cb_det_test ( )

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

    R8CB_DET_TEST tests R8CB_DET.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    13 March 2013

  Author:

    John Burkardt
*/
{
# define N 10
# define ML 2
# define MU 3

  double *a;
  double *a2;
  double det;
  int info;
  int pivot[N];
  int seed = 123456789;

  printf ( "\n" );
  printf ( "R8CB_DET_TEST\n" );
  printf ( "  R8CB_DET computes the determinant of a matrix factored by R8CB_NP_FA;\n" );
  printf ( "\n" );
  printf ( "  Matrix order N     = %d\n", N );
  printf ( "  Lower bandwidth ML = %d\n", ML );
  printf ( "  Upper bandwidth MU = %d\n", MU );
/*
  Set the matrix.
*/
  a = r8cb_random ( N, N, ML, MU, &seed );

  r8cb_print ( N, N, ML, MU, a, "  The compact band matrix:" );
/*
  Copy the matrix into a general array.
*/
  a2 = r8cb_to_r8ge ( N, N, ML, MU, a );
/*
  Factor the matrix.
*/
  info = r8cb_np_fa ( N, ML, MU, a );
/*
  Compute the determinant.
*/
  det = r8cb_det ( N, ML, MU, a );

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

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

  free ( a );
  free ( a2 );

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

void r8cb_dif2_test ( )

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

    R8CB_DIF2_TEST tests R8CB_DIF2.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 8;

  printf ( "\n" );
  printf ( "R8CB_DIF2_TEST\n" );
  printf ( "  R8CB_DIF2 computes an R8CB second difference matrix\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_dif2 ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB second difference matrix:" );

  free ( a );

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

void r8cb_indicator_test ( )

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

    R8CB_INDICATOR_TEST tests R8CB_INDICATOR.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    13 March 2013

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8CB_INDICATOR_TEST\n" );
  printf ( "  R8CB_INDICATOR sets up an R8CB indicator matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB indicator matrix:" );

  free ( a );

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

void r8cb_ml_test ( )

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

    R8CB_ML_TEST tests R8CB_ML.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 March 2013

  Author:

    John Burkardt
*/
{
# define N 10
# define ML 1
# define MU 2

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

  printf ( "\n" );
  printf ( "R8CB_ML_TEST\n" );
  printf ( "  R8CB_ML computes A*x or A'*X\n" );
  printf ( "  where A has been factored by R8CB_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N     = %d\n", N );
  printf ( "  Lower bandwidth ML = %d\n", ML );
  printf ( "  Upper bandwidth MU = %d\n", MU );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    a = r8cb_random ( N, N, ML, MU, &seed );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    if ( job == 0 )
    {
      b = r8cb_mv ( N, ML, MU, a, x );
    }
    else
    {
      b = r8cb_mtv ( N, ML, MU, a, x );
    }
/*
  Factor the matrix.
*/
    info = r8cb_np_fa ( N, ML, MU, a );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R8CB_ML_TEST - Fatal error!\n" );
      printf ( "  R8CB_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 = r8cb_ml ( N, ML, MU, a, x, job );

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

    free ( a );
    free ( b );
    free ( b2 );
    free ( x );
  }
  return;
# undef N
# undef ML
# undef MU
}
/******************************************************************************/

void r8cb_mtv_test ( )

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

    R8CB_MTV_TEST tests R8CB_MTV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int m = 8;
  int ml = 2;
  int mu = 1;
  int n = 8;
  double *x;

  printf ( "\n" );
  printf ( "R8CB_MTV_TEST\n" );
  printf ( "  R8CB_MTV computes b=A'*x, where A is an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( n, n, ml, mu );

  r8cb_print ( n, n, ml, mu, a, "  The R8CB matrix:" );

  x = r8vec_indicator1_new ( m );

  r8vec_print ( m, x, "  The vector x:" );

  b = r8cb_mtv ( n, ml, mu, a, x );

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

  free ( a );

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

void r8cb_mv_test ( )

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

    R8CB_MV_TEST tests R8CB_MV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int m = 8;
  int ml = 2;
  int mu = 1;
  int n = 8;
  double *x;

  printf ( "\n" );
  printf ( "R8CB_MV_TEST\n" );
  printf ( "  R8CB_MV computes b=A*x, where A is an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( n, n, ml, mu );

  r8cb_print ( n, n, ml, mu, a, "  The R8CB matrix:" );

  x = r8vec_indicator1_new ( n );

  r8vec_print ( n, x, "  The vector x:" );

  b = r8cb_mv ( n, ml, mu, a, x );

  r8vec_print ( m, b, "  The product b=A*x:" );

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

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

void r8cb_np_fa_test ( )

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

    R8CB_NP_FA_TEST tests R8CB_NP_FAL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 March 2013

  Author:

    John Burkardt
*/
{
# define N 10
# define ML 1
# define MU 2

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

  printf ( "\n" );
  printf ( "R8CB_NP_FA_TEST\n" );
  printf ( "  R8CB_NP_FA factors an R8CB matrix with no pivoting;\n" );
  printf ( "\n" );
  printf ( "  Matrix order N     = %d\n", N );
  printf ( "  Lower bandwidth ML = %d\n", ML );
  printf ( "  Upper bandwidth MU = %d\n", MU );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    a = r8cb_random ( N, N, ML, MU, &seed );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the right hand side.
*/
    if ( job == 0 )
    {
      b = r8cb_mv ( N, ML, MU, a, x );
    }
    else
    {
      b = r8cb_mtv ( N, ML, MU, a, x );
    }
/*
  Factor the matrix.
*/
    info = r8cb_np_fa ( N, ML, MU, a );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R8CB_NP_FA_TEST - Fatal error!\n" );
      printf ( "  R8CB_NP_FA claims the matrix is singular.\n" );
      printf ( "  The value of info is %d\n", info );
      return;
    }
/*
  Solve the system.
*/
    free ( x );

    x = r8cb_np_sl ( N, ML, MU, a, b, job );

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

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

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

void r8cb_np_sl_test ( )

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

    R8CB_NP_SL_TEST tests R8CB_NP_SL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 March 2013

  Author:

    John Burkardt
*/
{
# define N 10
# define ML 1
# define MU 2

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

  printf ( "\n" );
  printf ( "R8CB_NP_SL_TEST\n" );
  printf ( "  R8CB_NP_SL solves a linear system factored by R8CB_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N     = %d\n", N );
  printf ( "  Lower bandwidth ML = %d\n", ML );
  printf ( "  Upper bandwidth MU = %d\n", MU );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    a = r8cb_random ( N, N, ML, MU, &seed );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the right hand side.
*/
    if ( job == 0 )
    {
      b = r8cb_mv ( N, ML, MU, a, x );
    }
    else
    {
      b = r8cb_mtv ( N, ML, MU, a, x );
    }
/*
  Factor the matrix.
*/
    info = r8cb_np_fa ( N, ML, MU, a );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R8CB_NP_SL_TEST - Fatal error!\n" );
      printf ( "  R8CB_NP_FA claims the matrix is singular.\n" );
      printf ( "  The value of info is %d\n", info );
      return;
    }
/*
  Solve the system.
*/
    free ( x );

    x = r8cb_np_sl ( N, ML, MU, a, b, job );

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

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

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

void r8cb_print_test ( )

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

    R8CB_PRINT_TEST tests R8CB_PRINT.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    21 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8CB_PRINT_TEST\n" );
  printf ( "  R8CB_PRINT prints an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB matrix:" );

  free ( a );

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

void r8cb_print_some_test ( )

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

    R8CB_PRINT_SOME_TEST tests R8CB_PRINT_SOME.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8CB_PRINT_SOME_TEST\n" );
  printf ( "  R8CB_PRINT_SOME prints some of an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print_some ( m, n, ml, mu, a, 2, 2, 5, 5, "  Rows 2-5, Cols 2-5:" );

  free ( a );

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

void r8cb_random_test ( )

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

    R8CB_RANDOM_TEST tests R8CB_RANDOM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    27 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;
  int seed = 123456789;

  printf ( "\n" );
  printf ( "R8CB_RANDOM_TEST\n" );
  printf ( "  R8CB_RANDOM randomizes an R8CB matrix;\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_random ( m, n, ml, mu, &seed );

  r8cb_print ( m, n, ml, mu, a, "  The random R8CB matrix:" );

  free ( a );

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

void r8cb_to_r8ge_test ( )

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

    R8CB_TO_R8GE_TEST tests R8CB_TO_R8GE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *a_r8ge;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8CB_TO_R8GE_TEST\n" );
  printf ( "  R8CB_TO_R8GE converts an R8CB matrix to R8GE format.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB matrix:" );

  a_r8ge = r8cb_to_r8ge ( m, n, ml, mu, a );

  r8ge_print ( m, n, a_r8ge, "  The R8GE matrix:" );

  free ( a );
  free ( a_r8ge );

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

void r8cb_to_r8vec_test ( )

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

    R8CB_TO_R8VEC_TEST tests R8CB_TO_R8VEC.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *a_r8vec;
  double *a2;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8CB_TO_R8VEC_TEST\n" );
  printf ( "  R8CB_TO_R8VEC converts an R8CB matrix to an R8VEC.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB matrix:" );

  a_r8vec = r8cb_to_r8vec ( m, n, ml, mu, a );

  r8vec_print ( m * n, a_r8vec, "  The R8VEC:" );

  a2 = r8vec_to_r8cb ( m, n, ml, mu, a_r8vec );

  r8cb_print ( m, n, ml, mu, a2, "  The recovered R8CB matrix:" );

  free ( a );
  free ( a_r8vec );
  free ( a2 );

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

void r8cb_zeros_test ( )

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

    R8CB_ZEROS_TEST tests R8CB_ZEROS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    10 July 2016

  Author:

    John Burkardt
*/
{
# define M 8
# define N 10
# define ML 2
# define MU 3

  double *a;

  printf ( "\n" );
  printf ( "R8CB_ZEROS_TEST\n" );
  printf ( "  R8CB_ZEROS zeros an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );
  printf ( "  Lower bandwidth ML = %d\n", ML );
  printf ( "  Upper bandwidth MU = %d\n", MU );

  a = r8cb_zeros ( N, ML, MU );

  r8cb_print ( M, N, ML, MU, a, "  The R8CB zero matrix:" );

  free ( a );

  return;
# undef M
# undef N
# undef ML
# undef MU
}
/******************************************************************************/

void r8vec_to_r8cb_test ( )

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

    R8VEC_TO_R8CB_TEST tests R8VEC_TO_R8CB.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    26 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *a_r8vec;
  double *a2;
  int m = 8;
  int ml = 2;
  int mu = 3;
  int n = 10;

  printf ( "\n" );
  printf ( "R8VEC_TO_R8CB_TEST\n" );
  printf ( "  R8VEC_TO_R8CB converts an R8VEC to an R8CB matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", m );
  printf ( "  Matrix columns N = %d\n", n );
  printf ( "  Lower bandwidth ML = %d\n", ml );
  printf ( "  Upper bandwidth MU = %d\n", mu );

  a = r8cb_indicator ( m, n, ml, mu );

  r8cb_print ( m, n, ml, mu, a, "  The R8CB matrix:" );

  a_r8vec = r8cb_to_r8vec ( m, n, ml, mu, a );

  r8vec_print ( m * n, a_r8vec, "  The R8VEC:" );

  a2 = r8vec_to_r8cb ( m, n, ml, mu, a_r8vec );

  r8cb_print ( m, n, ml, mu, a2, "  The recovered R8CB matrix:" );

  free ( a );
  free ( a_r8vec );
  free ( a2 );

  return;
}