C.hr FACTOR
C@
C....*...1.........2.........3.........4.........5.........6.........7.*
C     FACTOR      12/20/87 01/01/01
C
C     PURPOSE
C     FACTOR A SYMMETRIC POSITIVE DEFINITE MATRIX A = RR' WHERE R IS
C     UPPER TRIANGULAR
C
C     USAGE
C     CALL FACTOR(A,DEL,M,EPS,IER)
C
C     ARGUMENTS
C     A   - ON INPUT CONTAINS THE UPPER TRIANGULAR PART OF A STORED
C           COLUMNWISE (STORAGE MODE 1).
C           ON OUTPUT CONTAINS THE UPPER TRIANGULAR PART OF R STORED
C           COLUMNWISE (STORAGE MODE 1).
C           THE (I,J) ELEMENT OF A IS IN PHYSICAL LOCATION J*(J-1)/2+I
C           WITH I.LE.J; THE TOTAL LENGTH IS M*(M+1)/2. THE SAME FOR R.
C           REAL*8,
C     DEL - ON OUTPUT CONTAINS THE JACOBIAN OF R WITH RESPECT TO A, AN
C           M*(M+1)/2 BY M*(M+1)/2 MATRIX STORED COLUMNWISE (STORAGE
C           MODE 0).  THE DERIVATIVE OF THE (I,J) ELEMENT OF R WITH
C           RESPECT TO THE (K,L) ELEMENT OF A WITH I.LE.J AND K.LE.L
C           IS IN CONCEPTUAL LOCATION (J*(J-1)/2+I,L*(L-1)/2+K) AND IN
C           PHYSICAL LOCATION (M*(M+1)/2)*(L*(L-1)/2+K-1)+J*(J-1)/2+I)
C           REAL*8
C     M   - INPUT, THE NUMBER OF ROWS (COLUMNS) IN A.
C           INTEGER*4
C     EPS - INPUT CONSTANT BETWEEN ZERO AND ONE WHICH IS USED AS
C           RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE, A
C           REASONABLE VALUE IS 1.D-13.
C           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
C     NOTE
C     THE IMPLICIT LOWER PART OF A IS THE TRANSPOSE OF THE UPPER PART.
C     THE IMPLICIT LOWER PART OF R IS ZEROES BELOW THE DIAGONAL.
C
      SUBROUTINE FACTOR(A,DEL,M,EPS,IER)
      IMPLICIT REAL*8 (A-H,O-Z)
      save
      REAL*8 A((M*M+M)/2),DEL((M*M+M)/2,(M*M+M)/2)
C
      IF (M.LE.0) THEN
      IER=-1
      RETURN
      ENDIF
C
      IF ((EPS.LE.0).OR.(EPS.GE.1.D0)) THEN
      IER=-1
      RETURN
      ENDIF
C
      LA=(M*M+M)/2
C
C     PUT DEL TO THE IDENTITY
C
      DO 10 I=1,LA
      DO 10 J=1,LA
10    DEL(I,J)=0.D0
      DO 20 I=1,LA
20    DEL(I,I)=1.D0
C
C     FACTORIZE; THIS IS A CHOLESKY FACTORIZATION GOTTEN BY WORKING FROM
C     SOUTHEAST TO NORTHWEST TO GET AN UPPER TRIANGULAR MATRIX INSTEAD
C     OF THE TRADITIONAL LOWER TRIANGULAR MATRIX GOTTEN BY WORKING FROM
C     NORTHWEST TO SOUTHEAST.
C
      IER=0
      DO 100 J1=0,M-1
      J=M-J1
      JJ=(J*J+J)/2
C
C     TAKE THE SQUARE ROOT OF AJJ, PUT IT IN AJJ
C
      R=A(JJ)
      IF (R.GT.0.D0) THEN
      R=DSQRT(R)
      A(JJ)=R
      DO 30 L=1,M
      DO 30 K=1,L
      KL=(L*L-L)/2+K
30    DEL(JJ,KL)=(0.5D0/R)*DEL(JJ,KL)
      ELSE
      IER=-1
      RETURN
      ENDIF
C
C     WE ARE FINISHED IF JJ.EQ.1, EQUIVALENTLY IF J.EQ.1
C
      IF (JJ.EQ.1) RETURN
C
C     DIVIDE THE ELEMENTS OF THE COLUMN ABOVE AJJ BY AJJ
C
      DO 50 I=1,J-1
      A((J*J-J)/2+I)=A((J*J-J)/2+I)/R
      DO 40 L=1,M
      DO 40 K=1,L
      KL=(L*L-L)/2+K
40    DEL((J*J-J)/2+I,KL)=DEL((J*J-J)/2+I,KL)/R
     &                    -(A((J*J-J)/2+I)/R)*DEL(JJ,KL)
50    CONTINUE
C
C     SWEEP THE NORTHEAST SUBMATRIX
C
      ASAVE=A((J*J-J)/2)
      DO 70 J0=1,J-1
      DO 70 I0=1,J0
      IJ=(J0*J0-J0)/2+I0
      A(IJ)=A(IJ)-A((J*J-J)/2+I0)*A((J*J-J)/2+J0)
      DO 60 L=1,M
      DO 60 K=1,L
      KL=(L*L-L)/2+K
60    DEL(IJ,KL)=DEL(IJ,KL)
     &           -A((J*J-J)/2+I0)*DEL((J*J-J)/2+J0,KL)
     &           -DEL((J*J-J)/2+I0,KL)*A((J*J-J)/2+J0)
70    CONTINUE
C
C     TOLERANCE CHECK
C
      IF(ASAVE.LE.0.0) THEN
      IER=-1
      RETURN
      ENDIF
      IF(A((J*J-J)/2).LT.DABS(EPS*ASAVE)) IER=J1
C
100   CONTINUE
C
      RETURN
      END