Télécharger la plaquette Cast3M


fr en

Télécharger choldc.eso

Retour à la liste

Numérotation des lignes : 

choldc

  1. C CHOLDC    SOURCE    CHAT      06/03/29    21:16:29     5360
  2.       SUBROUTINE CHOLDC(NP,N,A,P,INVA,INVP,iarr)
  3. C************************************************************************
  4. C
  5. C PROJET            :  CASTEM 2000
  6. C
  7. C NOM               :  CHOLDC
  8. C
  9. C DESCRIPTION       :  Appelle par RLECA1, RLENCF, RLENCO, RLEXCA.
  10. C
  11. C LANGAGE           :  FORTRAN 77
  12. C
  13. C AUTEUR            :  A. BECCANTINI
  14. C
  15. C************************************************************************
  16. C
  17. *
  18. ***** CHOLESKY DECOMPOSITION OF A: A=L.TRANPOSE(L)
  19. *     (FROM NUMERICAL RECIPIES)
  20. *
  21. *     INPUT
  22. *     A = SYMMETRIC AND POSITIVE DEFINITE MATRIX (ONLY THE UPPER
  23. *         TRIANGLE IS NEEDED)
  24. *     NP = PHYSICAL DIMENSION OF A
  25. *     N  = DIMENSION OF THE SUBMATRIX TO DECOMPOSE
  26. *
  27. *     OUTPUT
  28. *     A = (LOWER TRIANGLE OF A CONTAIN L, EXCEPT THE DIAGONAL
  29. *     P = IT STORES THE DIAGONAL OF L
  30. *     INVA = UPPER TRIANGLE CONTAIN INV(A),
  31. *            LOWER TRIANGLE CONTAIN INV(L), EXCEPT THE DIAGONAL
  32. *     INVP = IT STORES THE DIAGONAL OF INV(L)
  33. *     iarr = 0 NO PROBLEMS
  34. *
  35.       IMPLICIT INTEGER(I-N)
  36.       INTEGER NP,N, iarr
  37.       REAL*8 A(NP,NP),P(NP),INVA(NP,NP),INVP(NP)
  38. *
  39.       INTEGER I, J, K
  40.       REAL*8 SUM
  41.       iarr = 0
  42.       DO I = 1, N, 1
  43.          DO J = I, N, 1
  44.             SUM = A(I,J)
  45.             DO K = I - 1, 1, -1
  46.                SUM = SUM - A(I,K) * A (J,K)
  47.             ENDDO
  48.             IF ( I .EQ. J)THEN
  49.                IF(SUM .LE. 0.00)THEN
  50.                   iarr = 1
  51.                   GOTO 9999
  52.                ELSE
  53.                   P(I) = SQRT(SUM)
  54.                ENDIF
  55.             ELSE
  56.                A(J,I) = SUM / P(I)
  57.             ENDIF
  58.          ENDDO
  59.       ENDDO
  60. *
  61. ***** WE COMPUTE INV(L) IN WE STORE IT INTO INVP AND ON THE
  62. *     LOWER PART OF INVA
  63. *
  64.       DO J = 1, N, 1
  65.          INVP(J) = 1.0D0 / P(J)
  66.          DO I = J + 1, N, 1
  67.             SUM = -1.0D0 * INVP(J) * A(I,J)
  68.             DO K = J+1, I-1, 1
  69.                SUM = SUM - A(I,K) * INVA(K,J)
  70.             ENDDO
  71.             INVA(I,J) = SUM / P(I)
  72.          ENDDO
  73.       ENDDO
  74. *
  75. **** NOW WE COMPUTE INV(A)= TRANPOSE(INV(L)) . INV(L)
  76. *    AND WE STORE IT IN THE UPPER TRIANGLE OF INVA
  77. *
  78.       DO J = N, 1, -1
  79.          SUM = INVP(J)**2
  80.          DO K = J+1,N,1
  81.             SUM = SUM + INVA(K,J)**2
  82.          ENDDO
  83.          INVA(J,J) = SUM
  84.          DO I = J-1, 1, -1
  85.             SUM = INVA(J,I) * INVP(J)
  86.             DO K = J+1,N,1
  87.                SUM = SUM + INVA(K,I) * INVA(K,J)
  88.             ENDDO
  89.             INVA(I,J) = SUM
  90.          ENDDO
  91.       ENDDO
  92. *
  93.  9999 CONTINUE
  94.       RETURN
  95.       END
  96.  
  97.  
  98.  
  99.  
  100.  

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