Sources Esope | Cast3M

invere
C INVERE    SOURCE    GOUNAND   25/10/23    21:15:03     12386          
        SUBROUTINE INVERE(VP,N,IVP,K,DET)
C******************************************************************************
C       INVERSION D'UNE MATRICE NON SYMETRIQUE AVEC RECHERCHE DE PIVOT NON NUL
C       SUR UNE COLONNE
C       ENTREES
C       VP      MATRICE A INVERSER
C       N       ORDRE DE LA MATRICE A INVERSER
C       IVP     DIMENSION DE LA MATRICE DANS LE PROGRAMME D'APPEL
C       K       VECTEUR DE TRAVAIL ENTIER DE LONGUEUR N
C       SORTIES
C       VP      MATRICE INVERSEE
C       DET     DETERMINANT
C     SG: 21/10/2025 : recherche du plus grand pivot au lieu du premier
C                      non nul pour la precision
C******************************************************************************
        IMPLICIT INTEGER(I-N)
        IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL
        DIMENSION VP(IVP,IVP),K(N)
        DATA UN/1.D0/
C******************************************************************************
        DET=UN
        DO 5 I=1,N
           K(I)=I
 5      CONTINUE
        XNORM=XZERO
        DO J=1,N
           DO I=1,N
              XNORM=MAX(XNORM,ABS(VP(I,J)))
           ENDDO
        ENDDO
        XNORM=MAX(XNORM,XPETIT/XZPREC)*10.D0
C****************DEBUT DE L'INVERSION
        DO 80 II=1,N
C****************RECHERCHE D'UN PIVOT NON NUL SUR LA COLONNE II
           IMAX=II
           XMAX=ABS(VP(II,II))
           DO 10 I=II+1,N
              XCOU=ABS(VP(I,II))
              IF(ABS(XCOU).GT.XMAX) THEN
                 IMAX=I
                 XMAX=XCOU
              ENDIF
 10        CONTINUE
           IF(XMAX.GT.XNORM*XZPREC) THEN
              I=IMAX
              PIV=VP(IMAX,II)
              GOTO 20
           ENDIF
           DET=XZERO
           RETURN
 
C****************ECHANGER LA LIGNE II ET LA LIGNE I
 20        CONTINUE
           DET=DET*PIV
           IF(I.EQ.II) GOTO 40
           I1=K(II)
           K(II)=K(I)
           K(I)=I1
           DO 30 J=1,N
              C=VP(I,J)
              VP(I,J)=VP(II,J)
              VP(II,J)=C
 30        CONTINUE
           DET=-DET
C****************NORMALISATION DE LA LIGNE DU PIVOT
 40        CONTINUE
           C=UN/PIV
           VP(II,II)=UN
           DO 50 J=1,N
              VP(II,J)=VP(II,J)*C
 50        CONTINUE
C****************ELIMINATION
           DO 70 I=1,N
              IF(I.EQ.II) GOTO 70
              C=VP(I,II)
              VP(I,II)=XZERO
              DO 60 J=1,N
                 VP(I,J)=VP(I,J)-C*VP(II,J)
 60           CONTINUE
 70        CONTINUE
 80     CONTINUE
C****************REORDONNER LES COLONNES DE L'INVERSE
        DO 120 J=1,N
C****************CHERCHER J1 TEL QUE K(J1)=J
           DO 90 J1=J,N
              JJ=K(J1)
              IF(JJ.EQ.J) GOTO 100
 90        CONTINUE
 100       CONTINUE
           IF(J.EQ.J1) GOTO 120
C****************ECHANGER LES COLONNES J ET J1
           K(J1)=K(J)
           DO 110 I=1,N
              C=VP(I,J)
              VP(I,J)=VP(I,J1)
              VP(I,J1)=C
 110       CONTINUE
 120    CONTINUE
        RETURN
        END