Sources Esope | Cast3M

ivmat
C IVMAT     SOURCE    GOUNAND   06/08/01    21:15:03     5512
      SUBROUTINE IVMAT(NPN,PN,
     $     INVTMP,
     $     PNM1,DETPN,
     $     IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : IVMAT
C PROJET      : Noyau linéaire NLIN
C DESCRIPTION : Calcul de l'inverse d'une matrice pleine [PN] non
C               symétrique et de son déterminant.
C               NPN est son ordre.
C               INVTMP est un vecteur de travail temporaire
C               PNM1 est la matrice inversée.
C               DETPN est le déterminant de PN.
C
C LANGAGE     : Fortran 77 (sauf pour les E/S)
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : -
C APPELE PAR       : KFNREF
C***********************************************************************
C ENTREES            : NPN, PN
C ENTREES/SORTIES    : INVTMP
C SORTIES            : PNM1, DETPN
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v2, 01/08/06, recherche du plus grand pivot au lieu du
C                            premier non nul...
C VERSION    : v1, 29/07/99, version initiale
C HISTORIQUE : v1, 29/07/99, création
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
*
      INTEGER NPN
      REAL*8 PN(NPN,NPN)
      INTEGER INVTMP(NPN)
      REAL*8 PNM1(NPN,NPN)
      REAL*8 DETPN
*
      INTEGER IMPR,IRET
*
      INTEGER ITMP
      REAL*8 RTMP,PIV,UNSPIV
      INTEGER INPN,JNPN,KNPN
      LOGICAL PFOUND,JFOUND
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ivmat'
*
      DETPN=1.D0
      DO 1 INPN=1,NPN
         INVTMP(INPN)=INPN
         DO 12 JNPN=1,NPN
            PNM1(JNPN,INPN)=PN(JNPN,INPN)
 12      CONTINUE
 1    CONTINUE
*
* Parcourons les colonnes de la matrice
*
      DO 3 JNPN=1,NPN
* Cherchons le plus grand pivot sur la colonne JNPN
         INPN=JNPN
         XMAX=ABS(PNM1(INPN,JNPN))
         DO KNPN=JNPN+1,NPN
            XVAL=ABS(PNM1(KNPN,JNPN))
            IF (XVAL.GT.XMAX) THEN
               INPN=KNPN
               XMAX=XVAL
            ENDIF
         ENDDO
         PFOUND=(ABS(PNM1(INPN,JNPN)).GT.SQRT(XPETIT))
* On n'en a pas trouvé
         IF (.NOT.PFOUND) THEN
            DETPN=XZERO
            IF (IMPR.GT.0) THEN
               WRITE(IOIMP,*)
     $              'Je ne peux pas inverser la matrice fournie.'
            ENDIF
            GOTO 9999
         ENDIF
* On en a trouvé : on échange la ligne JNPN et la ligne INPN si besoin
         PIV=PNM1(INPN,JNPN)
         DETPN=DETPN*PIV
         IF (INPN.NE.JNPN) THEN
            ITMP=INVTMP(JNPN)
            INVTMP(JNPN)=INVTMP(INPN)
            INVTMP(INPN)=ITMP
            DO 34 KNPN=1,NPN
               RTMP=PNM1(INPN,KNPN)
               PNM1(INPN,KNPN)=PNM1(JNPN,KNPN)
               PNM1(JNPN,KNPN)=RTMP
 34         CONTINUE
            DETPN=-DETPN
         ENDIF
* On normalise la ligne du pivot
         UNSPIV=1.D0/PIV
         PNM1(JNPN,JNPN)=1.D0
         DO 36 KNPN=1,NPN
            PNM1(JNPN,KNPN)=PNM1(JNPN,KNPN)*UNSPIV
 36      CONTINUE
* Elimination
         DO 38 INPN=1,NPN
            IF (INPN.NE.JNPN) THEN
               RTMP=PNM1(INPN,JNPN)
               PNM1(INPN,JNPN)=XZERO
               DO 382 KNPN=1,NPN
                  PNM1(INPN,KNPN)=PNM1(INPN,KNPN)-RTMP*PNM1(JNPN,KNPN)
 382           CONTINUE
            ENDIF
 38      CONTINUE
 3    CONTINUE
* On réordonne les colonnes de l'inverse
      DO 5 KNPN=1,NPN
* Cherchons JNPN tel que INVTMP(JNPN)=KNPN
         JFOUND=.FALSE.
         JNPN=KNPN-1
* REPEAT..UNTIL
 52      CONTINUE
         IF (.NOT.JFOUND.AND.JNPN.LT.NPN) THEN
            JNPN=JNPN+1
            JFOUND=(INVTMP(JNPN).EQ.KNPN)
            GOTO 52
         ENDIF
* On n'a pas trouvé JNPN : ca n'est pas normal
         IF (.NOT.JFOUND) THEN
            WRITE(IOIMP,*) 'Je ne peux pas réordonner l''inverse de la'
            WRITE(IOIMP,*) 'matrice fournie.'
            GOTO 9999
         ENDIF
* Echangeons les colonnes JNPN et KNPN
         IF (KNPN.NE.JNPN) THEN
            INVTMP(JNPN)=INVTMP(KNPN)
            DO 54 INPN=1,NPN
               RTMP=PNM1(INPN,KNPN)
               PNM1(INPN,KNPN)=PNM1(INPN,JNPN)
               PNM1(INPN,JNPN)=RTMP
 54         CONTINUE
         ENDIF
 5    CONTINUE
      IF (IMPR.GT.3) THEN
         WRITE(IOIMP,*) 'On a créé [PNM1] (',NPN,'x',NPN,') :'
         DO 7 INPN=1,NPN
            WRITE(IOIMP,4004) (PNM1(INPN,JNPN),JNPN=1,NPN)
 7       CONTINUE
         WRITE(IOIMP,*) 'Le déterminant de PN vaut :',DETPN
      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Format handling
*
 4004 FORMAT (2X,6(1X,1PE13.5))
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine ivmat'
      RETURN
*
* End of subroutine IVMAT
*
      END