/* test_lsolver.c
 *
 * Test UMFPACK (V5.7.1) example
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_spmatrix.h>
#include <umfpack.h>


static int
test_lumfpack (gsl_spmatrix *A, gsl_vector *X, const gsl_vector *F)
{
  const long int n = (const long int) X->size;
  const long int m = (const long int) X->size;
  const long int *Ap = (const long int *) A->p;
  const long int *Ai = (const long int *) A->i;
  const double *Ax = (const double *) A->data;

  double *x = X->data;
  const double *b = F->data;

  double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL];
  void *Symbolic, *Numeric;

  umfpack_dl_defaults (Control);

  int status = umfpack_dl_symbolic (n, m, Ap, Ai, Ax, &Symbolic, Control, Info);
  if (status)
    {
      umfpack_dl_free_symbolic (&Symbolic);
      umfpack_dl_report_info (Control, Info);
      umfpack_dl_report_status (Control, status);
      return status;
    }

  status = umfpack_dl_numeric (Ap, Ai, Ax, Symbolic, &Numeric, Control, Info);
  if (status)
    {
      umfpack_dl_free_symbolic (&Symbolic);
      umfpack_dl_free_numeric (&Numeric);
      umfpack_dl_report_info (Control, Info);
      umfpack_dl_report_status (Control, status);
      return status;
    }

  status = umfpack_dl_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, Control, Info);
  if (status)
    {
      umfpack_dl_free_symbolic (&Symbolic);
      umfpack_dl_free_numeric (&Numeric);
      umfpack_dl_report_info (Control, Info);
      umfpack_dl_report_status (Control, status);
      return status;
    }

  umfpack_dl_free_symbolic (&Symbolic);
  umfpack_dl_free_numeric (&Numeric);

  return status;
}


/*=============================================================================*/
int 
main (void)
{
  int status;
  size_t Ap[] = { 0, 2, 5, 9, 10, 12 };
  size_t Ai[] = { 0, 1, 0, 2, 4, 1, 2, 3, 4, 2, 1, 4 };
  double Ax[] = { 2., 3., 3., -1., 4., 4., -3., 1., 2., 2., 6., 1. };
  double b[] = { 8., 45., -3., 3., 19. };
  double u[] = { [0] = 0.0, [4] = 0.0 };
  const size_t n = sizeof (b) / sizeof (double);
  size_t j, p;

  if (1)
    {
      // In practical finite element assembly this happens 
      // Swap row indices and corrosponding data 
      gsl_vector_ulong_view vi = gsl_vector_ulong_view_array (Ai, n);
      gsl_vector_ulong_swap_elements (&vi.vector, 2, 3);
      cblas_dswap (1, &Ax[2], 1, &Ax[3], 1); 
    }
  
  gsl_spmatrix *K = gsl_spmatrix_alloc (n, n);

  for (j = 0; j < n; j++)
    for (p = Ap[j]; p < Ap[j + 1]; p++)
      gsl_spmatrix_set (K, Ai[p], j, Ax[p]);
  
  gsl_spmatrix *A = gsl_spmatrix_ccs (K);

  gsl_vector_view x = gsl_vector_view_array (u, n);
  gsl_vector_const_view f = gsl_vector_const_view_array (b, n);
  
  // Solve system
  status = test_lumfpack (A, &x.vector, &f.vector);
  
  if (!status)
    gsl_vector_fprintf (stdout, &x.vector, "%f");

  if (0)
    gsl_spmatrix_fprintf (stderr, K, "%f");
  else if (1)
    gsl_spmatrix_fprintf (stderr, A, "%f"); 

  gsl_spmatrix_free (K);
  gsl_spmatrix_free (A);

  return 0;
}