Esope source files

   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 :