SUBROUTINE CHOLD(A,IRDA,NA,DETA,*)                                    0420
C                                                                           0430
C MATRIX INVERSION USING CHOLESKI DECOMPOSITION                             0440
C                                                                           0450
C INPUT ARGUMENTS                                                           0460
C A = ARRAY CONTAINING POSITIVE DEFINITE SYMMETRIC INPUT MATRIX             0470
C IRDA = ROW DIMENSION OF ARRAY CONTAINING INPUT MATRIX                     0480
C NA = SIZE OF INPUT MATRIX                                                 0490
C OUTPUT ARGUMENTS                                                          0500
C A = CONTAINS INVERSE OF INPUT MATRIX (INPUT DESTROYED)                    0510
C DETA = DETERMINANT OF INPUT MATRIX                                        0520
C * = ERROR RETURN (TAKEN IF NA .LT. 1 OR IF DETA .LT. SING)                0530
C                                                                           0540
      DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING                  0550
      DIMENSION A(IRDA,NA)                                                  0560
      SQRT(SUM) = DSQRT(SUM)                                                0570
      ABS(DETA) = DABS(DETA)                                                0580
      DATA SING/1D-10/                                                      0590
C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX             0600
      IF(NA .LT. 1) GO TO 18                                                0610
      DETA = A(1,1)                                                         0620
      A(1,1) = SQRT(A(1,1))                                                 0630
      IF(NA .EQ. 1) GO TO 6                                                 0640
      DO 1 I = 2,NA                                                         0650
    1   A(I,1) = A(I,1) / A(1,1)                                            0660
      DO 5 J = 2,NA                                                         0670
        SUM = 0.                                                            0680
        J1 = J - 1                                                          0690
        DO 2 K = 1,J1                                                       0700
    2     SUM = SUM + A(J,K) ** 2                                           0710
        DETA = DETA * (A(J,J) - SUM)                                        0720
        A(J,J) = SQRT(A(J,J) - SUM)                                         0730
        IF(J .EQ. NA) GO TO 5                                               0740
        J2 = J + 1                                                          0750
        DO 4 I = J2,NA                                                      0760
          SUM = 0.                                                          0770
          DO 3 K = 1,J1                                                     0780
    3       SUM = SUM + A(I,K) * A(J,K)                                     0790
    4     A(I,J) = (A(I,J) - SUM) / A(J,J)                                  0800
    5   CONTINUE                                                            0810
    6 IF(ABS(DETA) .LT. SING) GO TO 16                                      0820
C INVERSION OF LOWER TRIANGULAR MATRIX                                      0830
      DO 7 I = 1,NA                                                         0840
    7   A(I,I) = 1. / A(I,I)                                                0850
      IF(NA .EQ. 1) GO TO 10                                                0860
      N1 = NA - 1                                                           0870
      DO 9 J = 1,N1                                                         0880
        J2 = J + 1                                                          0890
        DO 9 I = J2,NA                                                      0900
          SUM = 0.                                                          0910
          I1 = I - 1                                                        0920
          DO 8 K = J,I1                                                     0930
    8       SUM = SUM + A(I,K) * A(K,J)                                     0940
    9     A(I,J) = - A(I,I) * SUM                                           0950
C CONSTRUCTION OF INVERSE OF INPUT MATRIX                                   0960
   10 DO 15 J = 1,NA                                                        0970
        IF(J .EQ. 1) GO TO 12                                               0980
        J1 = J - 1                                                          0990
        DO 11 I = 1,J1                                                      1000
   11     A(I,J) = A(J,I)                                                   1010
   12   DO 14 I = J,NA                                                      1020
          SUM = 0.                                                          1030
          DO 13 K = I,NA                                                    1040
   13       SUM = SUM + A(K,I) * A(K,J)                                     1050
   14     A(I,J) = SUM                                                      1060
   15   CONTINUE                                                            1070
      RETURN                                                                1080
   16 WRITE(6,17) DETA                                                      1090
   17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5)                  1100
      RETURN 1                                                              1110
   18 WRITE(6,19)                                                           1120
   19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD')                       1130
      RETURN 1                                                              1140
      END                                                                   1150