Sources Esope | Cast3M

invhok
C INVHOK    SOURCE    CHAT      05/01/13    00:43:00     5004
      SUBROUTINE INVHOK(G,GINV,N)
*
*    g est une matrice carre de dimension n
*    invg est le resultat la mtruice inverse de dimension n
*
*    calcul de l'inverse de g matrice symmetrique definie positive
*    par la méthode de cholevsky
*
*    utilisée pour la matrice tangente de la loi vsico-plastique
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 G(6,6),GINV(6,6),P(6),B(6,6)
      DATA B        /1.D0, 0.D0, 0.D0,0.D0,0.D0,0.D0,
     &               0.D0, 1.D0, 0.D0,0.D0,0.D0,0.D0,
     &               0.D0, 0.D0, 1.D0,0.D0,0.D0,0.D0,
     &               0.D0, 0.D0, 0.D0,1.D0,0.D0,0.D0,
     &               0.D0, 0.D0, 0.D0,0.D0,1.D0,0.D0,
     &               0.D0, 0.D0, 0.D0,0.D0,0.D0,1.D0/
*  cholevsky
*  la diagonale est stockee dans p
*  la partie inférieure est stockée dans g: g est modifiée
       DO 30 I=1,N
         DO 20 J=I ,N
            SUM=G(I,J)
            DO 10 K=I-1,1,-1
              SUM=SUM-G(I,K)*G(J,K)
 10         CONTINUE
            IF ( I .EQ. J) THEN
               IF ( SUM .LE. 0.D0 ) then
                  SUM = ABS(SUM)
                  IF ( SUM .EQ. 0.D0) SUM = 1.D-20
               ENDIF
               P(I)=SQRT(SUM)
            ELSE
              G(J,I)=SUM/P(I)
            ENDIF
 20      CONTINUE
 30    CONTINUE
*
*     boucle sur les colonnes de b
       DO 80 L=1,N
* premier subsitution l.y=b
         DO 50 I=1,N
           SUM=B(L,I)
           DO 40 K=I-1,1,-1
             SUM=SUM-G(I,K)*GINV(L,K)
 40        CONTINUE
           GINV(L,I)=SUM/P(I)
 50      CONTINUE
* seconde subsitution tl.x=y
         DO 70 I=N,1,-1
           SUM=GINV(L,I)
           DO 60 K=I+1,N
             SUM=SUM-G(K,I)*GINV(L,K)
 60        CONTINUE
           GINV(L,I)=SUM/P(I)
 70      CONTINUE
 80    CONTINUE
       RETURN
       END