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. -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.  
  183.  

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