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

# include "ppc_array.h"

# include <suitesparse/umfpack.h>

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

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

int main ( )

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

    umfpack_test3() tests umfpack().

  Discussion:

    Use sparse triplet format for initial data.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 September 2023

  Author:

    John Burkardt
*/
{
  int *Ai;
  int *Ap;
  double *Ax;
  double *b;
  int n = 5;
  void *Numeric;
  int nz;
  int status;
  void *Symbolic;
  int *Ti;
  int *Tj;
  double *Tx;
  double *x;

  printf ( "\n" );
  printf ( "umfpack_test3():\n" );
  printf ( "  C version\n" );
  printf ( "  Test umfpack().\n" );
/*
  Define the sparse matrix, right hand side, and solution.
*/
  nz = 12;
  make_vector ( Ti, nz );
  make_vector ( Tj, nz );
  make_vector ( Tx, nz );

  nz = 0;
  Ti[nz] = 0;
  Tj[nz] = 0;
  Tx[nz] = 2.0;
  nz = nz + 1;

  Ti[nz] = 0;
  Tj[nz] = 1;
  Tx[nz] = 3.0;
  nz = nz + 1;

  Ti[nz] = 1;
  Tj[nz] = 0;
  Tx[nz] = 3.0;
  nz = nz + 1;

  Ti[nz] = 1;
  Tj[nz] = 2;
  Tx[nz] = 4.0;
  nz = nz + 1;

  Ti[nz] = 1;
  Tj[nz] = 4;
  Tx[nz] = 6.0;
  nz = nz + 1;

  Ti[nz] = 2;
  Tj[nz] = 1;
  Tx[nz] = -1.0;
  nz = nz + 1;

  Ti[nz] = 2;
  Tj[nz] = 2;
  Tx[nz] = -3.0;
  nz = nz + 1;

  Ti[nz] = 2;
  Tj[nz] = 3;
  Tx[nz] = 2.0;
  nz = nz + 1;

  Ti[nz] = 3;
  Tj[nz] = 2;
  Tx[nz] = 1.0;
  nz = nz + 1;

  Ti[nz] = 4;
  Tj[nz] = 1;
  Tx[nz] = 4.0;
  nz = nz + 1;

  Ti[nz] = 4;
  Tj[nz] = 2;
  Tx[nz] = 2.0;
  nz = nz + 1;

  Ti[nz] = 4;
  Tj[nz] = 4;
  Tx[nz] = 1.0;
  nz = nz + 1;

  printf ( "  Number of nonzero matrix entries = %d\n", nz );

  printf ( "\n" );
  printf ( "  The vector Ti:\n" );
  printf ( "\n" );
  print_vector ( " %d", Ti, nz );

  printf ( "\n" );
  printf ( "  The vector Tj:\n" );
  printf ( "\n" );
  print_vector ( " %d", Tj, nz );

  printf ( "\n" );
  printf ( "  The vector Tx:\n" );
  printf ( "\n" );
  print_vector ( " %f", Tx, nz );

  make_vector ( b, 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 );
/*
  Allocate sparse arrays.
*/
  make_vector ( Ai, nz );
  make_vector ( Ap, n + 1 );
  make_vector ( Ax, nz );
/*
  Convert matrix to sparse form.
*/
  status = umfpack_di_triplet_to_col ( n, n, nz, Ti, Tj, Tx, Ap, Ai, Ax, NULL );
  if ( status != UMFPACK_OK )
  {
    error_and_exit ( status, __FILE__, __LINE__ );
  }
/*
  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__ );
  }

  make_vector ( x, n );
  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_vector ( Ai );
  free_vector ( Ap );
  free_vector ( Ax );
  free_vector ( Ti );
  free_vector ( Tj );
  free_vector ( Tx );
  free_vector ( x );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "umfpack_test3():\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
}