      SUBROUTINE DBCMAV(N,MA,MB,BCV,IHOM,NRTI,Q,NTI,PHIREC,NU,Z,KKK,
     1                  R,BB,WS)
C     --------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION MA(N,N),MB(N,N)
      DIMENSION BCV(N),Q(N,N,NTI),PHIREC(NU,NTI),Z(N,NTI),
     1          R(N,N),BB(N),WS(N,N)
      INTEGER KKK(N)
      LOGICAL HOM0,HOM1
C
      HOM0 = IHOM.EQ.0
      HOM1 = IHOM.NE.0
      DO 1000 I = 1 , N
        KKK(I) = I*(I+1)/2
 1000 CONTINUE
      DO 1010 I = 1 , N
        BB(I) = BCV(I)
 1010 CONTINUE
C     COMPUTATION OF MA.Q(1) STORED IN WS
      DO 1120 K = 1 , N
        CALL DMATVC(MA,N,N,N,N,Q(1,K,1),WS(1,K))
 1120 CONTINUE
      IF (HOM1) THEN
C     COMPUTATION OF MA.Q(1).Z(1) AND BCV-MA.Q(1).Z(1) STORED IN BB.
        DO 1210 I = 1 , N
          SOM = 0.D0
          DO 1220 J = 1 , N
            SOM = SOM + WS(I,J) * Z(J,1)
 1220     CONTINUE
          BB(I) = BB(I) - SOM
 1210   CONTINUE
      ENDIF
C     COMPUTATION OF MA.Q(1).PHIREC(1) STORED IN R
      DO 1420 K = 1 , N
        NK = KKK(K) - K
        DO 1410 I = 1 , N
          SOM = 0.D+0
          DO 1400 J = 1 , K
            NJ = NK + J
            SOM = SOM + WS(I,J) * PHIREC(NJ,1)
 1400     CONTINUE
          R(I,K) = SOM
 1410   CONTINUE
 1420 CONTINUE
C     COMPUTATION OF MB.Q(N) STORED IN WS
      DO 1520 K = 1 , N
        CALL DMATVC(MB,N,N,N,N,Q(1,K,NRTI),WS(1,K))
 1520 CONTINUE
      IF (HOM1) THEN
C     COMPUTATION OF MB.Q(N).Z(N) AND BCV-MB.Q(N).Z(N) STORED IN BB
        DO 1610 I = 1 , N
          SOM = 0.D0
          DO 1600 J = 1 , N
            SOM = SOM + WS(I,J) * Z(J,NRTI)
 1600     CONTINUE
          BB(I) = BB(I) - SOM
 1610   CONTINUE
      ENDIF
C     COMPUTATION OF MB.Q(N).PHIREC(N) AND
C     MA.Q(1).PHIREC(1) + MB.Q(N).PHIREC(N) STORED IN R
      DO 1820 K = 1 , N
        NK = KKK(K) - K
        DO 1810 I = 1 , N
          SOM = 0.D+0
          DO 1800 J = 1 , K
            NJ = NK + J
            SOM = SOM + WS(I,J) * PHIREC(NJ,NRTI)
 1800     CONTINUE
        R(I,K) = R(I,K) + SOM
 1810   CONTINUE
 1820 CONTINUE
      RETURN
C     END OF DBCMAV
      END
      SUBROUTINE DCROUT(A,N,M,B,EPS,P,V,IER)
C     --------------------------------------
C
      DOUBLE PRECISION A(M,1),B(1),V(1),EPS
      INTEGER P(1)
C
      IER = 0
      CALL DLUDEC(A,N,M,EPS,P,V,IER)
      IF (IER.NE.0) RETURN
      CALL DSOLDE(A,N,M,P,B)
      RETURN
C     END OF DCROUT
C
C     DCROUT IS AN DOUBLE PRECISION FORTRAN VERSION OF
C     T.DEKKER'S SOLDEC
      END
      SUBROUTINE DFQUS(PHI,N,IRM,IRN,DIAG,BET,K,KKK,U,NU,NTI,Q,
     1                 WPHI0,WSUM)
C     ----------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PHI(IRM,IRN),DIAG(IRN),BET(IRN),U(NU,NTI),
     1          Q(IRM,IRM,NTI),WPHI0(IRM,IRM),WSUM(IRN)
      INTEGER KKK(N)
C
      CALL DQEVAL(N,N,PHI,IRM,IRN,BET,WPHI0,WSUM)
      DO 1000 I = 1 , N
      DO 1000 J = 1 , N
      Q(I,J,K) = WPHI0(I,J)
 1000 CONTINUE
      U(1,K) = DIAG(1)
      IF (N.LT.2) GOTO 2000
      DO 1200 I = 2 , N
      NI = KKK(I)
      IMIN1 = I - 1
      NJ = NI - I
      U(NI,K) = DIAG(I)
      DO 1100 J = 1 , IMIN1
      NK = NJ + J
      U(NK,K) = PHI(J,I)
 1100 CONTINUE
 1200 CONTINUE
 2000 DO 2300 I = 1 , N
        NI = KKK(I)
        IF (U(NI,K).GE.0.D0) GOTO 2300
        DO 2100 J = 1 , N
          Q(J,I,K) = -Q(J,I,K)
 2100   CONTINUE
        DO 2200 J = I , N
          NJ = KKK(J) - J + I
          U(NJ,K) = -U(NJ,K)
 2200   CONTINUE
 2300 CONTINUE
      RETURN
C     END OF DFQUS
      END
      SUBROUTINE DFUNRC(N,KPART,NRTI,U,NU,NTI,KKK,PHIREC,W1,W2,
     1                  IERROR)
C     ----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION PHIREC(NU,NTI),U(NU,NTI),W1(N),W2(N)
      INTEGER KKK(N)
C
      IERROR = 0
      IF (KPART.LT.0) KPART=0
      IF (KPART.GT.N) KPART=N
      KPPL1 = KPART + 1
C     COMPUTATION OF THE PHI1'S FOR I <= KPART
      IF (KPART.GT.0) THEN
      DO 1400 K = 1 , KPART
      KK = KKK(K)
      NK = KK - K
      NK1 = NK + 1
      DO 1150 I = 1 , K
      PHIREC(NK + I,NRTI) = 0.0D+0
 1150 CONTINUE
      PHIREC(KK,NRTI) = 1.0D+0
      DO 1350 I = 2 , NRTI
      IN = NRTI + 1 - I
      IN1 = IN + 1
      DO 1250 IJ = 1 , K
      W2(IJ) = PHIREC(NK + IJ,IN1)
 1250 CONTINUE
      CALL DSOLUP(K,U(1,IN1),W2,W1,IERROR)
      IF (IERROR.GT.0) GOTO 3000
      DO 1300 IJ = 1 , K
      PHIREC(NK + IJ,IN) = W1(IJ)
 1300 CONTINUE
 1350 CONTINUE
 1400 CONTINUE
      ENDIF
C     THE REQUIRED PART OF PHIREC1 IS NOW COMPUTED.
C     THE REQUIRED PART OF PHIREC2 WILL BE COMPUTED NOW.
      IF (KPPL1.LE.N) THEN
C     COMPUTATION OF THE PHI2'S FOR I > KPART
      DO 1700 K = KPPL1 , N
      KK = KKK(K)
      NK = KK - K
      NK1 = NK + KPPL1
      DO 1500 I = NK1 , KK
      PHIREC(I,1) = 0.0D+0
 1500 CONTINUE
      PHIREC(KK,1) = 1.0D+0
      DO 1650 I = 2 , NRTI
      IN1 = I - 1
      DO 1600 IJ = KPPL1 , K
      SOM = 0.0D+0
      DO 1550 IN = IJ , K
      MK = KKK(IN) - IN
      SOM = SOM + PHIREC(NK + IN , IN1) * U(MK + IJ , I)
 1550 CONTINUE
      PHIREC(NK + IJ , I) = SOM
 1600 CONTINUE
 1650 CONTINUE
 1700 CONTINUE
C     COMPUTATION OF THE PHI1'S FOR I > KPART, IF KPART > 0
      IF (KPART.GT.0) THEN
      DO 2050 K = KPPL1 , N
      KK = KKK(K)
      NK = KK - K
      DO 1750 I = 1, KPART
      PHIREC(NK + I, NRTI) = 0.0D+0
 1750 CONTINUE
      DO 2000 I = 2 , NRTI
      IN = NRTI + 1 - I
      IN1 = IN + 1
      DO 1850 IJ = 1 , KPART
      SOM = 0.0D+0
      DO 1800 NK1 = KPPL1 , K
      MK = KKK(NK1 - 1)
      SOM = SOM + U(MK + IJ,IN1) * PHIREC(NK + NK1 ,IN)
 1800 CONTINUE
      W2(IJ) = - SOM + PHIREC(NK + IJ , IN1)
 1850 CONTINUE
      CALL DSOLUP(KPART,U(1,IN1),W2,W1,IERROR)
      IF (IERROR.GT.0) GOTO 3000
      DO 1950 IJ = 1 , KPART
      PHIREC(NK + IJ , IN) = W1(IJ)
 1950 CONTINUE
 2000 CONTINUE
 2050 CONTINUE
      ENDIF
      ENDIF
      RETURN
 3000 IERROR = 250
      RETURN
C     END OF DFUNRC
      END
      SUBROUTINE DGTUR(N,IHOM,NRTI,U,KU,NTI,Q,D,ER,KKK,KPART,IP1,IP2,
     1                 WDIA,WBET,WSUM,WPHI,WPHI0,IERROR)
C     ----------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(KU,NTI),Q(N,N,NTI),D(N,NTI),ER(5),WDIA(N),
     1          WBET(N),WSUM(N),WPHI(N,N),WPHI0(N,N)
      INTEGER KKK(N),IP1(N),IP2(N)
      LOGICAL NOSORT
C
      IERROR = 0
      ISORT = IP1(1)
      NKP = N
      IF (ISORT.NE.0) NKP = IP1(2)
      IF (NKP.LE.0) THEN
        KPART = 0
        RETURN
      ENDIF
      NTI1 = NTI - 1
      NOSORT = .FALSE.
C.....CHECK WHETHER TRANSFORMATION IS ALLOWED.
      TOL = DMAX1(ER(1),ER(2))
      TOL1 = TOL / ER(3)
      DO 1100 I = 1 , N
        IP1(I) = I
 1100 CONTINUE
      DO 1300 I = 1 , NKP
        IK = KKK(I)
        WDIA(I) = 0.0D+0
        DO 1200 J = 2 , NRTI
          SOM = DABS(U(IK,J))
          IF (SOM.GE.TOL1) NOSORT = .TRUE.
          WDIA(I) = WDIA(I) + DLOG10(SOM)
 1200 CONTINUE
 1300 CONTINUE
      KPART = 0
      DO 1400 I = 1 , NKP
        IF (WDIA(I).GT.0) KPART = KPART + 1
 1400 CONTINUE
C     PARTITIONING IS DETERMINED
      IF ((KPART.EQ.0).OR.(KPART.EQ.NKP)) RETURN
      IF (NOSORT) RETURN
      IFLAG = 0
      K = 0
      N1 = NKP - 1
C     DETERMINATION OF THE PERMUTATION WHICH GLOBALLY ORDERS THE MODES
      DO 1700 I = 1 , N1
        SORT = WDIA(I)
        IF (SORT.LE.0.D0) THEN
          IPL1 = I + 1
          DO 1500 J = IPL1 , NKP
            SOM = WDIA(J)
            IF ((SOM.GT.0.D0).AND.(SOM-SORT.GT.1.D0)) THEN
              IFLAG = 1
              WDIA(J) = SORT
              WDIA(I) = SOM
              JM = IP1(I)
              IP1(I) = IP1(J)
              IP1(J) = JM
              GOTO 1600
            ENDIF
 1500     CONTINUE
 1600     CONTINUE
        ENDIF
 1700 CONTINUE
      IF (IFLAG.EQ.0) RETURN
C.....THE PERMUTATION P IS NOW DETERMINED
C.....COMPUTATION OF Q1P,U2P,INV(P)D1,INV(P)W1
      IERROR = 1
 2000 DO 2100 I = 1 , N
      DO 2100 J = 1 , N
        Q(I,J,NTI) = Q(I,J,1)
        WPHI(I,J) = 0.0D+0
 2100 CONTINUE
      DO 2500 I = 1 , N
        IK = IP1(I)
        IK1 = KKK(IK) - IK
        DO 2300 J = 1 , N
           Q(J,I,1) = Q(J,IK,NTI)
 2300   CONTINUE
        DO 2400 J = 1 , IK
          WPHI(J,I) = U(IK1+J,2)
 2400   CONTINUE
 2500 CONTINUE
C.....TRANSFORMATION OF THE U'S,Q'S,D'S AND WI'S
      DO 3400 I = 2 , NRTI
        CALL DQUDEC(WPHI,N,N,N,N,WDIA,WBET,WSUM)
        CALL DFQUS(WPHI,N,N,N,WDIA,WBET,NTI,KKK,U,KU,NTI,Q,
     1             WPHI0,WSUM)
        DO 2700 J = 1 , KU
          U(J,I) = U(J,NTI)
 2700   CONTINUE
        IF (IHOM.EQ.0) GOTO 2950
        CALL DTAMVC(Q(1,1,NTI),N,N,N,N,D(1,I),WBET)
        DO 2900 J = 1 , N
          D(J,I) = WBET(J)
 2900   CONTINUE
 2950   DO 3000 J = 1 , N
          CALL DMATVC(Q(1,1,I),N,N,N,N,Q(1,J,NTI),WPHI0(1,J))
 3000   CONTINUE
        DO 3100 J = 1 , N
        DO 3100 L = 1 , N
          Q(J,L,I) = WPHI0(J,L)
 3100   CONTINUE
        I1 = I + 1
        IF (I1.GT.NRTI) GOTO 3400
C.....COMPUTATION OF UI+1.OI
        DO 3300 J = 1 , N
        DO 3300 L = 1 , N
          SOM = 0.0D+0
          DO 3200 K = L , N
            IK = KKK(K) - K
            SOM = SOM + U(IK+L,I1) * Q(K,J,NTI)
 3200     CONTINUE
          WPHI(L,J) = SOM
 3300  CONTINUE
 3400 CONTINUE
      RETURN
C     END OF DGTUR
      END
       SUBROUTINE DINPRO(A,NA,B,NB,N,C)
C      --------------------------------
C
       DOUBLE PRECISION A(1),B(1),C
C
C      LOCAL VARIABLES
       INTEGER NN,JL,JJ
       C = 0.0D+0
       IF (N.LE.0) GOTO 2000
       NN = 1 - NB
       JL = NA * N
       DO 1000 JJ = 1 , JL , NA
        NN = NN + NB
        C = C + A(JJ) * B(NN)
 1000  CONTINUE
 2000  CONTINUE
       RETURN
C      END OF DINPRO
C
C      DINPRO IS THE DOUBLE PRECISION VERSION OF SUBROUTINE
C      VVIPP FROM THE ACCU LIBRARY
C
      END
      SUBROUTINE DINTCH(A,NA,B,NB,N)
C     ------------------------------
C
      DOUBLE PRECISION A(1),B(1),C
C
      IF (N.LE.0) RETURN
      NN = 1 - NB
      JL = NA * N
      DO 1000 JJ = 1 , JL , NA
        NN = NN + NB
        C = B(NN)
        B(NN) = A(JJ)
        A(JJ) = C
 1000 CONTINUE
      RETURN
C     END OF DINTCH
C
C     DINTCH IS A DOUBLE PRECISION VERSION OF
C     J.M. VAN KATS SUBROUTINE INTCHA
      END
      SUBROUTINE DKPCH(N,U,NU,NTI,NRTI,KPART,KKK,ER,IER)
C     --------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NU,NTI),ER(5)
      INTEGER KKK(N)
C
      IER = 0
      PSI = DMAX1(ER(1),ER(2)) / ER(3)
      ER(5) = 0.0D+0
      BLOWUQ = 0.0D+0
      KPPL1 = KPART + 1
      IF (KPART.GE.N) GOTO 1300
      DO 1200 K = KPPL1 , N
        PROD = 1.0D+0
        KK = KKK(K)
        DO 1100 I = 2 , NRTI
          PROD = 1.0D+0 + PROD * DABS(U(KK,I))
          IF (BLOWUQ.LT.PROD) BLOWUQ = PROD
 1100   CONTINUE
 1200 CONTINUE
 1300 IF (KPART.LT.1) GOTO 2000
      DO 1500 K = 1 , KPART
        PROD = 1.0D+0
        KK = KKK(K)
        DO 1400 I = 2 , NRTI
          IN = NRTI + 2 - I
          PROD = 1.0D+0 + PROD / DABS(U(KK,IN))
          IF (ER(5).LT.PROD) ER(5) = PROD
 1400   CONTINUE
 1500 CONTINUE
      IF (KPART.GE.N) GOTO 2000
      BLOW1 = 0.0D+0
      BLOW2 = 0.0D+0
      DO 1700 K = 1 , KPART
        KK = KKK(K)
        AMX = 1.0D+0 / DABS(U(KK,NRTI))
        AMX2 = AMX
        IF (NRTI.LT.3) GOTO 1650
        DO 1600 I = 3 , NRTI
        IN = NRTI + 2 - I
        PROD = 1.0D+0 / DABS(U(KK,IN))
        AMX1 = PROD * AMX2
        AMX2 = DMAX1(AMX1,PROD)
        AMX = DMAX1(AMX,AMX2)
 1600   CONTINUE
 1650   IF (AMX.GT.BLOW1) BLOW1 = AMX
 1700 CONTINUE
      DO 1900 K = KPPL1 , N
        KK = KKK(K)
        AMX = DABS(U(KK,2))
        AMX2 = AMX
        DO 1800 I = 3 , NRTI
        PROD = DABS(U(KK,I))
        AMX1 = PROD * AMX2
        AMX2 = DMAX1(PROD,AMX1)
        AMX = DMAX1(AMX,AMX2)
 1800   CONTINUE
        IF (AMAX.GT.BLOW2) BLOW2 = AMX
 1900 CONTINUE
      BLOW = BLOW1 * BLOW2
      ER(5) = ER(5) + BLOW
 2000 ER(5) = DMAX1(ER(5),BLOWUQ)
      IF (ER(5).LT.PSI) RETURN
      IER = 240
      RETURN
C     END OF DKPCH
      END
      SUBROUTINE DLUDEC(A,N,M,EPS,P,V,IER)
C     ------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(M,1),V(1)
      INTEGER P(1),PK
C
      IER = 0
      R = -1.0D+0
      Z = 1.0D+0
      DO 1000 I = 1 , N
        CALL DINPRO(A(I,1),M,A(I,1),M,N,S)
        S = DSQRT(S)
        IF (S.GT.R) R = S
        V(I) = 1.0D+0 / S
 1000 CONTINUE
      E = EPS * R
      DO 5000 K = 1 , N
        R = -1.0D+0
        K1 = K - 1
        DO 2000 I = K , N
          CALL DINPRO(A(I,1),M,A(1,K),1,K1,S)
          A(I,K) = A(I,K) - S
          S = DABS(A(I,K)) * V(I)
          IF (S.LE.R) GOTO 2000
          R = S
          PK = I
 2000   CONTINUE
        P(K) = PK
        V(PK) = V(K)
        Z = A(PK,K)
        IF (DABS(Z).LT.E) GOTO 7000
        IF (PK.EQ.K) GOTO 3000
        CALL DINTCH(A(K,1),M,A(PK,1),M,N)
 3000   KK = K + 1
        IF (K.EQ.N) RETURN
        DO 4000 I = KK , N
          CALL DINPRO(A(K,1),M,A(1,I),1,K1,S)
          A(K,I) = (A(K,I) - S) / Z
 4000   CONTINUE
 5000 CONTINUE
 7000 IER = 1
      RETURN
C     END OF DLUDEC
C
C     DLUDEC IS AN FORTRAN VERSION OF T. DEKKER'S DEC
      END
      SUBROUTINE DMATVC(A,M,N,IRM,IRN,Y,X)
C     ------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(IRM,IRN),Y(IRN),X(IRM)
C
      IF (M.LE.0.OR.N.LE.0) GOTO 2000
      DO 1100 I=1,M
      S = 0.0D+0
      DO 1000 J=1,N
      S = S + A(I,J) * Y(J)
 1000 CONTINUE
      X(I) = S
 1100 CONTINUE
 2000 CONTINUE
      RETURN
C     END OF DMATVC
      END
      SUBROUTINE DPSR(N,KPART,NRTI,U,NU,NTI,D,KKK,V,W1,W2,IERROR)
C     -----------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(NU,NTI),D(N,NTI),V(N,NTI),W1(N),W2(N)
      INTEGER KKK(N)
C
      IERROR = 0
      KPPL1 = KPART + 1
C     COMPUTATION OF VI2
      IF (KPART.GE.N) GOTO 1500
      DO 1300 I = KPPL1 , N
      V(I,1) = 0.0D+0
 1300 CONTINUE
      DO 1450 I = 2 , NRTI
      IN = I - 1
      DO 1400 K = KPPL1 , N
      SOM = 0.0D+0
      DO 1350 L = K , N
      MK = KKK(L) - L
      SOM = SOM + U(MK+K,I) * V(L,IN)
      MK1 = MK + K
 1350 CONTINUE
      V(K,I) = SOM + D(K,I)
 1400 CONTINUE
 1450 CONTINUE
C     COMPUTATION OF VI1
 1500 IF (KPART.LT.1) RETURN
      DO 1600 I = 1 , KPART
      V(I,NRTI) = 0.0D+0
 1600 CONTINUE
      DO 1850 I = 2 , NRTI
      IN = NRTI + 1 - I
      IN1 = IN + 1
      DO 1700 K = 1 , KPART
      SOM = 0.0D+0
      IF (KPPL1.GT.N) GOTO 1660
      DO 1650 L = KPPL1 , N
      MK = KKK(L) - L
      SOM = SOM + U(MK+K,IN1) * V(L,IN)
 1650 CONTINUE
 1660 W2(K) = - SOM - D(K,IN1) + V(K,IN1)
 1700 CONTINUE
      CALL DSOLUP(KPART,U(1,IN1),W2,W1,IERROR)
      IF (IERROR.NE.0) GOTO 3000
      DO 1800 K = 1 , KPART
      V(K,IN) = W1(K)
 1800 CONTINUE
 1850 CONTINUE
      RETURN
 3000 IERROR = 250
      RETURN
C     END OF DPSR
      END
      SUBROUTINE DQEVAL(M,N,QR,IRM,IRN,BET,QREV,WD)
C     --------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION QR(IRM,IRN),QREV(IRM,IRM),BET(IRN),WD(IRM)
C
      DO 1100 I=N,M
      DO 1050 J=I,M
        QREV(I,J) = -QR(I,N) * QR(J,N) * BET(N)
        QREV(J,I) = QREV(I,J)
 1050 CONTINUE
      QREV(I,I) = QREV(I,I) + 1.0D+0
 1100 CONTINUE
      IF (N.LT.2) RETURN
      DO 1300 IK = 2 , N
      I = N + 2 - IK
C     THE NECESARY INNERPRODUCTS OF THE UK'S ARE STORED IN THE
C     ARRAY WD(1:IRM)
C
      IMIN1 = I - 1
      WD(IMIN1) = QR(IMIN1,IMIN1)
      IM = M - I + 1
      DO 1150 J = I , M
        CALL DINPRO(QR(I,IMIN1),1,QREV(I,J),1,IM,WD(J))
C
C     COMPUTATION OF (I - BET(I-1).U(I-1)') * PRODUCT
C     ( ( -BET(K).U(K).U(K)') K= I TO N)
C
 1150 CONTINUE
      DO 1200 K = I , M
      DO 1200 L = I , M
        QREV(K,L) = QREV(K,L) - BET(IMIN1) * QR(K,IMIN1) * WD(L)
 1200 CONTINUE
      DO 1250 J = I , M
        QREV(J,IMIN1) = - BET(IMIN1) * QR(J,IMIN1) * WD(IMIN1)
        QREV(IMIN1,J) = - BET(IMIN1) * QR(IMIN1,IMIN1) * WD(J)
 1250 CONTINUE
      QREV(IMIN1,IMIN1) = 1.0D+0 - BET(IMIN1)*WD(IMIN1)*WD(IMIN1)
 1300 CONTINUE
      RETURN
C     END OF DQEVAL
      END
      SUBROUTINE DQUDEC(A,M,N,IRM,IRN,DIAG,BET,WY)
C     ----------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(IRM,IRN),DIAG(IRN),BET(IRN),WY(IRN)
C
      IF (M.GT.IRM) GOTO 3000
      IF (N.GT.IRN) GOTO 3010
C
      DO 2000 K = 1 , N
       K1 = K + 1
       IK = M - K + 1
       CALL DINPRO(A(K,K),1,A(K,K),1,IK,SIGMA)
       IF (SIGMA.EQ.0.D0) GOTO 1400
       QRKK = A(K,K)
       ALPHAK = DSQRT(SIGMA)
       IF (QRKK.GE.0.0D+0) ALPHAK = - ALPHAK
       DIAG(K) = ALPHAK
       BETA = 1.0D+0 / (SIGMA - QRKK * ALPHAK)
       BET(K) = BETA
       A(K,K) = QRKK - ALPHAK
       IF (K1.GT.N) GOTO 2000
       DO 1200 J = K1 , N
        CALL DINPRO(A(K,K),1,A(K,J),1,IK,WY(J))
        WY(J) = WY(J) * BETA
 1200  CONTINUE
       DO 1300 J = K1 , N
       DO 1300 I = K  , M
        A(I,J) = A(I,J) - A(I,K) * WY(J)
 1300  CONTINUE
       GOTO 2000
 1400  BET(K) = 0.D0
       DIAG(K) = 0.D0
 2000 CONTINUE
      RETURN
 3000 WRITE(*,100)
      RETURN
 3010 WRITE(*,110)
      RETURN
C
  100 FORMAT(' ERROR: M > IRM ; NO DECOMPOSITION DONE')
  110 FORMAT(' ERROR: N > IRN ; NO DECOMPOSITION DONE')
C
C     END OF DQUDEC
      END
      SUBROUTINE DSBVP(N,IHOM,M0,MN,NRTI,Q,NTI,BCV,PHIREC,NU,V,
     1                 WI,EPS,COND,KKK,IP,WQ,WS,W1,W2,IER)
C    ------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION M0(N,N),MN(N,N)
      DIMENSION Q(N,N,NTI),BCV(N),PHIREC(NU,NTI),V(N,NTI),
     1          WI(N,NTI),WQ(N,N),WS(N,N),W1(N),W2(N)
      INTEGER KKK(N),IP(N)
C
  801 FORMAT(' ',5(D12.5,2X))
      CALL DBCMAV(N,M0,MN,BCV,IHOM,NRTI,Q,NTI,PHIREC,NU,WI,KKK,WQ,
     1            W1,WS)
      CALL DCROUT(WQ,N,N,W1,EPS,IP,W2,IER)
      IF (IER.NE.0) THEN
       IER = 260
       RETURN
      ENDIF
C     COMPUTATION OF THE CONDITION NUMBER
      COND = 0.0D+0
      DO 1300 I = 1 , N
        DO 1100 J = 1 , N
          W2(J) = 0.0D+0
 1100   CONTINUE
        W2(I) = 1.0D+0
        CALL DSOLDE(WQ,N,N,IP,W2)
        SOMIN = 0.0D+0
        DO 1200 J = 1 , N
          SOMIN = SOMIN + DABS(W2(J))
 1200   CONTINUE
        IF (SOMIN.GT.COND) COND = SOMIN
 1300 CONTINUE
C     COMPUTATION OF EI'=PHIREC(I)Y' AND EI=Q(I)EI'
      DO 1500 K=1,NRTI
      DO 1500 I=1,N
        SOMIN=0.0D+0
        DO 1400 J=I,N
          NJ=KKK(J) - J
          NI= NJ + I
          SOMIN = SOMIN + PHIREC(NI,K) * W1(J)
 1400   CONTINUE
        IF (IHOM.EQ.0) THEN
          WI(I,K) = SOMIN
        ELSE
          WI(I,K) = WI(I,K) + SOMIN
        ENDIF
 1500 CONTINUE
C     COMPUTATION OF Q(I)XI'
      DO 2000 K=1,NRTI
        CALL DMATVC(Q(1,1,K),N,N,N,N,WI(1,K),W2)
        DO 1800 J=1,N
          V(J,K) = W2(J)
 1800   CONTINUE
 2000 CONTINUE
      RETURN
C     END OF DSBVP
      END
      SUBROUTINE DSOLUP(M,U,B,X,IER)
C     ------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(1),X(1),B(1)
C
      IER = 0
      IF (M.LE.0) RETURN
      M1 = M * (M+1) / 2
      IF (U(M1).EQ.0.0D+0) GOTO 1200
      X(M) = B(M) / U(M1)
      IF (M.EQ.1) RETURN
      MMIN1 = M - 1
      DO 1100 I = 1 , MMIN1
      MI = M - I
      SOMIN = 0.0D+0
      MJ = (MI + 1) * MI / 2
      MIJ = MJ
      DO 1000 J = 1 , I
      MIJ = MIJ + MI - 1 + J
      SOMIN = SOMIN + U(MIJ) * X(MI+J)
 1000 CONTINUE
      IF (U(MJ).EQ.0.0D+0) GOTO 1200
      X(MI) = (B(MI) - SOMIN) / U(MJ)
 1100 CONTINUE
      RETURN
 1200 IER = 1
      RETURN
C     END OF DSOLUP
      END
      SUBROUTINE DSOLDE(A,N,M,P,B)
C     ----------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(M,1),B(1)
      INTEGER P(1),PK
C
      DO 1000 K = 1 , N
        R = B(K)
        PK = P(K)
        CALL DINPRO(A(K,1),M,B,1,K-1,C)
        B(K) = (B(PK) - C) / A(K,K)
        IF (PK.NE.K) B(PK) = R
 1000 CONTINUE
      IF (N.EQ.1) RETURN
      K = N
 2000 K = K - 1
      CALL DINPRO(A(K,K+1),M,B(K+1),1,N-K,C)
      B(K) = B(K) - C
      IF (K.GT.1) GOTO 2000
      RETURN
C     END OF DSOLDE
C
C     DSOLDE IS DOUBLE PRECISION VERSION OF H.ZWEERUS-VINK'S
C     AND J.M. VAN KAT'S SUBROUTINE SOLDEC
      END
      SUBROUTINE DTAMVC(A,M,N,IM,IN,Y,X)
C     ----------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(IM,IN),Y(IM),X(IN)
C
      IF (M.LE.0 .OR. N.LE.0) RETURN
      DO 1100 I = 1 , N
        SOM = 0.D0
        DO 1000 J = 1 , M
          SOM = SOM + A(J,I) * Y(J)
 1000   CONTINUE
        X(I) = SOM
 1100 CONTINUE
      RETURN
C     END OF DTAMVC
      END
      SUBROUTINE DSORTD(DIAG,N,NKP,KPART,IP,ISORT,IFLAG)
C     ---------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION DIAG(N)
      INTEGER IP(N)
C
      IFLAG = 0
      DO 1000 I = 1 , N
        IP(I) = I
 1000 CONTINUE
      NN = N
      IF (ISORT.EQ.1) THEN
C     ISORT = 1 ; CHECK WHETHER PERMUTATION IS NECESSARY.
        J = 0
        DO 1100 I = 1 , NKP
          IF (DABS(DIAG(I)).GT.1.D0) J = J + 1
 1100   CONTINUE
        IF (J.EQ.NKP) THEN
          KPART = J
          RETURN
        ENDIF
        NN = NKP
      ENDIF
C     ENDIF ISORT = 1
      IF (ISORT.LT.2) THEN
C     MONOTONIC ORDENING
        K = 0
        N1 = NN - 1
        DO 1400 I = 1 , N1
          SORT = 0.D0
          DO 1300 J = I , NN
            A = DABS(DIAG(J))
            IF (A.GT.SORT) THEN
              SORT = A
              JM = J
            ENDIF
 1300     CONTINUE
          IF (SORT.GT.1.D0) K = K+1
          IF (JM.NE.I) THEN
            IFLAG = IFLAG + 1
            A = DIAG(JM)
            DIAG(JM) = DIAG(I)
            DIAG(I) = A
            J = IP(JM)
            IP(JM) = IP(I)
            IP(I) = J
          ENDIF
 1400   CONTINUE
        KPART = K
        IF (DABS(DIAG(NN)).GT.1.D0) KPART = KPART + 1
        RETURN
      ENDIF
C     ENDIF ISORT < 2
C     ISORT = 2 OR 3; NO MONOTONIC ORDENING
      IF (ISORT.EQ.3) NN = NKP
      N1 = NN - 1
      DO 2200 I = 1 , N1
        A = DABS(DIAG(I))
        IF (A.LT.1.D0) THEN
          DO 2000 J = I+1 , NN
            B = DABS(DIAG(J))
            IF ((B.GT.1.D0).AND.(B/A.GE.10.D0)) THEN
              IFLAG = IFLAG + 1
              SORT = DIAG(J)
              DIAG(J) = DIAG(I)
              DIAG(I) = SORT
              JM = IP(J)
              IP(J) = IP(I)
              IP(I) = JM
              GOTO 2100
            ENDIF
 2000     CONTINUE
 2100     CONTINUE
        ENDIF
 2200 CONTINUE
      KPART = 0
      DO 2300 I = 1 , NN
        IF (DABS(DIAG(I)).GT.1.D0) KPART = KPART + 1
 2300 CONTINUE
      RETURN
C     END OF DSORTD
      END
      SUBROUTINE DRKFGS(F,NEQN,Y,T,HI,NHI,W)
C     --------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NEQN),HI(5),W(NEQN,7)
      EXTERNAL F
C
      CALL F(T,Y,W(1,1))
      DO 2000 I = 1 , NHI
        CALL DRKF1S(F,NEQN,Y,T,HI(I),W)
        DO 1000 J = 1 , NEQN
          Y(J) = W(J,7)
 1000   CONTINUE
        T = T + HI(I)
        IF (I.LT.NHI) CALL F(T,Y,W(1,1))
 2000 CONTINUE
      RETURN
C     END OF DRKFGS
      END
      SUBROUTINE DRKF1S(F,NEQN,Y,T,H,W)
C     ---------------------------------
C
C     FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NEQN),W(NEQN,7)
      EXTERNAL F
C
      CH=H/4.0D0
      DO 221 K=1,NEQN
  221   W(K,6)=Y(K)+CH*W(K,1)
      CALL F(T+CH,W(1,6),W(1,2))
C
      CH=3.0D0*H/32.0D0
      DO 222 K=1,NEQN
  222   W(K,6)=Y(K)+CH*(W(K,1)+3.0D0*W(K,2))
      CALL F(T+3.0D0*H/8.0D0,W(1,6),W(1,3))
C
      CH=H/2197.0D0
      DO 223 K=1,NEQN
  223   W(K,6) = Y(K) + CH * (1932.0D0 * W(K,1)
     1                + (7296.0D0 * W(K,3) - 7200.0D0 * W(K,2)))
      CALL F(T+12.0D0*H/13.0D0,W(1,6),W(1,4))
C
      CH=H/4104.0D0
      DO 224 K=1,NEQN
  224   W(K,6)=Y(K)+CH*((8341.0D0*W(K,1)-845.0D0*W(K,4))+
     1                    (29440.0D0*W(K,3)-32832.0D0*W(K,2)))
      CALL F(T+H,W(1,6),W(1,5))
C
      CH=H/20520.0D0
      DO 225 K=1,NEQN
  225   W(K,2)=Y(K)+CH*((-6080.0D0*W(K,1)+(9295.0D0*W(K,4)-
     1         5643.0D0*W(K,5)))+(41040.0D0*W(K,2)-28352.0D0*W(K,3)))
      CALL F(T+H/2.0D0,W(1,2),W(1,6))
C
C     COMPUTE APPROXIMATE SOLUTION AT T+H
C
      CH=H/7618050.0D0
      DO 230 K=1,NEQN
  230   W(K,7)=Y(K)+CH*((902880.0D0*W(K,1)+(3855735.0D0*W(K,4)-
     1        1371249.0D0*W(K,5)))+(3953664.0D0*W(K,3)+
     2        277020.0D0*W(K,6)))
C
      RETURN
C     END OF DRKF1S
      END
      SUBROUTINE ERRHAN(IERROR,ER,NEG)
C     --------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ER(5)
C
      IF (IERROR.EQ.100) WRITE(*,100)
      IF (IERROR.EQ.101) WRITE(*,101)
      IF (IERROR.EQ.103) WRITE(*,103)
      IF (IERROR.EQ.105) WRITE(*,105)
      IF (IERROR.EQ.106) WRITE(*,106)
      IF (IERROR.EQ.120) WRITE(*,120)
      IF (IERROR.EQ.121) WRITE(*,121)
      IF (IERROR.EQ.122) WRITE(*,122)
      IF (IERROR.EQ.123) WRITE(*,123) NEG
      IF (IERROR.EQ.200) WRITE(*,200)
      IF (IERROR.EQ.213) WRITE(*,213) ER(1)
      IF (IERROR.EQ.215) WRITE(*,215)
      IF (IERROR.EQ.218) WRITE(*,218)
      IF (IERROR.EQ.219) WRITE(*,219) NEG
      IF (IERROR.EQ.230) WRITE(*,230)
      IF (IERROR.EQ.231) WRITE(*,231)
      IF (IERROR.EQ.240) WRITE(*,240) ER(5)
      IF (IERROR.EQ.250) WRITE(*,250)
      IF (IERROR.EQ.260) WRITE(*,260)
      RETURN
  100 FORMAT(' 	NPUT ERROR: N<1 OR IHOM<0 OR NRTI<0 OR NTI<5 OR A=B OR',
     1 ' NU<N*(N+1)/2')
  101 FORMAT(' INPUT ERROR: ER(1), ER(2) OR ER(3) IS NEGATIVE')
  103 FORMAT(' INPUT ERROR: LW<8*N+2*N*N OR LIW<3*N')
  105 FORMAT(' INPUT ERROR: N<1 OR NRTI<0 OR NTI<3 OR A=B OR',
     1 ' NU<N*(N+1)/2')
  106 FORMAT(' INPUT ERROR: LW<7*N+3*N*NTI+4*N*N OR LIW<3*N+NTI')
  120 FORMAT(' INPUT ERROR: NRTI=1, BUT SUPPLIED OUTPUT POINTS IN TI ',
     1 /,' ARE NOT IN STRICT MONOTONE ORDER')
  121 FORMAT(' INPUT ERROR: NRTI=1, BUT A OR B ARE NOT INCLUDED IN THE',
     1 /,' GIVEN OUTPUTPOINTS IN TI')
  122 FORMAT(' INPUT ERROR: NTI TOO SMALL !')
  123 FORMAT(' INPUT ERROR: LWG IS LESS THAN THE NUMBER OF OUTPUT  ',
     1 'POINTS.',/,' LWG MUST BE GREATER THAN:',I10)
  200 FORMAT(' INCREMENT ON MINOR SHOOTING INTERVAL IS',
     1 ' GREATER THAN AMP')
  213 FORMAT(' ER(1) TOO SMALL, CHANGED INTO:',D16.9)
  215 FORMAT(' PURE RELATIVE ERROR TEST IMPOSSIBLE; MUST USE NON-ZERO',
     1 /,' ABSOLUTE TOLERANCE')
  218 FORMAT(' N, ER(1) OR ER(2) IS NEGATIVE')
  219 FORMAT(' ARRAY WG IS TOO SMALL. ESTIMATED VALUE FOR LWG:',I10)
  230 FORMAT(' NECESSARY DAMPING FACTOR IS TOO SMALL; YOU MAY TRY MORE',
     1 /,' OUTPUT POINTS OR START WITH A BETTER APPROXIMATION OF',/,
     2 ' THE SOLUTION.')
  231 FORMAT(' NUMBER OF ITERATIONS > ITLIM')
  240 FORMAT(' BAD DICHOTOMY; AMPLIFICATION FACTOR = ',D12.5)
  250 FORMAT(' ONE OF THE UPPER TRIANGULAR MATRIX OF THE MULTIPLE',
     1 /,' SHOOTING RECURSION IS SINGULAR')
  260 FORMAT(' PROBLEM TOO ILL-CONDITIONED WITH RESPECT TO THE BC')
C     END OF ERRHAN
      END
