#include <stdio.h>
#include "types.h"

int
ldu( Real *a,        /* coefficient matrix in ldu format */
     Real *y,        /* RHS on input, solution on output */
     const int neq,  /* No. of equations to solve        */
     const int ib,   /* bandwidth of the system          */
     const int lu    /* lu=1 factorize, lu=0 resolve     */ )
/***********************************************************************

Routine: ldu - factorize a non-symmetric, banded matrix and solve

C(TRAN)-version of original LDU solver by E. Thompson

Returns: 0 - normal execution
        -1 - failed due to an error in calling sequence

Author : Mark A. Christon

Date   : 01/08/2000

Date     |   Modification
=========|==============================================================
xx/xx/xx |

***********************************************************************/
{
  int i, j, k, irow, jrow, jcol, i1, j1, k1, k2; 
  int jbgn, jend, neqm1, idiag, ibd2;
  Real fac;

  /* Error trap */
  if (neq==0) {
    fprintf(stdout,"\n\tLDU: Number of equations is invalid ... ");
    fprintf(stdout,"Neq = %d\n", neq);
    return(-1);
  }
  if (ib==0) {
    fprintf(stdout,"\n\tLDU: bandwidth is invalid ... ");
    fprintf(stdout,"Ib = %d\n",ib);
    return(-1);
  }
  if (a==NULL) {
    fprintf(stdout,"\n\tLDU: Null pointer for coefficient matrix ...\n");
    return(-1);
  }
  if (y==NULL) {
    fprintf(stdout,"\n\tLDU: Null pointer for rhs vector ...\n");
    return(-1);
  }

  neqm1=neq-1;
  ibd2 =((ib-1)/2)-1;
  idiag=ibd2+1;

  /* Begin forward elimination */
  for (i=0;i<neqm1;i++) {
    jend=neqm1-i; 
    if (jend>ibd2) jend=ibd2;
    for (j=0;j<=jend;j++) {
      jrow=i+j+1;
      jcol=ibd2-j;
      if (lu==1) {
	fac=a[jrow+jcol*neq]/a[i+idiag*neq];
        a[jrow+jcol*neq]=fac;
	for (k=0;k<=jend;k++) {
	  k1=idiag+k+1;
          k2=jcol+k+1;
          a[jrow+k2*neq]=a[jrow+k2*neq]-fac*a[i+k1*neq];
	}
      }
      y[jrow]-=y[i]*a[jrow+jcol*neq];
    } /* end of j-loop */
  } /* end of i-loop */

  /* Back substitution */
  y[neqm1]/=a[neqm1+idiag*neq];
  jbgn=idiag+1;
  for (i=0;i<neqm1;i++) {
    irow=neqm1-i-1;
    jend=idiag+i+1;
    if (jend>ib-1) jend=ib-1;
    jrow=irow;
    for (j=jbgn;j<=jend;j++) {
      jrow=jrow+1;
      y[irow]-=a[irow+j*neq]*y[jrow];
    }
    y[irow]=y[irow]/a[irow+idiag*neq];
  }
}