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

# include "r8sd.h"

int main ( );
void r8sd_cg_test ( );
void r8sd_dif2_test ( );
void r8sd_indicator_test ( );
void r8sd_mv_test ( );
void r8sd_print_test ( );
void r8sd_print_some_test ( );
void r8sd_random_test ( );
void r8sd_res_test ( );
void r8sd_to_r8ge_test ( );
void r8sd_zeros_test ( );

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

int main ( )

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

    MAIN is the main program for R8SD_TEST.

  Discussion:

    R8SD_TEST tests R8SD.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 July 2016

  Author:

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

  r8sd_cg_test ( );
  r8sd_dif2_test ( );
  r8sd_indicator_test ( );
  r8sd_mv_test ( );
  r8sd_print_test ( );
  r8sd_print_some_test ( );
  r8sd_random_test ( );
  r8sd_res_test ( );
  r8sd_to_r8ge_test ( );
  r8sd_zeros_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "R8SD_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r8sd_cg_test ( )

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

    R8SD_CG_TEST tests R8SD_CG.

  Discussion:

    NX and NY are the number of grid points in the X and Y directions.
    N is the number of unknowns.
    NDIAG is the number of nonzero diagonals we will store.  We only
    store the main diagonal, and the superdiagonals.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    18 April 2013

  Author:

    John Burkardt
*/
{
# define NDIAG 3
# define NX 10
# define NY 10
# define N ( NX * NY )

  double a[N*NDIAG];
  double *b;
  double *b2;
  double err;
  int i;
  int j;
  int k;
  int offset[NDIAG] = { 0, 1, NX };
  double *x;
  double x_init[N];

  printf ( "\n" );
  printf ( "R8SD_CG_TEST\n" );
  printf ( "  R8SD_CG applies the conjugate gradient method\n" );
  printf ( "  to a symmetric positive definite linear\n" );
  printf ( "  system stored by diagonals.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
  printf ( "  Matrix diagonals NDIAG = %d\n", NDIAG );
  printf ( "\n" );
/*
  OFFSET tells us where the nonzero diagonals are.  It does this
  by recording how "high" or to the right the diagonals are from
  the main diagonal.

  Now we compute the numbers that go into the diagonals.  For this
  problem, we could simply store a column of 4's, and two columns of
  -1's, but I wanted to go through the motions of thinking about the
  value of each entry.  "K" counts the row of the original matrix
  that we are working on.
*/
  k = 0;

  for ( j = 1; j <= NY; j++ )
  {
    for ( i = 1; i <= NX; i++ )
    {
/*
  Central
*/
      a[k+0*N] = 4.0;
/*
  East ( = West )
*/
      if ( i == NX )
      {
        a[k+1*N] = 0.0;
      }
      else
      {
        a[k+1*N] = -1.0;
      }
/*
  North ( = South )
*/
      if ( j == NY )
      {
        a[k+2*N] = 0.0;
      }
      else
      {
        a[k+2*N] = -1.0;
      }
      k = k + 1;
    }
  }
/*
  Print some of the matrix.
*/
  r8sd_print_some ( N, NDIAG, offset, a, 1, 1, 10, 10, 
    "  First 10 rows and columns of matrix." );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );
  k = 0;
  for ( j = 1; j <= NY; j++ )
  {
    for ( i = 1; i <= NX; i++ )
    {
      x[k] = ( double ) ( 10 * i + j );
      k = k + 1;
    }
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8sd_mv ( N, N, NDIAG, offset, a, x );

  r8vec_print_some ( N, b, 10, "  Right hand side:" );
/*
  Set X to zero so no one accuses us of cheating.
*/
  for ( i = 0; i < N; i++ )
  {
    x_init[i] = 0.0;
  }
/*
  Call the conjugate gradient method.
*/
  free ( x );
  x = r8sd_cg ( N, NDIAG, offset, a, b, x_init );
/*
  Compute the residual, A*x-b
*/
  b2 = r8sd_mv ( N, N, NDIAG, offset, a, x );
 
  err = 0.0;
  for ( i = 0; i < N; i++ )
  {
    err = r8_max ( err, fabs ( b2[i] - b[i] ) );
  }
 
  r8vec_print_some ( N, x, 10, "  Solution:" );

  printf ( "\n" );
  printf ( "  Maximum residual = %g\n", err );
/*
  Note that if we're not satisfied with the solution, we can
  call again, using the computed X as our starting estimate.

  Call the conjugate gradient method AGAIN.
*/
  for ( i = 0; i < N; i++ )
  {
    x_init[i] = x[i];
  }
  free ( x );

  x = r8sd_cg ( N, NDIAG, offset, a, b, x_init );
/*
  Compute the residual, A*x-b
*/
  free ( b2 );
  b2 = r8sd_mv ( N, N, NDIAG, offset, a, x );
 
  err = 0.0;
  for ( i = 0; i < N; i++ )
  {
    err = r8_max ( err, fabs ( b2[i] - b[i] ) );
  }
 
  r8vec_print_some ( N, x, 10, "  Second attempt at solution:" );

  printf ( "\n" );
  printf ( "  Maximum residual of second attempt = %g\n", err );

  free ( b );
  free ( b2 );
  free ( x );

  return;
# undef N
# undef NDIAG
# undef NX
# undef NY
}
/******************************************************************************/

void r8sd_dif2_test ( )

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

    R8SD_DIF2_TEST tests R8SD_DIF2.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int ndiag = 2;
  int offset[2] = { 0, 1 };

  printf ( "\n" );
  printf ( "R8SD_DIF2_TEST\n" );
  printf ( "  R8SD_DIF2 sets up an R8SD second difference matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_dif2 ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD second difference matrix:" );

  free ( a );

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

void r8sd_indicator_test ( )

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

    R8SD_INDICATOR_TEST tests R8SD_INDICATOR.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    20 April 2013

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };

  printf ( "\n" );
  printf ( "R8SD_INDICATOR_TEST\n" );
  printf ( "  R8SD_INDICATOR sets up an R8SD indicator matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_indicator ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD indicator matrix:" );

  free ( a );

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

void r8sd_mv_test ( )

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

    R8SD_MV_TEST tests R8SD_MV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };
  double *x;

  printf ( "\n" );
  printf ( "R8SD_MV_TEST\n" );
  printf ( "  R8SD_MV computes b=A*x, where A is an R8SD matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_indicator ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD indicator matrix:" );

  x = r8vec_indicator1_new ( n );

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

  b = r8sd_mv ( n, n, ndiag, offset, a, x );

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

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

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

void r8sd_print_test ( )

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

    R8SD_PRINT_TEST tests R8SD_PRINT.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };

  printf ( "\n" );
  printf ( "R8SD_PRINT_TEST\n" );
  printf ( "  R8SD_PRINT prints an R8SD matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_indicator ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD matrix:" );

  free ( a );

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

void r8sd_print_some_test ( )

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

    R8SD_PRINT_SOME_TEST tests R8SD_PRINT_SOME.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 10;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };

  printf ( "\n" );
  printf ( "R8SD_PRINT_SOME_TEST\n" );
  printf ( "  R8SD_PRINT prints some of an R8SD matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_indicator ( n, ndiag, offset );

  r8sd_print_some ( n, ndiag, offset, a, 1, 2, 7, 6, "  Rows 1-7, Cols 2-6:" );

  free ( a );

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

void r8sd_random_test ( )

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

    R8SD_RANDOM_TEST tests R8SD_RANDOM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    20 April 2013

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };
  int seed;

  printf ( "\n" );
  printf ( "R8SD_RANDOM_TEST\n" );
  printf ( "  R8SD_RANDOM randomizes an R8SD matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  seed = 123456789;
  a = r8sd_random ( n, ndiag, offset, &seed );

  r8sd_print ( n, ndiag, offset, a, "  The random R8SD matrix:" );

  free ( a );

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

void r8sd_res_test ( )

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

    R8SD_RES_TEST tests R8SD_RES.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    17 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  double e;
  int i;
  int n = 10;
  int ndiag = 2;
  int offset[2] = { 0, 1 };
  double *r;
  int seed;
  double *x;
  double *x2;

  printf ( "\n" );
  printf ( "R8SD_RES_TEST\n" );
  printf ( "  R8SD_RES computes a residual R=b-A*x\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Number of diagonals is %d\n", ndiag );

  seed = 123456789;
  a = r8sd_random ( n, ndiag, offset, &seed );
/*
  Print some of the matrix.
*/
  r8sd_print ( n, ndiag, offset, a, "  The R8SD matrix:" );
  x = r8vec_indicator1_new ( n );

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

  b = r8sd_mv ( n, n, ndiag, offset, a, x );

  r8vec_print ( n, b, "  The product b=A*x" );
/*
  Make X2, a bad copy of X.
*/
  x2 = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    e = r8_uniform_01 ( &seed );
    x2[i] = x[i] + 0.1 * e;
  }
  r8vec_print ( n, x2, "  The defective vector X2:" );
/*
  Compute R = B-A*X2.
*/
  r = r8sd_res ( n, n, ndiag, offset, a, x2, b );
  
  r8vec_print ( n, r, "  Residual b-A*x2:" );

  free ( a );
  free ( b );
  free ( r );
  free ( x );
  free ( x2 );

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

void r8sd_to_r8ge_test ( )

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

    R8SD_TO_R8GE_TEST tests R8SD_TO_R8GE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *a_r8ge;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };

  printf ( "\n" );
  printf ( "R8SD_TO_R8GE_TEST\n" );
  printf ( "  R8SD_TO_R8GE converts an R8SD matrix to R8GE format.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_indicator ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD matrix:" );

  a_r8ge = r8sd_to_r8ge ( n, ndiag, offset, a );

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

  free ( a );
  free ( a_r8ge );

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

void r8sd_zeros_test ( )

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

    R8SD_ZEROS_TEST tests R8SD_ZEROS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 July 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int n = 5;
  int ndiag = 3;
  int offset[3] = { 0, 1, 3 };

  printf ( "\n" );
  printf ( "R8SD_ZEROS_TEST\n" );
  printf ( "  R8SD_ZEROS zeros an R8SD matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
  printf ( "  Matrix diagonals NDIAG = %d\n", ndiag );

  a = r8sd_zeros ( n, ndiag, offset );

  r8sd_print ( n, ndiag, offset, a, "  The R8SD matrix:" );

  free ( a );

  return;
}