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

# include "ppc_array.h"
# include "ppc_sparse_matrix.h"

# include <suitesparse/umfpack.h>

int main ( );
void error_and_exit ( int status, const char *file, int line );
void timestamp ( );

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

int main ( )

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

    umfpack_test2() tests umfpack().

  Discussion:

    Replace matrix of umfpack_test1() by a singular matrix and test error check.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    10 September 2023

  Author:

    John Burkardt
*/
{
  double **a;
  int *Ai;
  int *Ap;
  double *Ax;
  double *b;
  int i;
  int j;
  int n = 5;
  void *Numeric;
  int nz;
  int status;
  void *Symbolic;
  double *x;

  printf ( "\n" );
  printf ( "umfpack_test2():\n" );
  printf ( "  C version\n" );
  printf ( "  Test umfpack().\n" );
/*
  Define the sparse matrix, right hand side, and solution.
*/
  make_matrix ( a, n, n );
  make_vector ( b, n );
  make_vector ( x, n );

  a[0][0] = 2;
  a[0][1] = 3;
  a[0][2] = 0;
  a[0][3] = 0;
  a[0][4] = 0;

  a[1][0] = 3;
  a[1][1] = 0;
  a[1][2] = 4;
  a[1][3] = 0;
  a[1][4] = 6;

  a[2][0] = 0;
  a[2][1] = -1;
  a[2][2] = -3;
  a[2][3] = 2;
  a[2][4] = 0;

  a[3][0] = 0;
  a[3][1] = 0;
/*
  a[3][2] = 1;
*/
  a[3][2] = 0;
  a[3][3] = 0;
  a[3][4] = 0;

  a[4][0] = 0;
  a[4][1] = 4;
  a[4][2] = 2;
  a[4][3] = 0;
  a[4][4] = 1;

  printf ( "\n" );
  printf ( "  The matrix a:\n" );
  printf ( "\n" );
  print_matrix ( " %f", a, n, n );

  b[0] = 8;
  b[1] = 45;
  b[2] = -3;
  b[3] = 3;
  b[4] = 19;

  printf ( "\n" );
  printf ( "  The right hand side b:\n" );
  printf ( "\n" );
  print_vector ( " %f", b, n );
/*
  Count nonzeros.
*/
  nz = 0;
  for ( i = 0; i < n; i++ )
  {
    for ( j = 0; j < n; j++ )
    {
      if ( a[i][j] != 0.0 )
      {
        nz = nz + 1;
      }
    }
  }
  printf ( "  Number of nonzero matrix entries = %d\n", nz );
/*
  Allocate sparse arrays.
*/
  make_vector ( Ai, nz );
  make_vector ( Ap, n + 1 );
  make_vector ( Ax, nz );
/*
  Convert matrix to sparse form.
*/
  ccs_pack ( a, n, n, Ap, Ai, Ax );
/*
  Print sparse data.
*/
  printf ( "\n" );
  printf ( "  Ai:" );
  printf ( "\n" );
  print_vector ( " %d", Ai, nz );

  printf ( "\n" );
  printf ( "  Ap:" );
  printf ( "\n" );
  print_vector ( " %d", Ap, n + 1 );

  printf ( "\n" );
  printf ( "  Ax:" );
  printf ( "\n" );
  print_vector ( " %g", Ax, nz );
/*
  Call umfpack.
*/
  status = umfpack_di_symbolic ( n, n, Ap, Ai, Ax, &Symbolic, NULL, NULL );
  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }

  status = umfpack_di_numeric ( Ap, Ai, Ax, Symbolic, &Numeric, NULL, NULL );
  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }

  status = umfpack_di_solve ( UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, NULL, NULL );
  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }

  printf ( "\n" );
  printf ( "  Computed solution x:\n" );
  printf ( "\n" );
  print_vector ( " %g", x, n );

  umfpack_di_free_symbolic ( &Symbolic );

  umfpack_di_free_numeric ( &Numeric );
/*
  Free memory.
*/
  free_matrix ( a );
  free_vector ( b );
  free_vector ( Ai );
  free_vector ( Ap );
  free_vector ( Ax );
  free_vector ( x );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "umfpack_test2():\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

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

void error_and_exit ( int status, const char *file, int line )

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

    error_and_exit() reports an error and exits.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2014

  Author:

    Rouben Rostamian
*/
{
  fprintf ( stderr, "file %s, line %d: umfpack_di_symbolic() failed.\n",
    file, line );

  switch ( status )
  {
    case UMFPACK_ERROR_out_of_memory:
      fprintf ( stderr, "  Out of memory.\n" );
      break;
    case UMFPACK_WARNING_singular_matrix:
      fprintf ( stderr, "The matrix is singular.\n" );
    default:
      fprintf ( stderr, "  UMFPACK error code %d\n", status );
  }
  exit ( EXIT_FAILURE );
}
/******************************************************************************/

void timestamp ( )

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

    timestamp() prints the current YMDHMS date as a time stamp.

  Example:

    17 June 2014 09:45:54 AM

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    17 June 2014

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  printf ( "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}