      SUBROUTINE BPROD(M,N,L,A,KA,ML,MU,B,KB,NL,NU,C,KC,NC,MCL,MCU,IERR)
C     ******************************************************************
C     MULTIPLICATION OF REAL BANDED MATRICES
C     ******************************************************************
      REAL A(KA,N),B(KB,L),C(KC,NC)
      DOUBLE PRECISION DSUM
C
      IERR = 0
      MLP1 = ML + 1
      NLP1 = NL + 1
      NPML = N + ML
      NPNU = N + NU
      MCL = MIN0(M-1,ML+NL)
      IF (MCL .EQ. 0) GO TO 100
C
C     FIND THE FIRST NONZERO LOWER DIAGONAL
C
      MAXD = MCL
      DO 21 NDIAG = 1,MAXD
      IMJ = MAXD + 1 - NDIAG
      JMAX = MIN0(L,M-IMJ,NPML-IMJ)
        DO 11 J = 1,JMAX
        I = J + IMJ
        DSUM = 0.0
        IF (J .GT. NPNU) GO TO 11
        KMIN = MAX0(1,I-ML,J-NU)
        KMAX = MIN0(N,I+MU,J+NL)
        KK = MLP1 - I + KMIN
        JJ = NLP1 + J - KMIN
          DO 10 K = KMIN,KMAX
          DSUM = DSUM + DBLE(A(I,KK))*DBLE(B(K,JJ))
          KK = KK + 1
   10     JJ = JJ - 1
   11   C(I,1) = DSUM
C
      JMAX = MIN0(JMAX,NPNU)
        DO 20 J = 1,JMAX
        I = J + IMJ
        IF (C(I,1)) 30,20,30
   20   CONTINUE
   21 MCL = MCL - 1
      GO TO 100
C
   30 IF (MCL .GE. NC) GO TO 200
      C(1,1) = 0.0
      IF (MCL .EQ. 1) GO TO 100
C
C     COMPUTE THE REMAINING LOWER DIAGONALS
C
      JC = 1
      MIND = NDIAG + 1
      DO 42 NDIAG = MIND,MAXD
      JC = JC + 1
      IMJ = MAXD + 1 - NDIAG
      JMAX = MIN0(L,M-IMJ,NPML-IMJ)
        DO 41 J = 1,JMAX
        I = J + IMJ
        DSUM = 0.0
        IF (J .GT. NPNU) GO TO 41
        KMIN = MAX0(1,I-ML,J-NU)
        KMAX = MIN0(N,I+MU,J+NL)
        KK = MLP1 - I + KMIN
        JJ = NLP1 + J - KMIN
          DO 40 K = KMIN,KMAX
          DSUM = DSUM + DBLE(A(I,KK))*DBLE(B(K,JJ))
          KK = KK + 1
   40     JJ = JJ - 1
   41   C(I,JC) = DSUM
   42 CONTINUE
C
C     INSERT ZEROS IN THE UPPER LEFT CORNER
C
      IMAX = MCL
      DO 51 J = 1,MCL
        DO 50 I = 1,IMAX
   50   C(I,J) = 0.0
   51 IMAX = IMAX - 1
C
C     FIND THE LAST NONZERO UPPER DIAGONAL
C
  100 JC = MCL + 1
      MCU = MIN0(L-1,MU+NU)
      IF (MCU .EQ. 0) GO TO 140
C
      MAXD = MCU
      DO 121 NDIAG = 1,MAXD
      JMI = MAXD + 1 - NDIAG
      IMAX = MIN0(M,L-JMI,NPML)
        DO 111 I = 1,IMAX
        J = I + JMI
        DSUM = 0.0
        IF (J .GT. NPNU) GO TO 111
        KMIN = MAX0(1,I-ML,J-NU)
        KMAX = MIN0(N,I+MU,J+NL)
        KK = MLP1 - I + KMIN
        JJ = NLP1 + J - KMIN
          DO 110 K = KMIN,KMAX
          DSUM = DSUM + DBLE(A(I,KK))*DBLE(B(K,JJ))
          KK = KK + 1
  110     JJ = JJ - 1
  111   C(I,JC) = DSUM
C
      IMAX1 = MIN0(IMAX,NPNU-JMI)
        DO 120 I = 1,IMAX1
        IF (C(I,JC)) 130,120,130
  120   CONTINUE
  121 MCU = MCU - 1
      GO TO 140
C
  130 LAST = JC + MCU
      IF (LAST .GT. NC) GO TO 210
      DO 131 I = 1,IMAX
  131 C(I,LAST) = C(I,JC)
C
C     COMPUTE THE MAIN DIAGONAL AND THE REMAINING UPPER DIAGONALS
C
  140 MAXD = MAX0(1,MCU)
      DO 143 NDIAG = 1,MAXD
      JMI = NDIAG - 1
      IMAX = MIN0(M,L-JMI,NPML)
        DO 142 I = 1,IMAX
        J = I + JMI
        DSUM = 0.0
        IF (J .GT. NPNU) GO TO 142
        KMIN = MAX0(1,I-ML,J-NU)
        KMAX = MIN0(N,I+MU,J+NL)
        KK = MLP1 - I + KMIN
        JJ = NLP1 + J - KMIN
          DO 141 K = KMIN,KMAX
          DSUM = DSUM + DBLE(A(I,KK))*DBLE(B(K,JJ))
          KK = KK + 1
  141     JJ = JJ - 1
  142   C(I,JC) = DSUM
  143 JC = JC + 1
C
C     INSERT ZEROS IN THE LOWER RIGHT CORNER
C
      JMAX = MCL + MCU + 1
      IMIN = L - MCU + 1
      IMAX = MIN0(M,NPML)
      IF (IMIN .GT. IMAX) GO TO 160
C
      JMIN = MAX0(1,JMAX-IMAX+IMIN)
      J = JMAX
      DO 151 JJ = JMIN,JMAX
        DO 150 I = IMIN,IMAX
  150   C(I,J) = 0.0
      IMIN = IMIN + 1
  151 J = J - 1
C
C     STORE ZEROS IN THE FINAL M-IMAX ROWS
C
  160 IF (IMAX .EQ. M) RETURN
      IMIN = IMAX + 1
      DO 162 J = 1,JMAX
        DO 161 I = IMIN,M
  161   C(I,J) = 0.0
  162 CONTINUE
      RETURN
C
C     ERROR RETURN - C REQUIRES AT LEAST IERR COLUMNS
C
  200 IERR = MCL + 1
      RETURN
  210 IERR = LAST
      RETURN
      END
