C.hr DGMFAC
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     DGMFAC      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 UPPER TRIANGULAR
C
C     USAGE
C     CALL DGMFAC(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 DGMFAC(A,M,R,P,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 A(M*M),R(M*M),P(M*M)
C
C     PERMUTE THE ROWS AND COLUMNS OF A AND STORE IN P IN STORAGE MODE 1
C     THE PERMUTAIONS ARE BECAUSE WE FACTOR AS UPPER TRAINGULAR INSTEAD
C     OF THE STANDARD LOWER TRIANGULAR CHOLESKY FACTORIZATION.
C
      DO 10 J=1,M
      JJ=M+1-J
      DO 10 I=J,M
      II=M+1-I
10    P(JJ*(JJ-1)/2+II)=A(M*(J-1)+I)
C
C     COMPUTE THE CHOLESKY FACTORIZATION
C
      CALL DMFSD(P,M,EPS,IER)
      IF (IER.LT.0) RETURN
C
C     UNDO THE PERMUTATIONS AND PUT THE FACTOR IN R IN STORAGE MODE 1
C
      DO 20 J=1,M
      JJ=M+1-J
      DO 20 I=1,J
      II=M+1-I
20    R(J*(J-1)/2+I)=P(II*(II-1)/2+JJ)
C
C     COPY R TO P AND INVERT P
C
      DO 30 I=1,M*(M+1)/2
30    P(I)=R(I)
      CALL DRINV(P,M)
C
C     R AND P ARE IN COMPRESSED FORM, STORAGE MODE 1.  EXPAND TO
C     STORAGE MODE 0
C
      DO 40 JJ=1,M
      J=M+1-JJ
      DO 40 II=JJ,M
      I=M+1-II
      R(M*(J-1)+I)=R(J*(J-1)/2+I)
40    P(M*(J-1)+I)=P(J*(J-1)/2+I)
C
      DO 50 J=1,M
      DO 50 I=J,M
      IF(I.NE.J) THEN
         R(M*(J-1)+I)=0.D0
         P(M*(J-1)+I)=0.D0
         ENDIF
50    CONTINUE
C
      RETURN
      END