/* sparse.c     Sparse storage of a matrix
 *              specifically designed for solving PDE
 *              simultaneous equations.
 *              zero based subscripting
 *              column(nrow) is right hand side
 */

#include "sparse.h"
#include <stdio.h>
#define abs(x) ((x)<0.0?(-(x)):(x))

static int nrow;

/* newCell(j, val, next)  cell stuff, abstraction */
/* Cell operations, uses "private" part of header */

static link newCell(int aj, double aval, link anext)
{
  link next;

  next = (link)malloc(sizeof(struct cell));
  next->j = aj;
  next->val = aval;
  next->next = anext;
  return next;
} /* end newCell */

static void Setval(link* C, int aj, double aval) /* creates aj if not in list*/
{
  link* P;

  if(*C == NULL)
  {
    *C = newCell(aj, aval, NULL);
    return;
  }
  P = C; /* walk the list */
  while(1)
  {
    if(aj == (*P)->j)
    {
      (*P)->val = aval;
      return;
    }
    if(aj < (*P)->j) /* we passed it, insert before */
    {
      (*P)->next = newCell((*P)->j, (*P)->val, (*P)->next);
      (*P)->j = aj;
      (*P)->val = aval;
      return ;
    }
    if((*P)->next == NULL) /* insert at end */
    {
      (*P)->next = newCell(aj, aval, NULL);
      return;
    }
    P = &((*P)->next); /* move along list */
  } /* end while */
} /* end sparseSetval */

static void Addval(link* C, int aj, double aval) /* creates aj if not in list*/
{
  link* P;

  if(*C == NULL)
  {
    *C = newCell(aj, aval, NULL);
    return;
  }
  P = C; /* walk the list */
  while(1)
  {
    if(aj == (*P)->j)
    {
      (*P)->val = (*P)->val + aval;
      return;
    }
    if(aj < (*P)->j) /* we passed it, insert before */
    {
      (*P)->next = newCell((*P)->j, (*P)->val, (*P)->next);
      (*P)->j = aj;
      (*P)->val = aval;
      return ;
    }
    if((*P)->next == NULL) /* insert at end */
    {
      (*P)->next = newCell(aj, aval, NULL);
      return;
    }
    P = &((*P)->next); /* move along list */
  } /* end while */
} /* end Addval  aval added */

static double Getval(link C, int aj)
{
  link P;

  P = C; /* walk the list */
  while(P != NULL)
  {
    if(aj == P->j) return P->val;
    if(aj <  P->j) return 0.0;  /* we passed it */
    P = P->next; /* move along list */
  } /* end while */
  return 0.0;
} /* end Getval */


/* exported via sparse.h */

void sparseInit(struct sparse* A, int nrow)
{
  int i;

  A->nrow = nrow;
  A->A = (link *)malloc(nrow*sizeof(link));
  for(i=0; i<nrow; i++)
  {
    A->A[i] = NULL;
  }
} /* end sparseinit */

void sparseWrite_All(struct sparse A)
{
  link P;
  int i;

  nrow = A.nrow;
  printf("SparseWrite_All nrow=%d \n", nrow);
  for(i=0; i<nrow; i++)
  {
    P = A.A[i];
    while(P != NULL)
    {
      printf("i=%d, j=%d, val=%g \n", i, P->j, P->val );
      P = P->next;
    }
  }
} /* end sparseWrite_All */

void sparseImport(struct sparse A, double M[]) /* single index i*nrow+j */
{
  int i, j;

  nrow = A.nrow;
  for(i=0; i<nrow; i++)
    for(j=0; j<=nrow; j++)
      if(M[i*nrow+j] != 0.0) Setval(&(A.A[i]),j,M[i*nrow+j]);
} /* end sparseimport */

void sparseExport(struct sparse A, double M[]) /* size nrow*(nrow-1) */
{
  link P;
  int i, j;

  nrow = A.nrow;
  for(i=0; i<nrow; i++)
    for(j=0; j<=nrow; j++)
      M[i*nrow+j] = 0.0;

  for(i=0; i<nrow; i++)
  {
    P = A.A[i];
    while(P != NULL)
    {
      M[i*nrow+P->j] = P->val;
      P = P->next;
    }
  }
} /* end sparseExport */

double sparseGet(struct sparse A, int i, int j)
{
  nrow = A.nrow;
  if(i<0 || i>=nrow || j<0 || j>nrow)
  {
    printf("sparseGet error i=%d, j=%d \n", i, j);
    return 0.0;
  }
  return Getval(A.A[i],j);
} /* end sparseGet */

void sparsePut(struct sparse A, int i, int j, double val)
{
  nrow = A.nrow;
  if(i<0 || i>=nrow || j<0 || j>nrow)
  {
    printf("sparsePut error i=%d, j=%d \n", i, j);
    return;
  }
  Setval(&(A.A[i]), j, val);
} /* end sparse Put */

void sparseAdd(struct sparse A, int i, int j, double val)
{
  nrow = A.nrow;
  if(i<0 || i>=nrow || j<0 || j>nrow)
  {
    printf("sparseAdd error i=%d, j=%d \n", i, j);
    return;
  }
  Addval(&(A.A[i]), j, val);
} /* end sparseAdd */

void sparseGetRHS(struct sparse A, double Y[])
{
  int i;

  nrow = A.nrow;
  for(i=0; i<nrow; i++)
      Y[i] = Getval(A.A[i],nrow);
} /* end sparseGetRHS */

void sparseSetRHS(struct sparse A, double Y[])
{
  int i;

  nrow = A.nrow;
  for(i=0; i<nrow; i++) Setval(&(A.A[i]),nrow,Y[i]);
} /* end sparseSetRHS */

void sparseMul(struct sparse A, double X[], double Y[]) /* Y output */
{
  link P;
  int i;

  nrow = A.nrow;
  for(i=0; i<nrow; i++) Y[i] = 0.0;
  for(i=0; i<nrow; i++)
  {
    P = A.A[i];
    while(P != NULL)
    {
      if(P->j >= nrow) break;
      Y[i] = Y[i] + P->val*X[P->j];
      P = P->next;
    }
  }
} /* end sparseMul */

/* sparseSimeq      solve simultaneous equations AX=Y
 * solve real linear equations for X where Y = A * X
 * method: Gauss-Jordan elimination using maximum pivot
 *    Translated to Sparse by : Jon Squire , 25 March 2009
 *    First written by Jon Squire December 1959 for IBM 650, translated to
 *    other languages  e.g. Fortran converted to Ada converted to C
 *   ) converted to java,) to sparse.c
 */
void sparseSimeq(struct sparse A, double X[]) /* X returned solution */
{
  int Hold, i_Pivot, i, k;
  double   Pivot, Abs_Pivot;
  int  N  = nrow;
  int row[N]; /* row interchange indices */

  /* set up row interchange vectors */
  for(k=0; k<N; k++) row[k] = k;

  /*  begin main reduction loop */
  for(k=0; k<N; k++)
  {
    /* find largest element for pivot */
    Pivot = Getval(A.A[row[k]],k); /* B(row(k))(k);  */
    Abs_Pivot = abs(Pivot);
    i_Pivot = k;
    for(i=k+1; i<N; i++)
    {
      /* if A(row(i),k) = NULL) exit; end if; */
      Pivot = Getval(A.A[row[i]],k); /* B(row(i))(k)) */
      if(abs(Pivot) > Abs_Pivot) /* B(row(i))(k)) */
      {
        i_Pivot = i;
        Abs_Pivot = abs(Pivot);
      }
    }
    /* have pivot, interchange row indices
    Hold = row[k];
    row[k] = row[i_Pivot];
    row[i_Pivot] = Hold;
    /* check for near singular */
    if(Abs_Pivot < 1.0E-10)
    {
      printf("singular at k=%d redundant row=%d \n", k, row[k]);
    }
    else
    {
      /* can divide by pivot */
      /* reduce about pivot */
      /* for(int j=k+1; j<n+1; j++) */
      /* { */
      /* B(row(k))(j) = B(row(k))(j) / B(row(k))(k); */
      sparseDivide(A,row[k],k); /* j loop in here  should use pivot */
      /* } */
      /*  inner reduction loop */
      for(i=0; i<N; i++)
      {
        if(i!= k)
	{
          /* for(int j=k+1; j<n+1; j++) */
          /* { */
          /* B(row(i))(j) = B(row(i))(j) - B(row(i))(k) * B(row(k))(j); */
          sparseReduce(A,row[i],k,row[k]); /* j loop in here on both */
	                                   /* row(i) and row(k) */
        }
      }
      /*  finished inner reduction */
    }
  }  /*  end main reduction loop */
    /*  build  X  for return, unscrambling rows */
  for(i=0; i<N; i++) X[i] = Getval(A.A[row[i]],N); 
                            /* this could be direct if matrix shrunk */
} /* end sparseSimeq */

void sparseDivide(struct sparse A, int krow, int k) 
                  /* works on a krow j=k+1..n */
{
  link* P;
  double Piv;

  P = &(A.A[krow]);
  while((*P) != NULL) /* find pivot, position to pivot */
  {
    if((*P)->j == k) /* found, thus must exists */
    {
      Piv = (*P)->val;
      P = &((*P)->next);
      break;
    }
    P = &((*P)->next);
  } /* end while that gets Piv and positions to next */
  while((*P) != NULL)
  {
    (*P)->val = (*P)->val/Piv;
    P = &((*P)->next);
  }
} /*  end sparseDivide */

void sparseReduce(struct sparse A, int irow, int k, int krow)
{
  /* works on complete irow  j=k+1..n */
  /* using krow j=k+1..n  where it exists */
  /* B[row[i]][j] = B[row[i]][j] - B[row[i]][k] * B[row[k]][j]; */
  /*                               valik */
  link Pij, Pkj;
  double valik = 0.0;

  Pkj = A.A[krow];
  Pij = A.A[irow];
  /* find Pik and valik  constant for  while loop on j */
  while(Pij != NULL)
  {
    if(Pij->j == k)
    {
      valik = Pij->val; /* have multiplier */
      Pij = Pij->next;
      break;
    }
    else if(Pij->j > k)
    {
      return; /* no multiplier */
    }
    else
    {
      Pij = Pij->next;
    }
  }
  if(valik == 0.0) return; /* no multiplier */
  while(Pkj != NULL) /* position Pkj */
  {
    if(Pkj->j > k)
    {
      break;
    }
    Pkj = Pkj->next;
  }
  while(Pij!=NULL && Pkj!= NULL)
  {
    /* find Pij->j = Pkj,j (both move in j) */
    if(Pij->j < Pkj->j)
    {
      Pij = Pij->next;
    }
    else if(Pij->j == Pkj->j)
    {
      Pij->val = Pij->val - valik * Pkj->val;
      Pij = Pij->next;
      Pkj = Pkj->next;
    }
    else if(Pij->j > Pkj->j) /* must insert negative value, Pij->val = 0 */
    {
      Pij->next = newCell(Pij->j, Pij->val, Pij->next);
      Pij->j = Pkj->j;
      Pij->val = -valik*Pkj->val;
      Pij = Pij->next;
      Pkj = Pkj->next;
    }
    else
    {
      Pkj = Pkj->next;
    }
  }
} /* end sparseReduce */
/* end sparse.c */