/*
    Interface Routines nislvi, nifact and nisolv for NLEQ1S to call the 
    UMFPACK sparse linear solver. For a parameters description of these 
    routines see the file nleq1s.f.
    
    Written by:     L. Weimann
    Latest change:  March 2009

*/
#include <stdlib.h>
#include <stdio.h>
#include "umfpack.h"

int *Ap=NULL, *Ai=NULL;
double *x=NULL, *Aa=NULL;
int sys=UMFPACK_A;
double Control[UMFPACK_CONTROL], Info[UMFPACK_INFO];
void *Symbolic, *Numeric ;  

void nislvi_(int *n, int *nfmax, int *iopt, int *ifail, int *liwk, int *iwk,
             int *lpiwk, int *ltiwk, int *lrwk, double *rwk, int *lprwk, 
             int *ltrwk)
{
  int i;
  *ifail=0;  *ltiwk=0;  *lpiwk=1;  *lprwk=1;  *ltrwk=0;
  if (x) free(x);
  x = malloc( (*n)*sizeof(double) );
  if (!x) {printf("Allocation of memory for x failed!\n");
            *ifail=-99; return;};
  if (Ap) free(Ap);
  Ap = malloc( (*n+1)*sizeof(int) );
  if (!Ap) {printf("Allocation of memory for Ap failed!\n");
            *ifail=-99; return;};
  if (Ai) free(Ai);
  Ai = malloc( (*nfmax)*sizeof(int) );
  if (!Ai) {printf("Allocation of memory for Ai failed!\n");
            *ifail=-99; return;};
  if (Aa) free(Aa);
  Aa = malloc( (*nfmax)*sizeof(double) );
  if (!Aa) {printf("Allocation of memory for Aa failed!\n");
            *ifail=-99; return;};
  return;
}

void nifact_(int *n, int *nfill, double *a, int *irow, int *icol, int *iopt,
             int *qstruc, int *ifail, int *liwk, int *iwk, int *lrwk,
             double *rwk)
{
  int i;
  for (i=0;i<*nfill;i++) { irow[i]--; icol[i]--;};
  umfpack_di_triplet_to_col( *n, *n, *nfill, irow, icol, a, Ap, Ai, Aa, NULL );
  *ifail = umfpack_di_symbolic( *n, *n, Ap, Ai, Aa, &Symbolic, Control, Info );
  if ( *ifail != UMFPACK_OK ) {printf("umfpack_di_symbolic failed\n"); return;}; 
  *ifail = umfpack_di_numeric( Ap, Ai, Aa, Symbolic, &Numeric, Control, Info );
  if ( *ifail != UMFPACK_OK ) printf("umfpack_di_numeric failed\n"); 
  return;
}

void nisolv_(int *n, int *nfill, double *a, int *irow, int *icol, double *b,
             int *iopt, int *ifail, int *liwk, int *iwk, int *lrwk,
             double *rwk)
{
  int i;
  *ifail = umfpack_di_solve( UMFPACK_A, Ap, Ai, Aa, x, b, Numeric, 
                             Control, Info );
  if ( *ifail != UMFPACK_OK ) return; 
  for (i=0;i<*n;i++) b[i]=x[i];
  return;
}