SUBROUTINE UDU( 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 SYMMETRIC 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 WIDHT 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
C
C     BEGIN FORWARD ELIMINATION
C
      DO 450 I=1,NEM1
      JEND=NUMEQ-I+1
      IF(JEND.GT.IB) JEND=IB
      DO 440 J=2,JEND
      J1=I+J-1
      IF(LU.EQ.0) GO TO 435
      FAC=A(I,J)/A(I,1)
      K1=0
      DO 430 K=J,JEND
      K1=K1+1
      A(J1,K1)=A(J1,K1)-A(I,K)*FAC
  430 CONTINUE
      A(I,J)=FAC
  435 CONTINUE
      Y(J1)=Y(J1)-Y(I)*A(I,J)
  440 CONTINUE
  450 CONTINUE
C
C     BEGIN BACK SUBSTITUTION
C
      Y(NUMEQ)=Y(NUMEQ)/A(NUMEQ,1)
      DO 470 I=1,NEM1
      I1=NUMEQ-I
      Y(I1)=Y(I1)/A(I1,1)
      JEND=NUMEQ-I1+1
      IF(JEND.GT.IB) JEND=IB
      DO 460 J=2,JEND
      J1=I1+J-1
      Y(I1)=Y(I1)-A(I1,J)*Y(J1)
  460 CONTINUE
  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