Télécharger ivmat.eso

Retour à la liste

Numérotation des lignes :

  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.  
  42. -INC PPARAM
  43. -INC CCOPTIO
  44. -INC CCREEL
  45. *
  46. INTEGER NPN
  47. REAL*8 PN(NPN,NPN)
  48. INTEGER INVTMP(NPN)
  49. REAL*8 PNM1(NPN,NPN)
  50. REAL*8 DETPN
  51. *
  52. INTEGER IMPR,IRET
  53. *
  54. INTEGER ITMP
  55. REAL*8 RTMP,PIV,UNSPIV
  56. INTEGER INPN,JNPN,KNPN
  57. LOGICAL PFOUND,JFOUND
  58. *
  59. * Executable statements
  60. *
  61. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ivmat'
  62. *
  63. DETPN=1.D0
  64. DO 1 INPN=1,NPN
  65. INVTMP(INPN)=INPN
  66. DO 12 JNPN=1,NPN
  67. PNM1(JNPN,INPN)=PN(JNPN,INPN)
  68. 12 CONTINUE
  69. 1 CONTINUE
  70. *
  71. * Parcourons les colonnes de la matrice
  72. *
  73. DO 3 JNPN=1,NPN
  74. * Cherchons le plus grand pivot sur la colonne JNPN
  75. INPN=JNPN
  76. XMAX=ABS(PNM1(INPN,JNPN))
  77. DO KNPN=JNPN+1,NPN
  78. XVAL=ABS(PNM1(KNPN,JNPN))
  79. IF (XVAL.GT.XMAX) THEN
  80. INPN=KNPN
  81. XMAX=XVAL
  82. ENDIF
  83. ENDDO
  84. PFOUND=(ABS(PNM1(INPN,JNPN)).GT.SQRT(XPETIT))
  85. * On n'en a pas trouvé
  86. IF (.NOT.PFOUND) THEN
  87. DETPN=XZERO
  88. IF (IMPR.GT.0) THEN
  89. WRITE(IOIMP,*)
  90. $ 'Je ne peux pas inverser la matrice fournie.'
  91. ENDIF
  92. GOTO 9999
  93. ENDIF
  94. * On en a trouvé : on échange la ligne JNPN et la ligne INPN si besoin
  95. PIV=PNM1(INPN,JNPN)
  96. DETPN=DETPN*PIV
  97. IF (INPN.NE.JNPN) THEN
  98. ITMP=INVTMP(JNPN)
  99. INVTMP(JNPN)=INVTMP(INPN)
  100. INVTMP(INPN)=ITMP
  101. DO 34 KNPN=1,NPN
  102. RTMP=PNM1(INPN,KNPN)
  103. PNM1(INPN,KNPN)=PNM1(JNPN,KNPN)
  104. PNM1(JNPN,KNPN)=RTMP
  105. 34 CONTINUE
  106. DETPN=-DETPN
  107. ENDIF
  108. * On normalise la ligne du pivot
  109. UNSPIV=1.D0/PIV
  110. PNM1(JNPN,JNPN)=1.D0
  111. DO 36 KNPN=1,NPN
  112. PNM1(JNPN,KNPN)=PNM1(JNPN,KNPN)*UNSPIV
  113. 36 CONTINUE
  114. * Elimination
  115. DO 38 INPN=1,NPN
  116. IF (INPN.NE.JNPN) THEN
  117. RTMP=PNM1(INPN,JNPN)
  118. PNM1(INPN,JNPN)=XZERO
  119. DO 382 KNPN=1,NPN
  120. PNM1(INPN,KNPN)=PNM1(INPN,KNPN)-RTMP*PNM1(JNPN,KNPN)
  121. 382 CONTINUE
  122. ENDIF
  123. 38 CONTINUE
  124. 3 CONTINUE
  125. * On réordonne les colonnes de l'inverse
  126. DO 5 KNPN=1,NPN
  127. * Cherchons JNPN tel que INVTMP(JNPN)=KNPN
  128. JFOUND=.FALSE.
  129. JNPN=KNPN-1
  130. * REPEAT..UNTIL
  131. 52 CONTINUE
  132. IF (.NOT.JFOUND.AND.JNPN.LT.NPN) THEN
  133. JNPN=JNPN+1
  134. JFOUND=(INVTMP(JNPN).EQ.KNPN)
  135. GOTO 52
  136. ENDIF
  137. * On n'a pas trouvé JNPN : ca n'est pas normal
  138. IF (.NOT.JFOUND) THEN
  139. WRITE(IOIMP,*) 'Je ne peux pas réordonner l''inverse de la'
  140. WRITE(IOIMP,*) 'matrice fournie.'
  141. GOTO 9999
  142. ENDIF
  143. * Echangeons les colonnes JNPN et KNPN
  144. IF (KNPN.NE.JNPN) THEN
  145. INVTMP(JNPN)=INVTMP(KNPN)
  146. DO 54 INPN=1,NPN
  147. RTMP=PNM1(INPN,KNPN)
  148. PNM1(INPN,KNPN)=PNM1(INPN,JNPN)
  149. PNM1(INPN,JNPN)=RTMP
  150. 54 CONTINUE
  151. ENDIF
  152. 5 CONTINUE
  153. IF (IMPR.GT.3) THEN
  154. WRITE(IOIMP,*) 'On a créé [PNM1] (',NPN,'x',NPN,') :'
  155. DO 7 INPN=1,NPN
  156. WRITE(IOIMP,4004) (PNM1(INPN,JNPN),JNPN=1,NPN)
  157. 7 CONTINUE
  158. WRITE(IOIMP,*) 'Le déterminant de PN vaut :',DETPN
  159. ENDIF
  160. *
  161. * Normal termination
  162. *
  163. IRET=0
  164. RETURN
  165. *
  166. * Format handling
  167. *
  168. 4004 FORMAT (2X,6(1X,1PE13.5))
  169. *
  170. * Error handling
  171. *
  172. 9999 CONTINUE
  173. IRET=1
  174. WRITE(IOIMP,*) 'An error was detected in subroutine ivmat'
  175. RETURN
  176. *
  177. * End of subroutine IVMAT
  178. *
  179. END
  180.  
  181.  
  182.  
  183.  
  184.  
  185.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales