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

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

Routine: udu - factorize a symmetric, banded matrix and solve

C-version of original UDU solver by E. Thompson

Returns: 0 - normal executing
        -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, i1, j1, k1, jend, neqm1;
  Real fac;

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

  neqm1=neq-1;

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

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