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

# include "r8sm.h"

int main ( );
void r8sm_indicator_test ( );
void r8sm_ml_test ( );
void r8sm_mtv_test ( );
void r8sm_mv_test ( );
void r8sm_print_test ( );
void r8sm_print_some_test ( );
void r8sm_random_test ( );
void r8sm_sl_test ( );
void r8sm_to_r8ge_test ( );
void r8sm_zeros_test ( );

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

int main ( )

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

    MAIN is the main program for R8SM_PRB.

  Discussion:

    R8SM_PRB tests R8SM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    01 June 2016

  Author:

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

  r8sm_indicator_test ( );
  r8sm_ml_test ( );
  r8sm_mtv_test ( );
  r8sm_mv_test ( );
  r8sm_print_test ( );
  r8sm_print_some_test ( );
  r8sm_random_test ( );
  r8sm_sl_test ( );
  r8sm_to_r8ge_test ( );
  r8sm_zeros_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "R8SM_PRB\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r8sm_indicator_test ( )

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

    R8SM_INDICATOR_TEST tests R8SM_INDICATOR.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 5;
  int n = 4;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_INDICATOR_TEST\n" );
  printf ( "  R8SM_INDICATOR sets up an R8SM indicator matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  r8sm_indicator ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM indicator matrix:" );

  free ( a );
  free ( u );
  free ( v );

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

void r8sm_ml_test ( )

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

    R8SM_ML_TEST tests R8SM_ML.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    20 April 2013

  Author:

    John Burkardt
*/
{
# define M 7
# define N 7

  double a[M*N];
  double *b;
  double *b2;
  int info;
  int i;
  int job;
  int pivot[N];
  int seed = 123456789;
  double u[M];
  double v[N];
  double *x;

  printf ( "\n" );
  printf ( "R8SM_ML_TEST\n" );
  printf ( "  R8SM_ML computes A*x or A'*X\n" );
  printf ( "  where A is a Sherman Morrison matrix.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );

  for ( job = 0; job <= 1; job++ )
  {
/*
  Set the matrix.
*/
    r8sm_random ( M, N, &seed, a, u, v );

    r8sm_print ( M, N, a, u, v, "  The Sherman Morrison matrix:" );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    if ( job == 0 )
    {
      b = r8sm_mv ( M, N, a, u, v, x );
    }
    else
    {
      b = r8sm_mtv ( M, N, a, u, v, x );
    }
/*
  Factor the matrix.
*/
    info = r8ge_fa ( N, a, pivot );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "  Fatal error!\n" );
      printf ( "  R8GE_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 = r8sm_ml ( N, a, u, v, pivot, 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 ( b );
    free ( b2 );
    free ( x );
  }

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

void r8sm_mtv_test ( )

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

    R8SM_MTV_TEST tests R8SM_MTV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int m = 5;
  int n = 4;
  double *u;
  double *v;
  double *x;

  printf ( "\n" );
  printf ( "R8SM_MTV_TEST\n" );
  printf ( "  R8SM_MTV computes b=A'*x, where A is an R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );
  
  r8sm_indicator ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM matrix:" );

  x = r8vec_indicator1_new ( m );

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

  b = r8sm_mtv ( m, n, a, u, v, x );

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

  free ( a );
  free ( b );
  free ( u );
  free ( v );
  free ( x );

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

void r8sm_mv_test ( )

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

    R8SM_MV_TEST tests R8SM_MV.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *b;
  int m = 5;
  int n = 4;
  double *u;
  double *v;
  double *x;

  printf ( "\n" );
  printf ( "R8SM_MV_TEST\n" );
  printf ( "  R8SM_MV computes b=A*x, where A is an R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );
  
  r8sm_indicator ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM matrix:" );

  x = r8vec_indicator1_new ( n );

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

  b = r8sm_mv ( m, n, a, u, v, x );

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

  free ( a );
  free ( b );
  free ( u );
  free ( v );
  free ( x );

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

void r8sm_print_test ( )

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

    R8SM_PRINT_TEST tests R8SM_PRINT.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 5;
  int n = 4;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_PRINT_TEST\n" );
  printf ( "  R8SM_PRINT prints an R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  r8sm_indicator ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM matrix:" );

  free ( a );
  free ( u );
  free ( v );

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

void r8sm_print_some_test ( )

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

    R8SM_PRINT_SOME_TEST tests R8SM_PRINT_SOME.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    02 June 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 9;
  int n = 9;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_PRINT_SOME_TEST\n" );
  printf ( "  R8SM_PRINT_SOME prints some of an R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  r8sm_indicator ( m, n, a, u, v );

  r8sm_print_some ( m, n, a, u, v, 2, 3, 5, 7, "  Rows 2-5, Cols 3-7:" );

  free ( a );
  free ( u );
  free ( v );

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

void r8sm_random_test ( )

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

    R8SM_RANDOM_TEST tests R8SM_RANDOM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 5;
  int n = 4;
  int seed;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_RANDOM_TEST\n" );
  printf ( "  R8SM_RANDOM sets up a random R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  seed = 123456789;
  r8sm_random ( m, n, &seed, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The random R8SM matrix:" );

  free ( a );
  free ( u );
  free ( v );

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

void r8sm_sl_test ( )

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

    R8SM_SL_TEST tests R8SM_SL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    19 April 2013

  Author:

    John Burkardt
*/
{
# define M 5
# define N 5

  double a[M*N];
  double *b;
  int i;
  int ierror;
  int info;
  int job;
  int pivot[N];
  int seed = 123456789;
  double u[M];
  double v[N];
  double *x;

  printf ( "\n" );
  printf ( "R8SM_SL_TEST\n" );
  printf ( "  R8SM_SL implements the Sherman-Morrison method \n" );
  printf ( "  for solving a perturbed linear system.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );

  for ( job = 1; 0 <= job; job-- )
  {
/*
  Set the matrix.
*/
    r8sm_random ( M, N, &seed, a, u, v );

    r8sm_print ( M, N, a, u, v, "  The Sherman-Morrison matrix A:" );
/*
  Set the desired solution.
*/
    x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
    if ( job == 0 )
    {
      b = r8sm_mv ( M, N, a, u, v, x );
    }
    else
    {
      b = r8sm_mtv ( M, N, a, u, v, x );
    }

    r8vec_print ( N, b, "  The right hand side vector B:" );
/*
  Factor the matrix.
*/
    info = r8ge_fa ( N, a, pivot );

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "  Fatal error!\n" );
      printf ( "  R8GE_FA declares the matrix is singular!\n" );
      printf ( "  The value of INFO is %d\n", info );
      continue;
    }
/*
  Solve the linear system.
*/
    b = r8sm_sl ( N, a, u, v, b, pivot, job );
 
    if ( job == 0 )
    {
      r8vec_print ( N, b, "  Solution to A * X = B:" );
    }
    else
    {
      r8vec_print ( N, b, "  Solution to A' * X = B:" );
    }

    free ( b );
  }

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

void r8sm_to_r8ge_test ( )

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

    R8SM_TO_R8GE_TEST tests R8SM_TO_R8GE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  double *a_r8ge;
  int m = 5;
  int n = 4;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_TO_R8GE_TEST\n" );
  printf ( "  R8SM_TO_R8GE converts an R8SM matrix to R8GE format;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  r8sm_indicator ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM matrix:" );

  a_r8ge = r8sm_to_r8ge ( m, n, a, u, v );

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

  free ( a );
  free ( a_r8ge );
  free ( u );
  free ( v );

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

void r8sm_zeros_test ( )

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

    R8SM_ZEROS_TEST tests R8SM_ZEROS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    30 May 2016

  Author:

    John Burkardt
*/
{
  double *a;
  int m = 5;
  int n = 4;
  double *u;
  double *v;

  printf ( "\n" );
  printf ( "R8SM_ZEROS_TEST\n" );
  printf ( "  R8SM_ZEROS sets up a zero R8SM matrix;\n" );
  printf ( "\n" );
  printf ( "  M = %d\n", m );
  printf ( "  N = %d\n", n );

  a = ( double * ) malloc ( m * n * sizeof ( double ) );
  u = ( double * ) malloc ( m * sizeof ( double ) );
  v = ( double * ) malloc ( n * sizeof ( double ) );

  r8sm_zeros ( m, n, a, u, v );

  r8sm_print ( m, n, a, u, v, "  The R8SM zero matrix:" );

  free ( a );
  free ( u );
  free ( v );

  return;
}