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

# include "r8ge_np.h"

int main ( );
void r8ge_np_det_test ( );
void r8ge_np_fa_test ( );
void r8ge_np_inverse_test ( );
void r8ge_np_ml_test ( );
void r8ge_np_sl_test ( );
void r8ge_np_trf_test ( );
void r8ge_np_trm_test ( );
void r8ge_np_trs_test ( );

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

int main ( )

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

    MAIN is the main program for R8GE_NP_TEST.

  Discussion:

    R8GE_NP_TEST tests R8GE_NP.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    09 June 2016

  Author:

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

  r8ge_np_det_test ( );
  r8ge_np_fa_test ( );
  r8ge_np_inverse_test ( );
  r8ge_np_ml_test ( );
  r8ge_np_sl_test ( );
  r8ge_np_trf_test ( );
  r8ge_np_trm_test ( );
  r8ge_np_trs_test ( );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "R8GE_NP_TEST\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void r8ge_np_det_test ( )

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

    R8GE_NP_DET_TEST tests R8GE_NP_DET.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 June 2016

  Author:

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

  printf ( "\n" );
  printf ( "R8GE_NP_DET_TEST\n" );
  printf ( "  R8GE_NP_DET computes the determinant of a matrix\n" );
  printf ( "  that was factored by R8GE_NP_FA,\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", n );
/*
  Set the matrix.
*/
  a = r8ge_dif2 ( n, n );
/*
  Factor the matrix.
*/
  info = r8ge_np_fa ( n, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_DET_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_FA declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Get the determinant.
*/
  det = r8ge_np_det ( n, a );
  printf ( "\n" );
  printf ( "  Determinant of -1, 2, -1 matrix is %g\n", det );
  printf ( "  Exact value is %g\n", ( double ) ( n + 1 ) );

  free ( a );

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

void r8ge_np_fa_test ( )

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

    R8GE_NP_FA_TEST tests R8GE_NP_FA.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 April 2013

  Author:

    John Burkardt
*/
{
# define N 10

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

  printf ( "\n" );
  printf ( "R8GE_NP_FA_TEST\n" );
  printf ( "  R8GE_NP_FA LU factors an R8GE matrix without pivoting,\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( N, N, &seed );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );

  for ( i = 0; i < N; i++ )
  {
    x[i] = 1.0;
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8ge_mv ( N, N, a, x );
/*
  Factor the matrix.
*/
  info = r8ge_np_fa ( N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_FA_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_FA declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Solve the linear system.
*/
  job = 0;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );
 
  r8vec_print_some ( N, x, 10, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 0;
  free ( b );
  b = r8ge_np_ml ( N, a, x, job );
/*
  Solve the system
*/
  job = 0;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );

  r8vec_print_some ( N, x, 10, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 1;
  free ( b );
  b = r8ge_np_ml ( N, a, x, job );
/*
  Solve the system
*/
  job = 1;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );

  r8vec_print ( N, x, "  Solution of transposed system:" );

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

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

void r8ge_np_inverse_test ( )

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

    R8GE_NP_INVERSE_TEST tests R8GE_NP_INVERSE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    06 April 2013

  Author:

    John Burkardt
*/
{
# define N 5

  double *a;
  double *a_lu;
  double *b;
  double *c;
  int i;
  int info;
  int seed = 123456789;
  int j;

  printf ( "\n" );
  printf ( "R8GE_NP_INVERSE_TEST\n" );
  printf ( "  R8GE_NP_INVERSE computes the inverse of a matrix\n" );
  printf ( "  that was factored by R8GE_NP_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( N, N, &seed );

  r8ge_print ( N, N, a, "  The random matrix:" );
/*
  Factor and invert the matrix.
*/
  a_lu = ( double * ) malloc ( N * N * sizeof ( double ) );

  for ( i = 0; i < N; i++ )
  {
    for ( j = 0; j < N; j++ )
    {
      a_lu[i+j*N] = a[i+j*N];
    }
  }

  info = r8ge_np_fa ( N, a_lu );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_INVERSE_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_FA declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  b = r8ge_np_inverse ( N, a_lu );

  r8ge_print ( N, N, b, "  The inverse matrix:" );
/*
  Compute A * B = C.
*/
  c = r8ge_mm ( N, N, N, a, b );

  r8ge_print ( N, N, c, "  The product:" );

  free ( a );
  free ( a_lu );
  free ( b );
  free ( c );

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

void r8ge_np_ml_test ( )

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

    R8GE_NP_ML_TEST tests R8GE_NP_ML.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    09 April 2013

  Author:

    John Burkardt
*/
{
# define N 10

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

  printf ( "\n" );
  printf ( "R8GE_NP_ML_TEST\n" );
  printf ( "  For a matrix in general storage,\n" ); 
  printf ( "  R8GE_NP_ML computes A*x or A'*X\n" );
  printf ( "  where A has been factored by R8GE_NP_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );

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

    if ( info != 0 )
    {
      printf ( "\n" );
      printf ( "R8GE_NP_ML_TEST - Fatal error!\n" );
      printf ( "  R8GE_NP_FA declares the matrix is singular!\n" );
      printf ( "  The value of INFO is %d\n", info );
      continue;
    }
/*
  Now multiply factored matrix times solution to get right hand side again.
*/
    b2 = r8ge_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 r8ge_np_sl_test ( )

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

    R8GE_NP_SL_TEST tests R8GE_NP_FA, R8GE_NP_SL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 April 2013

  Author:

    John Burkardt
*/
{
# define N 10

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

  printf ( "\n" );
  printf ( "R8GE_NP_SL_TEST\n" );
  printf ( "  R8GE_NP_SL solves a system factored by R8GE_NP_FA.\n" );
  printf ( "\n" );
  printf ( "  Matrix order N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( N, N, &seed );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );

  for ( i = 0; i < N; i++ )
  {
    x[i] = 1.0;
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8ge_mv ( N, N, a, x );
/*
  Factor the matrix.
*/
  info = r8ge_np_fa ( N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_SL_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_FA declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Solve the linear system.
*/
  job = 0;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );
 
  r8vec_print_some ( N, x, 10, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 0;
  free ( b );
  b = r8ge_np_ml ( N, a, x, job );
/*
  Solve the system
*/
  job = 0;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );

  r8vec_print_some ( N, x, 10, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 1;
  free ( b );
  b = r8ge_np_ml ( N, a, x, job );
/*
  Solve the system
*/
  job = 1;
  free ( x );
  x = r8ge_np_sl ( N, a, b, job );

  r8vec_print ( N, x, "  Solution of transposed system:" );

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

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

void r8ge_np_trf_test ( )

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

    R8GE_NP_TRF_TEST tests R8GE_NP_TRF.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    12 April 2013

  Author:

    John Burkardt
*/
{
# define M 10
# define N 10
# define NRHS 1

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

  printf ( "\n" );
  printf ( "R8GE_NP_TRF_TEST\n" );
  printf ( "  R8GE_NP_TRF factors an R8GE matrix without pivoting,\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( M, N, &seed );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );
  for ( i = 0; i < N; i++ )
  {
    x[i] = 1.0;
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8ge_mv ( M, N, a, x );
/*
  Factor the matrix.
*/
  info = r8ge_np_trf ( M, N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRF_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRF declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Solve the linear system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );
 
  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRF_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 0;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system
*/
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRF_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 1;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'T', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRF_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution of transposed system:" );

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

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

void r8ge_np_trm_test ( )

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

    R8GE_NP_TRM_TEST tests R8GE_NP_TRM.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    12 April 2013

  Author:

    John Burkardt
*/
{
# define M 10
# define N 10
# define NRHS 1

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

  printf ( "\n" );
  printf ( "R8GE_NP_TRM_TEST\n" );
  printf ( "  R8GE_NP_TRM computes A*x when A was factored by R8GE_NP_TRF.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( M, N, &seed );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );
  for ( i = 0; i < N; i++ )
  {
    x[i] = 1.0;
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8ge_mv ( M, N, a, x );
/*
  Factor the matrix.
*/
  info = r8ge_np_trf ( M, N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRM_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRF declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Solve the linear system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );
 
  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRM_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 0;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system
*/
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRM_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 1;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'T', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRM_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution of transposed system:" );

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

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

void r8ge_np_trs_test ( )

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

    R8GE_NP_TRS_TEST tests R8GE_NP_TRS.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    12 April 2013

  Author:

    John Burkardt
*/
{
# define M 10
# define N 10
# define NRHS 1

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

  printf ( "\n" );
  printf ( "R8GE_NP_TRS_TEST\n" );
  printf ( "  R8GE_NP_TRS solves a linear system factored by R8GE_NP_TRF.\n" );
  printf ( "\n" );
  printf ( "  Matrix rows M =    %d\n", M );
  printf ( "  Matrix columns N = %d\n", N );
/*
  Set the matrix.
*/
  a = r8ge_random ( M, N, &seed );
/*
  Set the desired solution.
*/
  x = ( double * ) malloc ( N * sizeof ( double ) );
  for ( i = 0; i < N; i++ )
  {
    x[i] = 1.0;
  }
/*
  Compute the corresponding right hand side.
*/
  b = r8ge_mv ( M, N, a, x );
/*
  Factor the matrix.
*/
  info = r8ge_np_trf ( M, N, a );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRS_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRF declares the matrix is singular!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
/*
  Solve the linear system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );
 
  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRS_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }
  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 0;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system
*/
  x = r8ge_np_trs ( N, NRHS, 'N', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRS_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution:" );
/*
  Set the desired solution.
*/
  free ( x );
  x = r8vec_indicator1_new ( N );
/*
  Compute the corresponding right hand side.
*/
  job = 1;
  free ( b );
  b = r8ge_np_trm ( M, N, a, x, job );
/*
  Solve the system.
*/
  free ( x );
  x = r8ge_np_trs ( N, NRHS, 'T', a, b );

  if ( info != 0 )
  {
    printf ( "\n" );
    printf ( "R8GE_NP_TRS_TEST - Fatal error!\n" );
    printf ( "  R8GE_NP_TRS returned an error condition!\n" );
    printf ( "  The value of INFO is %d\n", info );
    return;
  }

  r8vec_print ( N, x, "  Solution of transposed system:" );

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

  return;
# undef M
# undef N
# undef NRHS
}