Sources Esope | Cast3M

inver1
C INVER1    SOURCE    CHAT      05/01/13    00:42:36     5004
      SUBROUTINE INVER1(TRAVAI,NZ,ICRIT,EPS)
C
C  ====================================================================
C    SOUS-PROGRAMME APPELE PAR RAYE3 (RAYONNEMENT, oper. RAYN)
C    -->  derive de inver.eso:
C    -->  version esope
C    -->  on travaille sur les matrices transposees
C    -->  traitement des tableaux par colonne
C    RAPPEL:
C    INVERSION DE MATRICE PLEINE NON SYMETRIQUE PAR RESOLUTION
C    SUCCESSIVE DE NZ SYSTEMES LINEAIRES
C    A  MATRICE (NZ*NZ) A INVERSER EN ENTREE, MATRICE INVERSEE EN SORTIE
C    ICRIT=1 SI MATRICE SINGULIERE, 0 SINON
C    BT  TABLEAU DE REELS DE DIMENSION NZ*NZ
C    IS  TABLEAU D'ENTIERS DE DIMENSION NZ
C    EPS PRECISION
C
C    Le segment TRAVAI contenant A et BT est envoye actif par RAYE3
C    et retourné actif
C
C    var locales : BTT transposee de BT   dans segment DIVERS
C  ====================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC PPARAM
-INC CCOPTIO
 
      SEGMENT TRAVAI
         REAL *8 A(NEL,NEL), BT(NEL,NEL)
         INTEGER IS(NEL)
      ENDSEGMENT
 
 
      SEGMENT DIVERS
         REAL*8 BTT(NZ,NZ)
      ENDSEGMENT
      SEGINI DIVERS
C
C     on peut stocker ABS(A) au lieu de le recalculer dans la boucle
C     40 mais d'une part le gain en temps calcul n'est pas substantiel
C     et d'autre part cela implique une place memoire supplementaire
C
C     DO 160 J=1,NZ
C        DO 160 I=1,NZ
C        ABA(I,J)=ABS(A(I,J))
C 160 CONTINUE
C
C      INITIALISATIONS
C      IS SUITE REPRESENTANT L'INDICE J DE X(J) SOLU CORRESPONDANT
C      A LA I EME COLONNE DE LA MATRICE A TRIANGULARISEE
C
      DO 20 I=1,NZ
         DO 10 J=1,NZ
         BT(J,I)=0.D0
  10     CONTINUE
      BT(I,I)=1.D0
      IS(I)=I
  20  CONTINUE
      ICRIT=0
C
C  1-  TRIANGULARISATION
C
      NZM1=NZ-1
      DO 100 NR=1,NZM1
C
C      CHOIX DU PIVOT
C
         PIVOT=0.D0
         DO 40 K=NR,NZ
            DO 40 L=NR,NZ
C              ABSKL=ABA(L,K)
               ABSKL=ABS(A(L,K))
               IF(ABSKL.GT.PIVOT) THEN
                 I=K
                 J=L
                 PIVOT=ABSKL
               ENDIF
  40     CONTINUE
C
C     LE PIVOT EST-IL NUL?
C
         IF(PIVOT.LE.EPS) THEN
           ICRIT=1
           RETURN
         ENDIF
C
C      CHANGEMENT DE LIGNE : PLACE LE PIVOT EN NR EME LIGNE
C
         DO 50 L=1,NZ
            D=A(L,NR)
            A(L,NR)=A(L,I)
            A(L,I)=D
   50    CONTINUE
         DO 60 L=1,NZ
            E=BT(L,NR)
            BT(L,NR)=BT(L,I)
            BT(L,I)=E
   60    CONTINUE
C
C      CHANGEMENT DE COLONNE : PLACE LE PIVOT EN R EME COLONNE
C
         DO 70 M=1,NZ
            CC=A(NR,M)
            A(NR,M)=A(J,M)
            A(J,M)=CC
   70    CONTINUE
C
C      INDICE DES VARIABLES CORRESPONDANT A LA J EME ET A LA R EME COLON
C
         ISR=IS(NR)
         IS(NR)=IS(J)
         IS(J)=ISR
C
C      CALCUL DE LA NOUVELLE MATRICE A
C
         NRP1=NR+1
         DO 90 I=NRP1,NZ
            IF(A(NR,I).NE.0.D0)THEN
               G=A(NR,I)/A(NR,NR)
               DO 80 J=1,NZ
                  A(J,I)=A(J,I)-G*A(J,NR)
                  BT(J,I)=BT(J,I)-G*BT(J,NR)
   80          CONTINUE
            ENDIF
   90    CONTINUE
  100 CONTINUE
C
C  2-  RESOLUTION DU SYSTEME TRIANGULARISE
C
* On introduit la transposée de B pour performance calcul fortran
*
      DO 150 I=1,NZ
         DO 150 J=1,NZ
         BTT(J,I)=BT(I,J)
  150 CONTINUE
 
      DO 130 J=1,NZ
      BTT(NZ,J)=BTT(NZ,J)/A(NZ,NZ)
         DO 120 I= NZM1,1,-1
         F=0.D0
         I1=I+1
            DO 110 JJ=I1,NZ
            F=F-A(JJ,I)*BTT(JJ,J)
  110       CONTINUE
         BTT(I,J)=(BTT(I,J)+F)/A(I,I)
  120    CONTINUE
  130 CONTINUE
C
      DO 140 L=1,NZ
      IL=IS(L)
      DO 140 J=1,NZ
      A(J,IL)=BTT(L,J)
  140 CONTINUE
 
      SEGSUP DIVERS
      RETURN
      END