SUBROUTINE LDU( A, Y, LU, NUMEQ, IB, IDIM, JDIM, LLPT)
C***********************************************************************
C                                                                      *
C     CE 665 FINITE ELEMENT METHOD                                     *
C     ERIK THOMPSON                                                    *
C     COLORADO STATE UNIVERSITY                                        *
C     FORT COLLINS, COLORADO  80523                                    *
C                                                                      *
C     A SUBROUTINE FOR THE FACTORIZATION OF A SYMMETRIC BANDED         *
C     MATRIX AND THE SOLUTION OF THE MATRIX EQUATION                   *
C                                                                      *
C                        AX=F                                          *
C                                                                      *
C                                                                      *
C                                                                      *
C     A     = NAME OF MATRIX ( UNFACTORED IF LU=1, FACTORED IF LU=0 )  *
C             NOTE: THE A-MATRIX MUST BE IN THE PROPER BANDED FOR      *
C             USED FOR NONSYMMETRIC MATRICES.                          *
C     Y     = RIGHT HAND SIDE WHEN SUBROUTINE IS CALLED FROM           *
C             MAIN PROGRAM. (I.E.  Y=F)                                *
C           = SOLUTION OF MATRIX EQUATION WHEN SUBROUTINE IS           *
C             RETURNED TO MAIN PROGRAM.                                *
C     LU    = 1   MEANS THAT THE A-MATRIX IS IN ORIGINAL FORM.         *
C                 FACTORIZATION IS DESIRED AS WELL AS THE              *
C                 SOLUTION VECTOR, Y=X.                                *
C     LU    = 0   MEANS THAT THE A-MATRIX AS ENTERED IS ALREADY        *
C                 FACTORED AND ALL THAT IS DESIRED IS THE SOLUTION     *
C                 VECTOR, Y=X.                                         *
C     NUMEQ = THE NUMBER OF EQUATIONS REPRESENTED BY THE MATRIX        *
C             EQUATION.  NOTE:  THIS MUST BE EQUAL TO OR LESS          *
C             THAN THE FIRST DIMENSION OF THE A-MATRIX.                *
C     IB    = THE 'HALF' BAND WIDTH OF THE SET OF EQAUTIONS.           *
C             NOTE:  THIS VALUE MUST BE EQUAL TO OR LESS THAN THE      *
C             SECOND DIMENSION OF THE A-MATRIX.                        *
C     IDIM  = FIRST DIMENSION OF THE A-MATRIX AS SPECIFIED IN THE      *
C             DIMENSION STATEMENT OF THE CALLING PROGRAM.              *
C     JDIM  = SECOND DIMENSION OF THE A-MATRIX AS SPECIFIED IN THE     *
C             DIMENSION STATEMENT OF THE CALLING PROGRAM.              *
C     LLPT  = TAPE NUMBER USED FOR OUTPUT (USUALLY EQUALS 6 )          *
C                                                                      *
C***********************************************************************
#ifdef DBL
      implicit double precision (a-h,o-z)
#endif
      DIMENSION A(IDIM,JDIM), Y(IDIM)
C
      IF(NUMEQ.GT.IDIM) GO TO 7001
      IF(IB   .GT.JDIM) GO TO 7002
C
      NEM1=NUMEQ-1
      IBD2=(IB-1)/2
C
C     BEGIN FORWARD ELIMINATION
C
      IDIAG=IBD2+1
C
      DO 450 I=1,NEM1
        JEND=NUMEQ-I
        IF(JEND.GT.IBD2) JEND=IBD2
        DO 440 J=1,JEND
          JROW=I+J
          JCOL=IBD2-J+1
          IF(LU.EQ.0) GO TO 435
          FAC=A(JROW,JCOL)/A(I,IDIAG)
          A(JROW,JCOL)=FAC
          DO 430 K=1,JEND
            K1=IDIAG+K
            K2=JCOL+K
            A(JROW,K2)=A(JROW,K2)-FAC*A(I,K1)
 430      CONTINUE
 435      CONTINUE
          Y(JROW)=Y(JROW)-A(JROW,JCOL)*Y(I)
 440    CONTINUE
 450  CONTINUE
C
C     BEGIN BACK SUBSTITUTION
C
      Y(NUMEQ)=Y(NUMEQ)/A(NUMEQ,IDIAG)
C
      JBGN=IDIAG+1
      DO 470 I=1,NEM1
        IROW=NUMEQ-I
        JEND=IDIAG+I
        IF(JEND.GT.IB) JEND=IB
        JROW=IROW
        DO 460 J=JBGN,JEND
          JROW=JROW+1
          Y(IROW)=Y(IROW)-A(IROW,J)*Y(JROW)
 460    CONTINUE
        Y(IROW)=Y(IROW)/A(IROW,IDIAG)
 470  CONTINUE
C
      RETURN
C
C     ERROR MESAAGES
C
 7001 WRITE(LLPT,1) NUMEQ,IDIM
      STOP
 7002 WRITE(LLPT,2) IB,JDIM
      STOP
C
    1 FORMAT('0 NUMBER OF EQUATIONS EXCEEDS IDIM',/,
     2       ' NUMEQ = ',I5,3X,'IDIM = ',I5)
    2 FORMAT('0 BANDWIDTH EXCEEDS JDIM',/,
     2       ' IB = ',I5,3X,'JDIM = ',I5)
C
      END