[Help-gsl] finite element assembly and calling UMFPACK for linear system solution

[Brijesh Upadhaya]

Hi everyone,

I am working on a finite element problem and wanted to ask if anyone of you
have tried using UMFPACK for the linear algebra. I have following findings
when using sparse matrix module from gsl-2.5.

1. In CCS the row indices are not sorted so UMFPACK gives error (see
attached example)
2. Adding boundary conditions were also difficult once compressed, so I
couldn't utilize gsl_spmatrix_add() to add mass matrix with stiffness
matrix. Instead I had to do it with triplet and calling the get/set
routines. Could only add boundary condition at triplet level and compressed
at the end and passed the pointers to UMFPACK routines.
3. Another thing (not related to gsl) was that UMFPACK don't use (size_t
*), instead it has (int *) and (long *) for column pointers and row
indices. I tried to cast (size_t *) to (long *) to utilize the pointers of
gsl_spmatrix object.

I am bit curious to know if anyone have a better experience. I wanted to
use direct solver instead of iterative solver.

King Regards,
Brijesh

/* 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;
}