Numérotation des lignes :

choldc
C CHOLDC    SOURCE    CHAT      06/03/29    21:16:29     5360      SUBROUTINE CHOLDC(NP,N,A,P,INVA,INVP,iarr)C************************************************************************CC PROJET            :  CASTEM 2000CC NOM               :  CHOLDCCC DESCRIPTION       :  Appelle par RLECA1, RLENCF, RLENCO, RLEXCA.CC LANGAGE           :  FORTRAN 77CC AUTEUR            :  A. BECCANTINICC************************************************************************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     

© Cast3M 2003 - Tous droits réservés.
Mentions légales