1 : C IVMAT SOURCE GOUNAND 06/08/01 21:15:03 5512 2 : SUBROUTINE IVMAT(NPN,PN, 3 : $ INVTMP, 4 : $ PNM1,DETPN, 5 : $ IMPR,IRET) 6 : IMPLICIT INTEGER(I-N) 7 : IMPLICIT REAL*8 (A-H,O-Z) 8 : C*********************************************************************** 9 : C NOM : IVMAT 10 : C PROJET : Noyau linéaire NLIN 11 : C DESCRIPTION : Calcul de l'inverse d'une matrice pleine [PN] non 12 : C symétrique et de son déterminant. 13 : C NPN est son ordre. 14 : C INVTMP est un vecteur de travail temporaire 15 : C PNM1 est la matrice inversée. 16 : C DETPN est le déterminant de PN. 17 : C 18 : C LANGAGE : Fortran 77 (sauf pour les E/S) 19 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 20 : C mél : gounand@semt2.smts.cea.fr 21 : C*********************************************************************** 22 : C APPELES : - 23 : C APPELE PAR : KFNREF 24 : C*********************************************************************** 25 : C ENTREES : NPN, PN 26 : C ENTREES/SORTIES : INVTMP 27 : C SORTIES : PNM1, DETPN 28 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 29 : C*********************************************************************** 30 : C VERSION : v2, 01/08/06, recherche du plus grand pivot au lieu du 31 : C premier non nul... 32 : C VERSION : v1, 29/07/99, version initiale 33 : C HISTORIQUE : v1, 29/07/99, création 34 : C HISTORIQUE : 35 : C HISTORIQUE : 36 : C*********************************************************************** 37 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 38 : C en cas de modification de ce sous-programme afin de faciliter 39 : C la maintenance ! 40 : C*********************************************************************** 41 : -INC CCOPTIO 42 : -INC CCREEL 43 : * 44 : INTEGER NPN 45 : REAL*8 PN(NPN,NPN) 46 : INTEGER INVTMP(NPN) 47 : REAL*8 PNM1(NPN,NPN) 48 : REAL*8 DETPN 49 : * 50 : INTEGER IMPR,IRET 51 : * 52 : INTEGER ITMP 53 : REAL*8 RTMP,PIV,UNSPIV 54 : INTEGER INPN,JNPN,KNPN 55 : LOGICAL PFOUND,JFOUND 56 : * 57 : * Executable statements 58 : * 59 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ivmat' 60 : * 61 : DETPN=1.D0 62 : DO 1 INPN=1,NPN 63 : INVTMP(INPN)=INPN 64 : DO 12 JNPN=1,NPN 65 : PNM1(JNPN,INPN)=PN(JNPN,INPN) 66 : 12 CONTINUE 67 : 1 CONTINUE 68 : * 69 : * Parcourons les colonnes de la matrice 70 : * 71 : DO 3 JNPN=1,NPN 72 : * Cherchons le plus grand pivot sur la colonne JNPN 73 : INPN=JNPN 74 : XMAX=ABS(PNM1(INPN,JNPN)) 75 : DO KNPN=JNPN+1,NPN 76 : XVAL=ABS(PNM1(KNPN,JNPN)) 77 : IF (XVAL.GT.XMAX) THEN 78 : INPN=KNPN 79 : XMAX=XVAL 80 : ENDIF 81 : ENDDO 82 : PFOUND=(ABS(PNM1(INPN,JNPN)).GT.SQRT(XPETIT)) 83 : * On n'en a pas trouvé 84 : IF (.NOT.PFOUND) THEN 85 : DETPN=XZERO 86 : IF (IMPR.GT.0) THEN 87 : WRITE(IOIMP,*) 88 : $ 'Je ne peux pas inverser la matrice fournie.' 89 : ENDIF 90 : GOTO 9999 91 : ENDIF 92 : * On en a trouvé : on échange la ligne JNPN et la ligne INPN si besoin 93 : PIV=PNM1(INPN,JNPN) 94 : DETPN=DETPN*PIV 95 : IF (INPN.NE.JNPN) THEN 96 : ITMP=INVTMP(JNPN) 97 : INVTMP(JNPN)=INVTMP(INPN) 98 : INVTMP(INPN)=ITMP 99 : DO 34 KNPN=1,NPN 100 : RTMP=PNM1(INPN,KNPN) 101 : PNM1(INPN,KNPN)=PNM1(JNPN,KNPN) 102 : PNM1(JNPN,KNPN)=RTMP 103 : 34 CONTINUE 104 : DETPN=-DETPN 105 : ENDIF 106 : * On normalise la ligne du pivot 107 : UNSPIV=1.D0/PIV 108 : PNM1(JNPN,JNPN)=1.D0 109 : DO 36 KNPN=1,NPN 110 : PNM1(JNPN,KNPN)=PNM1(JNPN,KNPN)*UNSPIV 111 : 36 CONTINUE 112 : * Elimination 113 : DO 38 INPN=1,NPN 114 : IF (INPN.NE.JNPN) THEN 115 : RTMP=PNM1(INPN,JNPN) 116 : PNM1(INPN,JNPN)=XZERO 117 : DO 382 KNPN=1,NPN 118 : PNM1(INPN,KNPN)=PNM1(INPN,KNPN)-RTMP*PNM1(JNPN,KNPN) 119 : 382 CONTINUE 120 : ENDIF 121 : 38 CONTINUE 122 : 3 CONTINUE 123 : * On réordonne les colonnes de l'inverse 124 : DO 5 KNPN=1,NPN 125 : * Cherchons JNPN tel que INVTMP(JNPN)=KNPN 126 : JFOUND=.FALSE. 127 : JNPN=KNPN-1 128 : * REPEAT..UNTIL 129 : 52 CONTINUE 130 : IF (.NOT.JFOUND.AND.JNPN.LT.NPN) THEN 131 : JNPN=JNPN+1 132 : JFOUND=(INVTMP(JNPN).EQ.KNPN) 133 : GOTO 52 134 : ENDIF 135 : * On n'a pas trouvé JNPN : ca n'est pas normal 136 : IF (.NOT.JFOUND) THEN 137 : WRITE(IOIMP,*) 'Je ne peux pas réordonner l''inverse de la' 138 : WRITE(IOIMP,*) 'matrice fournie.' 139 : GOTO 9999 140 : ENDIF 141 : * Echangeons les colonnes JNPN et KNPN 142 : IF (KNPN.NE.JNPN) THEN 143 : INVTMP(JNPN)=INVTMP(KNPN) 144 : DO 54 INPN=1,NPN 145 : RTMP=PNM1(INPN,KNPN) 146 : PNM1(INPN,KNPN)=PNM1(INPN,JNPN) 147 : PNM1(INPN,JNPN)=RTMP 148 : 54 CONTINUE 149 : ENDIF 150 : 5 CONTINUE 151 : IF (IMPR.GT.3) THEN 152 : WRITE(IOIMP,*) 'On a créé [PNM1] (',NPN,'x',NPN,') :' 153 : DO 7 INPN=1,NPN 154 : WRITE(IOIMP,4004) (PNM1(INPN,JNPN),JNPN=1,NPN) 155 : 7 CONTINUE 156 : WRITE(IOIMP,*) 'Le déterminant de PN vaut :',DETPN 157 : ENDIF 158 : * 159 : * Normal termination 160 : * 161 : IRET=0 162 : RETURN 163 : * 164 : * Format handling 165 : * 166 : 4004 FORMAT (2X,6(1X,1PE13.5)) 167 : * 168 : * Error handling 169 : * 170 : 9999 CONTINUE 171 : IRET=1 172 : WRITE(IOIMP,*) 'An error was detected in subroutine ivmat' 173 : RETURN 174 : * 175 : * End of subroutine IVMAT 176 : * 177 : END 178 : 179 : 180 : 181 : 182 :
© Cast3M 2003 - All rights reserved.
Disclaimer