#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define max(a,b) ((a>=b)?a:b)
#define min(a,b) ((a>=b)?b:a)

typedef int   ITYPE;
typedef float DTYPE;
typedef float ATYPE;
typedef float VTYPE;


void lu_decomp(neqn,nrow,level,ierr,A)

    ITYPE neqn,nrow,level,*ierr;
    ATYPE A[][nrow];

{
    ITYPE i,j,k,nm1,kb,nband,kpi;
    VTYPE zero,pick,mult;

    if(level>0) printf("LU_DECOMP: neqn=%d nrow=%d\n",neqn,nrow);

    if((nrow%2)==0){
        printf("lu_decomp.c: Wrong nrow=%d\n",nrow);
        *ierr=1; return;
    }

    nband=nrow/2;
    nm1=neqn-1;
    zero=0.0;

    /*------------------*/
    /* LU Factorization */
    /*------------------*/

    for(k=0;k<neqn;k++){

        kb=min((nm1-k),nband);

        for(i=1;i<=kb;i++){
            kpi=k+i;
            pick=A[kpi][nband-i];
            /* if(pick.eq.zero) break; */
            mult=pick/A[k][nband];
            A[kpi][nband-i]=mult;
            for(j=1;j<=kb;j++)
                A[kpi][nband-i+j]-=mult*A[k][nband+j];
        }

    }
}


void substitute(neqn,nrow,level,ierr,A,b)

    ITYPE neqn,nrow,level,*ierr;
    ATYPE A[][nrow],b[];

{
    ITYPE i,j,k,nm1,kb,nband,kpi;
    VTYPE zero,pick,mult;

    if(level>0) printf("SUBSTITUTE: neqn=%d nrow=%d\n",neqn,nrow);

    if((nrow%2)==0){
        printf("substitute.c: Wrong nrow=%d\n",nrow);
        *ierr=1; return;
    }

    nband=nrow/2;
    nm1=neqn-1;
    zero=0.0;

   /*---------------------*/
   /* Forward Elimination */
   /*---------------------*/

    for(k=1;k<neqn;k++){
        kb=min(k,nband);
        for(j=1;j<=kb;j++)
            b[k]-=A[k][nband-j]*b[k-j];
    }

   /*-------------------*/
   /* Back Substitution */
   /*-------------------*/

    b[nm1]/=A[nm1][nband];

    for(k=neqn-2;k>=0;k--){
        kb=min((nm1-k),nband);
        for(j=1;j<=kb;j++)
            b[k]-=A[k][nband+j]*b[k+j];
        b[k]/=A[k][nband];
    }

}