GMFADS Subroutine

      SUBROUTINE GMFADS(NN,A,NWK,MAXA)
C
C     This subroutine computes the LDU decomposition of a symmetric positive
C     definite matrix B where only the upper triangular skyline structure
C     is stored.  The decomposition is done by the Gill-Murray
C     strategy from P.E. Gill and W. Murray, Newton type Methods
C     for Unconstrained and Linearly Constrained Optimization,
C     Mathematical Programming, 7, 311-350 (1974) and gives an
C     approximate decomposition in the case of a nonpositive
C     definite or ill-conditioned matrix.
C
C     Input variables:
C
C        NN -- dimension of B.
C
C        A -- one dimensional real array containing the upper 
C             triangular skyline portion of a symmetric matrix B in
C             packed skyline storage format.
C
C        NWK -- number of elements in A.
C
C        MAXA -- an integer array of dimension NN+1 containing the 
C                locations of the diagonal elements of B in A.  
C                By convention, MAXA(NN+1)=NWK+1.  
C
C     Output variables:
C
C        A -- the upper triangular skyline portion of the LDU 
C             decomposition of the symmetric matrix B (or B + E if B
C             was not sufficiently positive definite).
C
C
C     No working storage is required by this routine.
C
C
      USE REAL_PRECISION
      INTEGER, INTENT(IN):: NN,NWK,MAXA(NN+1)
      REAL (KIND=R8), INTENT(IN OUT):: A(NWK)
      INTEGER:: I,I0,I1,I2,I3,I4,J,J1,K,K1,K2,KH,KL,KN,KU,KZ,L,L1,
     &   L2,L3,M,M1,N1,NNN
      REAL (KIND=R8):: BET,DEL,DJ,G,GAM,GAM1,PHI,
     &   THE,THE1,XT1,XT2,ZET,ZET1
      G=0.0
      GAM=0.0
      DO I=1,NN
         K=MAXA(I)
         G=G+A(K)*A(K)
         GAM1=ABS(A(K))
         IF(GAM1.GT.GAM)GAM=GAM1
      END DO
      ZET=0.0
      DO I=1,NN
         K=MAXA(I)
         K1=MAXA(I+1)-1
         K2=K1-K
         IF (K2.EQ.0) CYCLE
         L=K+1
         DO J=L,K1
            G=G+2.0*A(J)*A(J)
            ZET1=ABS(A(J))
            IF(ZET1.GT.ZET)ZET=ZET1
         END DO
      END DO
      ZET=ZET/NN
      DEL=EPSILON(1.0_R8)
      BET=DEL
      IF (ZET .GT. BET) BET=ZET
      IF (GAM .GT. BET) BET=GAM
      G=SQRT(G)
      IF (G .GT. 1.0) DEL=DEL*G
      DO I=1,NN
         N1=I-1
         KN=MAXA(I)
         KL=KN+1
         KU=MAXA(I+1)-1
         KH=KU-KL
         PHI=A(KN)
         IF (KH .GE. 0) THEN
           K1=KN+1
           K2=I
           DO J=K1,KU
              K2=K2-1
              KZ=MAXA(K2)
              PHI=PHI-A(J)*A(J)*A(KZ)
           END DO
         END IF
         PHI=ABS(PHI)
         L=I+1
         THE=0.0
         NNN=NN+1
         IF (L .NE. NNN) THEN
           DO J=L,NN
              L1=MAXA(J)
              L2=MAXA(J+1)
              L3=L2-L1-1
              M=J-I
              IF (L3 .LT. M) CYCLE
              M1=L1+M
              IF (N1 .NE. 0) THEN
                DO J1=1,N1
                  I0=MAXA(J1)
                  I1=MAXA(L)
                  I2=I-J1
                  I3=I1-KN-1
                  I4=J-J1
                  IF (I3 .LT. I2) CYCLE
                  IF (L3 .LT. I4) CYCLE
                  XT1=A(KN+I2)
                  XT2=A(L1+I4)
                  A(M1)=A(M1)-XT1*XT2*A(I0)
                END DO
              END IF
            THE1=ABS(A(M1))
            IF (THE .LT. THE1) THE=THE1
            END DO
         END IF
         THE=THE*THE/BET
         DJ=DEL
         IF (PHI .GT. DJ) DJ=PHI
         IF (THE .GT. DJ) DJ=THE
         A(KN)=DJ
         IF (L .EQ. NNN) CYCLE
         DO J=L,NN
            L1=MAXA(J)
            L2=MAXA(J+1)
            L3=L2-L1-1
            M=J-I
            IF (L3 .LT. M) CYCLE
            M1=L1+M
            A(M1)=A(M1)/A(KN)
         END DO
      END DO
      RETURN
      END SUBROUTINE GMFADS