C
C     ..................................................................
C
C        SUBROUTINE DMCHB
C
C        PURPOSE
C           FOR A GIVEN POSITIVE-DEFINITE M BY M MATRIX A WITH SYMMETRIC
C           BAND STRUCTURE AND - IF NECESSARY - A GIVEN GENERAL M BY N
C           MATRIX R, THE FOLLOWING CALCULATIONS (DEPENDENT ON THE
C           VALUE OF THE DECISION PARAMETER IOP) ARE PERFORMED
C           (1) MATRIX A IS FACTORIZED (IF IOP IS NOT NEGATIVE), THAT
C               MEANS BAND MATRIX TU WITH UPPER CODIAGONALS ONLY IS
C               GENERATED ON THE LOCATIONS OF A SUCH THAT
C               TRANSPOSE(TU)*TU=A.
C           (2) MATRIX R IS MULTIPLIED ON THE LEFT BY INVERSE(TU)
C               AND/OR INVERSE(TRANSPOSE(TU)) AND THE RESULT IS STORED
C               IN THE LOCATIONS OF R.
C           THIS SUBROUTINE ESPECIALLY CAN BE USED TO SOLVE THE SYSTEM
C           OF SIMULTANEOUS LINEAR EQUATIONS A*X=R WITH POSITIVE-
C           DEFINITE COEFFICIENT MATRIX A OF SYMMETRIC BAND STRUCTURE.
C
C        USAGE
C           CALL DMCHB (R,A,M,N,MUD,IOP,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           R      - INPUT IN CASES IOP=-3,-2,-1,1,2,3  DOUBLE PRECISION
C                          M BY N RIGHT HAND SIDE MATRIX,
C                          IN CASE IOP=0  IRRELEVANT.
C                    OUTPUT IN CASES IOP=1,-1  INVERSE(A)*R,
C                           IN CASES IOP=2,-2  INVERSE(TU)*R,
C                           IN CASES IOP=3,-3  INVERSE(TRANSPOSE(TU))*R,
C                           IN CASE  IOP=0     UNCHANGED.
C           A      - INPUT IN CASES IOP=0,1,2,3  DOUBLE PRECISION M BY M
C                          POSITIVE-DEFINITE COEFFICIENT MATRIX OF
C                          SYMMETRIC BAND STRUCTURE STORED IN
C                          COMPRESSED FORM (SEE REMARKS),
C                          IN CASES IOP=-1,-2,-3 DOUBLE PRECISION M BY M
C                          BAND MATRIX TU WITH UPPER CODIAGONALS ONLY,
C                          STORED IN COMPRESSED FORM (SEE REMARKS).
C                    OUTPUT IN ALL CASES  BAND MATRIX TU WITH UPPER
C                           CODIAGONALS ONLY, STORED IN COMPRESSED FORM
C                           (THAT MEANS UNCHANGED IF IOP=-1,-2,-3).
C           M      - INPUT VALUE SPECIFYING THE NUMBER OF ROWS AND
C                    COLUMNS OF A AND THE NUMBER OF ROWS OF R.
C           N      - INPUT VALUE SPECIFYING THE NUMBER OF COLUMNS OF R
C                    (IRRELEVANT IN CASE IOP=0).
C           MUD    - INPUT VALUE SPECIFYING THE NUMBER OF UPPER
C                    CODIAGONALS OF A.
C           IOP    - ONE OF THE VALUES -3,-2,-1,0,1,2,3 GIVEN AS INPUT
C                    AND USED AS DECISION PARAMETER.
C           EPS    - SINGLE PRECISION INPUT VALUE USED AS RELATIVE
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANT DIGITS.
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C                     IER=0  - NO ERROR,
C                     IER=-1 - NO RESULT BECAUSE OF WRONG INPUT
C                              PARAMETERS M,MUD,IOP (SEE REMARKS),
C                              OR BECAUSE OF A NONPOSITIVE RADICAND AT
C                              SOME FACTORIZATION STEP,
C                              OR BECAUSE OF A ZERO DIAGONAL ELEMENT
C                              AT SOME DIVISION STEP.
C                     IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
C                              CANCE INDICATED AT FACTORIZATION STEP K+1
C                              WHERE RADICAND WAS NO LONGER GREATER
C                              THAN EPS*A(K+1,K+1).
C
C        REMARKS
C           UPPER PART OF SYMMETRIC BAND MATRIX A CONSISTING OF MAIN
C           DIAGONAL AND MUD UPPER CODIAGONALS (RESP. BAND MATRIX TU
C           CONSISTING OF MAIN DIAGONAL AND MUD UPPER CODIAGONALS)
C           IS ASSUMED TO BE STORED IN COMPRESSED FORM, I.E. ROWWISE
C           IN TOTALLY NEEDED M+MUD*(2M-MUD-1)/2 SUCCESSIVE STORAGE
C           LOCATIONS. ON RETURN UPPER BAND FACTOR TU (ON THE LOCATIONS
C           OF A) IS STORED IN THE SAME WAY.
C           RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE
C           IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN RESULT MATRIX
C           INVERSE(A)*R OR INVERSE(TU)*R OR INVERSE(TRANSPOSE(TU))*R
C           IS STORED COLUMNWISE TOO ON THE LOCATIONS OF R.
C           INPUT PARAMETERS M, MUD, IOP SHOULD SATISFY THE FOLLOWING
C           RESTRICTIONS     MUD NOT LESS THAN ZERO,
C                            1+MUD NOT GREATER THAN M,
C                            ABS(IOP) NOT GREATER THAN 3.
C           NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE
C           RESTRICTIONS ARE NOT SATISFIED.
C           THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT
C           PARAMETERS ARE SATISFIED, IF RADICANDS AT ALL FACTORIZATION
C           STEPS ARE POSITIVE AND/OR IF ALL DIAGONAL ELEMENTS OF
C           UPPER BAND FACTOR TU ARE NONZERO.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           FACTORIZATION IS DONE USING CHOLESKY-S SQUARE-ROOT METHOD,
C           WHICH GENERATES THE UPPER BAND MATRIX TU SUCH THAT
C           TRANSPOSE(TU)*TU=A. TU IS RETURNED AS RESULT ON THE
C           LOCATIONS OF A. FURTHER, DEPENDENT ON THE ACTUAL VALUE OF
C           IOP, DIVISION OF R BY TRANSPOSE(TU) AND/OR TU IS PERFORMED
C           AND THE RESULT IS RETURNED ON THE LOCATIONS OF R.
C           FOR REFERENCE, SEE H. RUTISHAUSER, ALGORITHMUS 1 - LINEARES
C           GLEICHUNGSSYSTEM MIT SYMMETRISCHER POSITIV-DEFINITER
C           BANDMATRIX NACH CHOLESKY - , COMPUTING (ARCHIVES FOR
C           ELECTRONIC COMPUTING), VOL.1, ISS.1 (1966), PP.77-78.
C
C     ..................................................................
C
      SUBROUTINE DMCHB(R,A,M,N,MUD,IOP,EPS,IER)
C
C
      DIMENSION R(1),A(1)
      DOUBLE PRECISION TOL,SUM,PIV,R,A
C
C        TEST ON WRONG INPUT PARAMETERS
      IF(IABS(IOP)-3)1,1,43
    1 IF(MUD)43,2,2
    2 MC=MUD+1
      IF(M-MC)43,3,3
    3 MR=M-MUD
      IER=0
C
C        MC IS THE MAXIMUM NUMBER OF ELEMENTS IN THE ROWS OF ARRAY A
C        MR IS THE INDEX OF THE LAST ROW IN ARRAY A WITH MC ELEMENTS
C
C     ******************************************************************
C
C        START FACTORIZATION OF MATRIX A
      IF(IOP)24,4,4
    4 IEND=0
      LLDST=MUD
      DO 23 K=1,M
      IST=IEND+1
      IEND=IST+MUD
      J=K-MR
      IF(J)6,6,5
    5 IEND=IEND-J
    6 IF(J-1)8,8,7
    7 LLDST=LLDST-1
    8 LMAX=MUD
      J=MC-K
      IF(J)10,10,9
    9 LMAX=LMAX-J
   10 ID=0
      TOL=A(IST)*EPS
C
C        START FACTORIZATION-LOOP OVER K-TH ROW
      DO 23 I=IST,IEND
      SUM=0.D0
      IF(LMAX)14,14,11
C
C        PREPARE INNER LOOP
   11 LL=IST
      LLD=LLDST
C
C        START INNER LOOP
      DO 13 L=1,LMAX
      LL=LL-LLD
      LLL=LL+ID
      SUM=SUM+A(LL)*A(LLL)
      IF(LLD-MUD)12,13,13
   12 LLD=LLD+1
   13 CONTINUE
C        END OF INNER LOOP
C
C        TRANSFORM ELEMENT A(I)
   14 SUM=A(I)-SUM
      IF(I-IST)15,15,20
C
C        A(I) IS DIAGONAL ELEMENT. ERROR TEST.
   15 IF(SUM)43,43,16
C
C        TEST ON LOSS OF SIGNIFICANT DIGITS AND WARNING
   16 IF(SUM-TOL)17,17,19
   17 IF(IER)18,18,19
   18 IER=K-1
C
C        COMPUTATION OF PIVOT ELEMENT
   19 PIV=DSQRT(SUM)
      A(I)=PIV
      PIV=1.D0/PIV
      GO TO 21
C
C        A(I) IS NOT DIAGONAL ELEMENT
   20 A(I)=SUM*PIV
C
C        UPDATE ID AND LMAX
   21 ID=ID+1
      IF(ID-J)23,23,22
   22 LMAX=LMAX-1
   23 CONTINUE
C
C        END OF FACTORIZATION-LOOP OVER K-TH ROW
C        END OF FACTORIZATION OF MATRIX A
C
C     ******************************************************************
C
C        PREPARE MATRIX DIVISIONS
      IF(IOP)24,44,24
   24 ID=N*M
      IEND=IABS(IOP)-2
      IF(IEND)25,35,25
C
C     ******************************************************************
C
C        START DIVISION BY TRANSPOSE OF MATRIX TU (TU IS STORED IN
C        LOCATIONS OF A)
   25 IST=1
      LMAX=0
      J=-MR
      LLDST=MUD
      DO 34 K=1,M
      PIV=A(IST)
      IF(PIV)26,43,26
   26 PIV=1.D0/PIV
C
C        STA-T BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
      DO 30 I=K,ID,M
      SUM=0.D0
      IF(LMAX)30,30,27
C
C        PREPARE INNER LOOP
   27 LL=IST
      LLL=I
      LLD=LLDST
C
C        START INNER LOOP
      DO 29 L=1,LMAX
      LL=LL-LLD
      LLL=LLL-1
      SUM=SUM+A(LL)*R(LLL)
      IF(LLD-MUD)28,29,29
   28 LLD=LLD+1
   29 CONTINUE
C        END OF INNER LOOP
C
C        TRANSFORM ELEMENT R(I)
   30 R(I)=PIV*(R(I)-SUM)
C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
C
C        UPDATE PARAMETERS LMAX, IST AND LLDST
      IF(MC-K)32,32,31
   31 LMAX=K
   32 IST=IST+MC
      J=J+1
      IF(J)34,34,33
   33 IST=IST-J
      LLDST=LLDST-1
   34 CONTINUE
C
C        END OF DIVISION BY TRANSPOSE OF MATRIX TU
C
C     ******************************************************************
C
C        START DIVISION BY MATRIX TU (TU IS STORED ON LOCATIONS OF A)
      IF(IEND)35,35,44
   35 IST=M+(MUD*(M+M-MC))/2+1
      LMAX=0
      K=M
   36 IEND=IST-1
      IST=IEND-LMAX
      PIV=A(IST)
      IF(PIV)37,43,37
   37 PIV=1.D0/PIV
      L=IST+1
C
C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
      DO 40 I=K,ID,M
      SUM=0.D0
      IF(LMAX)40,40,38
   38 LLL=I
C
C        START INNER LOOP
      DO 39 LL=L,IEND
      LLL=LLL+1
   39 SUM=SUM+A(LL)*R(LLL)
C        END OF INNER LOOP
C
C        TRANSFORM ELEMENT R(I)
   40 R(I)=PIV*(R(I)-SUM)
C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R
C
C        UPDATE PARAMETERS LMAX AND K
      IF(K-MR)42,42,41
   41 LMAX=LMAX+1
   42 K=K-1
      IF(K)44,44,36
C
C        END OF DIVISION BY MATRIX TU
C
C     ******************************************************************
C
C        ERROR EXIT IN CASE OF WRONG INPUT PARAMETERS OR PIVOT ELEMENT
C        LESS THAN OR EQUAL TO ZERO
   43 IER=-1
   44 RETURN
      END
