Sources Esope | Cast3M

inver
C INVER     SOURCE    SP204843  26/02/16    21:15:06     12477          
      SUBROUTINE INVER(A,NZ,ICRIT,B,IS,EPS)
C
C  ====================================================================
C    SOUS-PROGRAMME FORTRAN APPELE PAR ELFINV
C    INVERSION DE MATRICE PAR RESOLUTION SUCCESSIVE DE NZ SYSTEMES
C    LINEAIRES
C    A  MATRICE (NZ*NZ) A INVERSER EN ENTREE, MATRICE INVERSEE EN SORTIE
C    ICRIT=1 SI MATRICE SINGULIERE, 0 SINON
C    B  TABLEAU DE REELS DE DIMENSION NZ*NZ
C    IS  TABLEAU D'ENTIERS DE DIMENSION NZ
C    EPS PRECISION
C  ====================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
      DIMENSION A(NZ,1),B(NZ,1)
      DIMENSION IS(1)
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
         B(J,I)=0.D0
  10     CONTINUE
      B(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 41 L=NR,NZ
          ABSKL=ABS(A(K,L))
          IF(ABSKL.GT.PIVOT) THEN
            I=K
            J=L
            PIVOT=ABSKL
          ENDIF
  41    CONTINUE
  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(NR,L)
      A(NR,L)=A(I,L)
      A(I,L)=D
   50 CONTINUE
      DO 60 L=1,NZ
      E=B(NR,L)
      B(NR,L)=B(I,L)
      B(I,L)=E
   60 CONTINUE
C
C      CHANGEMENT DE COLONNE : PLACE LE PIVOT EN R EME COLONNE
C
      DO 70 M=1,NZ
      CC  =A(M,NR)
      A(M,NR)=A(M,J)
      A(M,J)=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(I,NR).NE.0.D0)THEN
        G=A(I,NR)/A(NR,NR)
        DO 80 J=1,NZ
        A(I,J)=A(I,J)-G*A(NR,J)
        B(I,J)=B(I,J)-G*B(NR,J)
   80   CONTINUE
      ENDIF
   90 CONTINUE
  100 CONTINUE
C
C  2-  RESOLUTION DU SYSTEME TRIANGULARISE
C
      DO 130 J=1,NZ
         IF (ABS(B(NZ,J)).GT.(1.-XZPREC)*ABS(B(NZ,J))) THEN
           IF (ABS(A(NZ,NZ)).LE.((1.-XZPREC)*ABS(A(NZ,NZ)))) THEN
              ICRIT = 1
              RETURN
           ENDIF
           B(NZ,J)=B(NZ,J)/A(NZ,NZ)
         ENDIF
         DO 120 I= NZM1,1,-1
         F=0.D0
         I1=I+1
            DO 110 JJ=I1,NZ
            F=F-A(I,JJ)*B(JJ,J)
  110       CONTINUE
         B(I,J)=(B(I,J)+F)/A(I,I)
  120    CONTINUE
  130 CONTINUE
C
      DO 140 L=1,NZ
         IL=IS(L)
         DO 141 J=1,NZ
           A(IL,J)=B(L,J)
  141   CONTINUE
  140 CONTINUE
      RETURN
      END