SUBROUTINE DBCDEC(W,MP1,N, AUX)  ! decomposition, bandwidth M
*     Symmetric band matrix: decomposition, solution, complete
*                                inverse and band part of the inverse
*     M=MP1-1                                    N*M(M-1) dot operations
      DOUBLE PRECISION W(MP1,N),V(MP1,N),B(N),X(N),AUX(N),VS(*),RXW
*     ...
      DO I=1,N
       AUX(I)=16.0D0*W(1,I) ! save diagonal elements
      END DO
      DO I=1,N
       IF(W(1,I)+AUX(I).NE.AUX(I)) THEN
          W(1,I)=1.0/W(1,I)
       ELSE
          W(1,I)=0.0D0  ! singular
       END IF
       DO J=1,MIN(MP1-1,N-I)
        RXW=W(J+1,I)*W(1,I)            
        DO K=1,MIN(MP1-1,N-I)+1-J
         W(K,I+J)=W(K,I+J)-W(K+J,I)*RXW 
        END DO
        W(J+1,I)=RXW
       END DO
      END DO
      RETURN

      ENTRY DBCSLV(W,MP1,N,B, X)  ! solution B -> X
*                                                N*(2M-1) dot operations
      DO I=1,N
       X(I)=B(I)
      END DO 
      DO I=1,N ! forward substitution
       DO J=1,MIN(MP1-1,N-I)
        X(J+I)=X(J+I)-W(J+1,I)*X(I)   
       END DO
      END DO
      DO I=N,1,-1 ! backward substitution
       RXW=X(I)*W(1,I)                 
       DO J=1,MIN(MP1-1,N-I)
        RXW=RXW-W(J+1,I)*X(J+I)         
       END DO
       X(I)=RXW
      END DO
      RETURN

      ENTRY DBCIEL(W,MP1,N, V)    ! V = inverse band matrix elements
*                                               N*M*(M-1) dot operations
      DO I=N,1,-1
       RXW=W(1,I)
       DO J=I,MAX(1,I-MP1+1),-1
        DO K=J+1,MIN(N,J+MP1-1)
         RXW=RXW-V(1+ABS(I-K),MIN(I,K))*W(1+K-J,J)
        END DO
        V(1+I-J,J)=RXW
        RXW=0.0D0
       END DO  
      END DO
      RETURN

      ENTRY DBCINV(W,MP1,N, VS)   ! V = inverse symmetric matrix
*                                             N*N/2*(M-1) dot operations
      DO I=N,1,-1 
       RXW=W(1,I) 
       DO J=I,1,-1
        DO K=J+1,MIN(N,J+MP1-1)
         RXW=RXW-VS((MAX(I,K)*(MAX(I,K)-1))/2+MIN(I,K))*W(1+K-J,J)
        END DO 
        VS((I*I-I)/2+J)=RXW
        RXW=0.0D0
       END DO
      END DO
      END