Télécharger mdilun.eso

Retour à la liste

Numérotation des lignes :

  1. C MDILUN SOURCE GOUNAND 06/03/06 21:18:10 5319
  2. SUBROUTINE MDILUN(N,NNZ,
  3. $ ROWPTR,COLIND,VAL,
  4. $ DINV,
  5. $ IRET)
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. C***********************************************************************
  9. C NOM : MDILUN
  10. C DESCRIPTION :
  11. C Calcul du préconditionneur D-ILU d'une matrice Morse
  12. C NON symétrique.
  13. C D-ILU : Diagonal Incomplete LU factorization
  14. C
  15. C Ce sous-programme est appelé par prdilu.eso.
  16. C Il est écrit en Fortran car il est censé etre rapide
  17. C
  18. C LANGAGE : FORTRAN 77 + chouilla ESOPE (E/S)
  19. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  20. C mél : gounand@semt2.smts.cea.fr
  21. C REFERENCE (bibtex-like) :
  22. C @BOOK{templates,
  23. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  24. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  25. C H. Van der Vorst},
  26. C TITLE={Templates for the Solution of Linear Systems :
  27. C Building Blocks for Iterative Methods},
  28. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  29. C -> URL : http://www.netlib.org/templates/Templates.html
  30. C***********************************************************************
  31. C APPELES : -
  32. C***********************************************************************
  33. C ENTREES : N, NNZ,
  34. C ROWPTR, COLIND, VAL
  35. C ENTREES/SORTIES : -
  36. C SORTIES : DINVPIV, IRET
  37. C CODE RETOUR (IRET) : 0 si ok
  38. C <0 si problème (diagonale nulle...)
  39. C N : nombre de degrés de liberté
  40. C NNZ : nombre de valeurs non nulles de la matrice Morse
  41. C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
  42. C et valeurs de la matrice Morse
  43. C DINVPIV : vecteur contenant l'inverse de la factorisation
  44. C D-ILU de la matrice morse ci-dessus.
  45. C***********************************************************************
  46. C VERSION : v1, 01/04/98, version initiale
  47. C HISTORIQUE : v1, 01/04/98, création
  48. C HISTORIQUE :
  49. C HISTORIQUE :
  50. C***********************************************************************
  51. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  52. C en cas de modification de ce sous-programme afin de faciliter
  53. C la maintenance !
  54. C***********************************************************************
  55. -INC CCOPTIO
  56. * .. Scalar Arguments ..
  57. INTEGER N,NNZ
  58. INTEGER IRET
  59. * ..
  60. * .. Array Arguments ..
  61. * .. Matrices stockées en Morse
  62. INTEGER ROWPTR(N+1)
  63. INTEGER COLIND(NNZ)
  64. REAL*8 VAL(NNZ)
  65. * .. Vecteurs
  66. C INTEGER TMPCOL(N)
  67. C REAL*8 TMPVAL(N)
  68. REAL*8 DINV(N)
  69. *
  70. * .. Variables locales
  71. * .. Parameters
  72. REAL*8 ZERO ,PETIT ,ONE
  73. PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
  74. * ..
  75. C Nombre de pivots petits
  76. INTEGER NBPIVP
  77. C Nombre de pivots inférieurs à 0
  78. INTEGER NBPIVI
  79. REAL*8 VALPIV,VALDIA
  80. C***
  81. IRET=0
  82. NBPIVP=0
  83. NBPIVI=0
  84. C On boucle sur les lignes
  85. DO 1 IN=1,N
  86. IDEB=ROWPTR(IN)
  87. IFIN=ROWPTR(IN+1)-1
  88. ICOL=IDEB
  89. IF (IDEB.LE.IFIN) THEN
  90. C On cherche le terme Aii
  91. 11 CONTINUE
  92. IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
  93. ICOL=ICOL+1
  94. GOTO 11
  95. ENDIF
  96. C On ne l'a pas trouvé
  97. IF (COLIND(ICOL).NE.IN) THEN
  98. WRITE(IOIMP,*) 'diag.',IN,'inexistante'
  99. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  100. WRITE(IOIMP,*) 'est impossible.'
  101. IRET=-1
  102. GOTO 9999
  103. ELSE
  104. C On l'a trouvé
  105. C Dii <- Dii+Aii
  106. C Dii <- 1/Dii
  107. VALDIA=VAL(ICOL)
  108. VALPIV=DINV(IN)+VALDIA
  109. IF (VALPIV.EQ.ZERO) THEN
  110. WRITE(IOIMP,*) 'Pivot',IN,'nul'
  111. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  112. WRITE(IOIMP,*) 'est impossible.'
  113. IRET=-2
  114. GOTO 9999
  115. ELSE
  116. IF (ABS(VALPIV).LT.PETIT) THEN
  117. NBPIVP=NBPIVP+1
  118. ENDIF
  119. IF (VALPIV.LT.ZERO) THEN
  120. NBPIVI=NBPIVI+1
  121. ENDIF
  122. ENDIF
  123. VALPIV=ONE/VALPIV
  124. DINV(IN)=VALPIV
  125. JCOL=ICOL+1
  126. JFIN=IFIN
  127. IF (JCOL.LE.JFIN) THEN
  128. C On parcourt les valeurs Aij (j>i).
  129. C Pour chaque valeur, on regarde ligne j s'il y a une
  130. C valeur Aji correspondante.
  131. C Si oui, Djj <- Djj - AjiDiiAij
  132. C Sinon, rien !
  133. 12 CONTINUE
  134. JN=COLIND(JCOL)
  135. C Recherche de Aji
  136. KDEB=ROWPTR(JN)
  137. KFIN=ROWPTR(JN+1)-1
  138. KCOL=KDEB
  139. 121 CONTINUE
  140. IF (COLIND(KCOL).LT.IN.AND.KCOL.LT.KFIN) THEN
  141. KCOL=KCOL+1
  142. GOTO 121
  143. ENDIF
  144. C On a trouve Aji !
  145. IF (COLIND(KCOL).EQ.IN) THEN
  146. DINV(JN)=DINV(JN)-
  147. $ (VAL(KCOL)*VALPIV*VAL(JCOL))
  148. ENDIF
  149. C On boucle...
  150. JCOL=JCOL+1
  151. IF (JCOL.LE.JFIN) GOTO 12
  152. ENDIF
  153. ENDIF
  154. ELSE
  155. WRITE(IOIMP,*) 'Ligne',IN,'vide'
  156. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  157. WRITE(IOIMP,*) 'est impossible.'
  158. IRET=-3
  159. GOTO 9999
  160. ENDIF
  161. 1 CONTINUE
  162. *
  163. * Warning(s)
  164. *
  165. IF (NBPIVP.GT.0) THEN
  166. WRITE(IOIMP,*) 'WARNING !'
  167. WRITE(IOIMP,*) 'mdilun :',NBPIVP,' |pivot(s)|<',PETIT
  168. ENDIF
  169. *! IF (NBPIVI.GT.0) THEN
  170. *! WRITE(IOIMP,*) 'WARNING !'
  171. *! WRITE(IOIMP,*) 'mdilun :',NBPIVI,' pivot(s) < 0'
  172. *! ENDIF
  173. *
  174. * Normal termination
  175. *
  176. RETURN
  177. *
  178. * Error handling
  179. *
  180. 9999 CONTINUE
  181. WRITE(IOIMP,*) 'An error was detected in mdilun.eso'
  182. RETURN
  183. *
  184. * End of MDILUN
  185. *
  186. END
  187.  
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  

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