C CHOLDC    SOURCE    PV090527  25/11/25    21:15:01     12407          
      SUBROUTINE CHOLDC(NP,N,A,P,INVA,INVP,iarr)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CHOLDC
C
C DESCRIPTION       :  Appelle par RLECA1, RLENCF, RLENCO, RLEXCA.
C
C LANGAGE           :  FORTRAN 77
C
C AUTEUR            :  A. BECCANTINI
C
C************************************************************************
C
*
***** CHOLESKY DECOMPOSITION OF A: A=L.TRANPOSE(L)
*     (FROM NUMERICAL RECIPIES)
*
*     INPUT
*     A = SYMMETRIC AND POSITIVE DEFINITE MATRIX (ONLY THE UPPER
*         TRIANGLE IS NEEDED)
*     NP = PHYSICAL DIMENSION OF A
*     N  = DIMENSION OF THE SUBMATRIX TO DECOMPOSE
*
*     OUTPUT
*     A = (LOWER TRIANGLE OF A CONTAIN L, EXCEPT THE DIAGONAL
*     P = IT STORES THE DIAGONAL OF L
*     INVA = UPPER TRIANGLE CONTAIN INV(A),
*            LOWER TRIANGLE CONTAIN INV(L), EXCEPT THE DIAGONAL
*     INVP = IT STORES THE DIAGONAL OF INV(L)
*     iarr = 0 NO PROBLEMS
*
      IMPLICIT INTEGER(I-N)
      INTEGER NP,N, iarr
      REAL*8 A(NP,NP),P(NP),INVA(NP,NP),INVP(NP)
*
      INTEGER I, J, K
      REAL*8 SUM
      iarr = 0
      DO I = 1, N, 1
         DO J = I, N, 1
            SUM = A(I,J)
            DO K = I - 1, 1, -1
               SUM = SUM - A(I,K) * A (J,K)
            ENDDO
            IF ( I .EQ. J)THEN
               IF(SUM .LE. 0.00)THEN
                  iarr = 1
                  GOTO 9999
               ELSE
                  P(I) = SQRT(SUM)
               ENDIF
            ELSE
               A(J,I) = SUM / P(I)
            ENDIF
         ENDDO
      ENDDO
*
***** WE COMPUTE INV(L) IN WE STORE IT INTO INVP AND ON THE
*     LOWER PART OF INVA
*
      DO J = 1, N, 1
         INVP(J) = 1.0D0 / P(J)
         DO I = J + 1, N, 1
            SUM = -1.0D0 * INVP(J) * A(I,J)
            DO K = J+1, I-1, 1
               SUM = SUM - A(I,K) * INVA(K,J)
            ENDDO
            INVA(I,J) = SUM / P(I)
         ENDDO
      ENDDO
*
**** NOW WE COMPUTE INV(A)= TRANPOSE(INV(L)) . INV(L)
*    AND WE STORE IT IN THE UPPER TRIANGLE OF INVA
*
      DO J = N, 1, -1
         SUM = INVP(J)**2
         DO K = J+1,N,1
            SUM = SUM + INVA(K,J)**2
         ENDDO
         INVA(J,J) = SUM
         DO I = J-1, 1, -1
            SUM = INVA(J,I) * INVP(J)
            DO K = J+1,N,1
               SUM = SUM + INVA(K,I) * INVA(K,J)
            ENDDO
            INVA(I,J) = SUM
         ENDDO
      ENDDO
*
 9999 CONTINUE
      RETURN
      END