Télécharger choldc.eso

Retour à la liste

Numérotation des lignes :

  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