C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMCHL      1/16/89 (Drew's twelfth birthday)
C
C     PURPOSE
C     FACTOR A SYMMETRIC, POSITIVE DEFINITE MATRIX AS A = RR' AND FACTOR
C     ITS INVERSE AS INV(A)=P'P WHERE R AND P ARE LOWER TRIANGULAR
C
C     USAGE
C     CALL DGMCHL(A,M,R,P,EPS,IER)
C
C     SUBROUTINES CALLED
C     DMFSD, DRINV
C
C     ARGUMENTS
C     A   - AN M BY M SYMMETRIC, POSITIVE DEFINITE MATRIX STORED
C           COLUMNWISE (STORAGE MODE 0).
C           INPUT, REAL*8
C     R   - FACTOR OF A, A=RR', AN M BY M MATRIX STORED COLUMNWISE
C           (STORAGE MODE 0).
C           OUTPUT, REAL*8
C     P   - FACTOR OF A INVERSE, INV(A)=P'P, AN M BY M MATRIX STORED
C           COLUMNWISE (STORAGE MODE 0).
C           OUTPUT, REAL*8
C     M   - THE NUMBER OF ROWS (COLUMNS) IN A.
C           INPUT, INTEGER*4
C     EPS - CONSTANT BETWEEN ZERO AND ONE WHICH IS USED AS RELATIVE
C           TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE, A REASONABLE
C           VALUE IS 1.D-13.
C           INPUT, REAL*8
C     IER - OUTPUT ERROR PARAMETER, CODED AS FOLLOWS:
C           IER=0  - NO ERROR.
C           IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAMETER M,EPS OR
C                    BECAUSE SOME RADICAND IS NOT POSITIVE (MATRIX A IS
C                    NOT POSITIVE DEFINITE, POSSIBLY DUE TO LOSS OF
C                    SIGNIFICANCE).
C           IER=J  - WARNING WHICH INDICATES LOSS OF SIGNIFICANCE.  THE
C                    RADICAND FORMED AT FACTORIZATION STEP J+1 WAS STILL
C                    POSITIVE BUT NO LONGER GREATER THAN
C                    DABS(EPS*A(J*(J+1)/2)).
C
      SUBROUTINE DGMCHL(A,M,R,P,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 A(1),R(1),P(1)
C
C     COMPUTE R' OF THE CHOLESKY FACTORIZATION, PUT IT IN P
C
      DO 10 J=1,M
      DO 10 I=1,J
10    P(J*(J-1)/2+I)=A(M*(J-1)+I)
C
      CALL DMFSD(P,M,EPS,IER)
      IF (IER.LT.0) RETURN
C
C     TRANSPOSE AND PUT THE FACTOR IN R IN STORAGE MODE 0
C
      DO 20 J=1,M
      DO 20 I=1,J
20    R(M*(J-1)+I)=0.D0
C
      DO 30 J=1,M
      DO 30 I=1,J
30    R(M*(I-1)+J)=P(J*(J-1)/2+I)
C
C     INVERT R' TO GET P' AND PUT IT IN P
C
      CALL DRINV(P,M)
C
C     EXPAND P' TO STORAGE MODE 0
C
      DO 40 JJ=1,M
      J=M+1-JJ
      DO 40 II=JJ,M
      I=M+1-II
40    P(M*(J-1)+I)=P(J*(J-1)/2+I)
C
C     FILL IN LOWER TRIANGLE
C
      DO 50 J=1,M
      DO 50 I=1,J
50    P(M*(I-1)+J)=P(M*(J-1)+I)
C
C     ZERO OUT UPPER TRIANGLE
C
      DO 60 J=1,M
      DO 60 I=1,J
      IF(I.NE.J) THEN
         P(M*(J-1)+I)=0.D0
         ENDIF
60    CONTINUE
C
      RETURN
      END