SUBROUTINE FACTOR(A, B, C, G, H, N)                               FAC   10
C THE SUBROUTINE IS A NORMALISED FACTORISATION OF A SQUARE MATRIX OF
C ORDER GREATER THAN 4.
C THE COEFFICIENT MATRIX P IS SYMMETRIC,POSITIVE DEFINITE AND
C QUINDIAGONAL WITH NON-ZERO ELEMENTS ALSO IN THE LAST TWO COLUMNS OF
C THE FIRST ROW, THE LAST COLUMN OF THE SECOND ROW, THE FIRST COLUMN
C OF ROW (N-1) AND THE FIRST TWO COLUMNS OF THE LAST ROW. A,B,C ARE
C VECTORS EACH OF N ELEMENTS CONTAINING RESPECTIVELY THE SUPER-SUPER-
C DIAGONAL, SUPER-DIAGONAL AND DIAGONAL ELEMENTS OF P I.E. P(I,I+2),
C P(I,I+1) AND P(I,I). THE MATRIX ELEMENTS P(1,N-1), P(1,N), P(2,N)
C ARE STORED IN A(N-1), B(N), A(N) RESPECTIVELY, AS ARE P(N-1,1),
C P(N,1) AND P(N,2)
C THE MATRIX P IS FACTORISED INTO P=DRTD WHERE D IS A DIAGONAL MATRIX
C AND T IS A REAL, UPPER TRIANGULAR MATRIX WITH UNIT DIAGONAL ELEMENTS
C AND NON-ZERO ELEMENTS IN THE LAST TWO COLUMNS AND ON THE SUPER- AND
C SUPER-SUPER-DIAGONALS. R DENOTES THE TRANSPOSE OF T.
C THE INPUT VECTORS ARE OVERWRITTEN DURING THE FACTORISATION STAGE SO
C THAT VECTOR C CONTAINS THE DIAGONAL ELEMENTS OF D, VECTOR B CONTAINS
C THE ELEMENTS  T(I,I+1) FOR I=1 TO N-3  AND VECTOR A CONTAINS THE
C ELEMENTS  T(I,I+2)  FOR I=1 TO N-4.
C TWO RESULT VECTORS G,H ARE USED TO STORE THE ELEMENTS  T(I,N-1) FOR
C I=1 TO N-4  AND  T(I,N) FOR I=1 TO N-3 RESPECTIVELY. THE ELEMENTS
C T(N-3,N-1), T(N-2,N-1), T(N-2,N) AND T(N-1,N) ARE THEN GIVEN BY
C A(N-3)+G(N-3), B(N-2)+G(N-2), A(N-2)+H(N-2), B(N-1)+H(N-1)
C RESPECTIVELY.
      DIMENSION A(N), B(N), C(N), G(N), H(N)
      N1 = N - 1
      N2 = N - 2
      N3 = N - 3
      N4 = N - 4
      C(1) = SQRT(C(1))
      V = B(1)/C(1)
      C(2) = SQRT(C(2)-V*V)
      B(1) = V/C(2)
      DO 10 I=3,N2
        I1 = I - 1
        I2 = I - 2
        U = A(I2)/C(I2)
        V = B(I1)/C(I1) - B(I2)*U
        C(I) = SQRT(C(I)-U*U-V*V)
        A(I2) = U/C(I)
        B(I1) = V/C(I)
   10 CONTINUE
      G(1) = A(N1)/C(1)
      G(2) = -B(1)*G(1)
      U = A(N3)/C(N3)
      V = B(N2)/C(N2) - B(N3)*U
      DO 20 I=3,N2
        G(I) = -B(I-1)*G(I-1) - A(I-2)*G(I-2)
   20 CONTINUE
      W = 0.0
      DO 30 I=1,N4
        W = W + G(I)*G(I)
   30 CONTINUE
      C(N1) = SQRT(C(N1)-W-(G(N3)+U)**2-(G(N2)+V)**2)
      W = 1.0/C(N1)
      A(N3) = U*W
      B(N2) = V*W
      DO 40 I=1,N2
        G(I) = G(I)*W
   40 CONTINUE
      H(1) = B(N)/C(1)
      H(2) = A(N)/C(2) - B(1)*H(1)
      U = A(N2)/C(N2)
      V = B(N1)/C(N1) - B(N2)*U
      DO 50 I=3,N1
        H(I) = -B(I-1)*H(I-1) - A(I-2)*H(I-2)
   50 CONTINUE
      W = 0.0
      DO 60 I=1,N2
        W = W + G(I)*H(I)
   60 CONTINUE
      H(N1) = H(N1) - W - G(N2)*U
      W = 0.0
      DO 70 I=1,N3
        W = W + H(I)*H(I)
   70 CONTINUE
      C(N) = SQRT(C(N)-W-(H(N2)+U)**2-(H(N1)+V)**2)
      W = 1.0/C(N)
      A(N2) = U*W
      B(N1) = V*W
      DO 80 I=1,N1
        H(I) = H(I)*W
   80 CONTINUE
      RETURN
      END
      SUBROUTINE SOLVE(E, F, G, H, Q, N)                                SOL   10
C THE SUBROUTINE SOLVES THE SET OF N LINEAR EQUATIONS RTZ=Q WHERE T IS
C AN UPPER TRIANGULAR MATRIX WITH UNIT ELEMENTS ON THE DIAGONAL AND R
C DENOTES THE TRANSPOSE OF T. THE OTHER NON-ZERO ELEMENTS OF T OCCUR ON
C THE SUPER-DIAGONAL I.E. T(I,I+1), THE SUPER-SUPER-DIAGONAL I.E.
C T(I,I+2) AND IN THE LAST TWO COLUMNS I.E. T(I,N-1) AND T(I,N).
C THE STORAGE IS SUCH THAT THE FOLLOWING ARE EQUIVALENT----
C        T(I,I+1)   AND E(I)   FOR I=1 TO N-3
C        T(I,I+2)   AND F(I)   FOR I=1 TO N-4
C        T(I,N-1)   AND G(I)   FOR I=1 TO N-4
C        T(I,N)     AND H(I)   FOR I=1 TO N-3
C        T(N-3,N-1) AND F(N-3)+G(N-3)
C        T(N-2,N-1) AND E(N-2)+G(N-2)
C        T(N-2,N)   AND F(N-2)+H(N-2)
C        T(N-1,N)   AND E(N-1)+H(N-1)
C THE SOLUTION IS EFFECTED BY A FORWARD SUBSTITUTION PROCESS FOLLOWED
C BY A BACKWARD SUBSTITUTION AND THE INPUT VECTOR Q IS SUCCESSIVELY
C OVERWRITTEN BY THE INTERMEDIATE SOLUTION (OF THE FORWARD ELIMINATION)
C AND THEN THE FINAL SOLUTION.
      DIMENSION E(N), F(N), G(N), H(N), Q(N)
      N1 = N - 1
      N2 = N - 2
      N3 = N - 3
      Q(2) = Q(2) - E(1)*Q(1)
      DO 10 I=3,N2
        Q(I) = Q(I) - E(I-1)*Q(I-1) - F(I-2)*Q(I-2)
   10 CONTINUE
      U = 0.0
      V = 0.0
      DO 20 I=1,N2
        U = U + G(I)*Q(I)
        V = V + H(I)*Q(I)
   20 CONTINUE
      Q(N1) = Q(N1) - E(N2)*Q(N2) - F(N3)*Q(N3) - U
      Q(N) = Q(N) - (E(N1)+H(N1))*Q(N1) - F(N2)*Q(N2) - V
      Q(N1) = Q(N1) - (E(N1)+H(N1))*Q(N)
      Q(N2) = Q(N2) - (E(N2)+G(N2))*Q(N1) - (F(N2)+H(N2))*Q(N)
      DO 30 II=1,N3
        I = N3 - II + 1
        Q(I) = Q(I) - E(I)*Q(I+1) - F(I)*Q(I+2) - G(I)*Q(N1) -
     *   H(I)*Q(N)
   30 CONTINUE
      RETURN
      END
      SUBROUTINE RHS(D, S, N)                                           RHS   10
C THE SUBROUTINE FORMS Q WHERE BQ=S AND B IS A DIAGONAL MATRIX OF ORDER
C N WHOSE NON-ZERO DIAGONAL ENTRIES ARE STORED IN A VECTOR D OF N
C ELEMENTS. THE INPUT VECTOR S IS OVERWRITTEN BY THE RESULT Q
      DIMENSION D(N), S(N)
      DO 10 I=1,N
        S(I) = S(I)/D(I)
   10 CONTINUE
      RETURN
      END