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.  
  56. -INC PPARAM
  57. -INC CCOPTIO
  58. * .. Scalar Arguments ..
  59. INTEGER N,NNZ
  60. INTEGER IRET
  61. * ..
  62. * .. Array Arguments ..
  63. * .. Matrices stockées en Morse
  64. INTEGER ROWPTR(N+1)
  65. INTEGER COLIND(NNZ)
  66. REAL*8 VAL(NNZ)
  67. * .. Vecteurs
  68. C INTEGER TMPCOL(N)
  69. C REAL*8 TMPVAL(N)
  70. REAL*8 DINV(N)
  71. *
  72. * .. Variables locales
  73. * .. Parameters
  74. REAL*8 ZERO ,PETIT ,ONE
  75. PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
  76. * ..
  77. C Nombre de pivots petits
  78. INTEGER NBPIVP
  79. C Nombre de pivots inférieurs à 0
  80. INTEGER NBPIVI
  81. REAL*8 VALPIV,VALDIA
  82. C***
  83. IRET=0
  84. NBPIVP=0
  85. NBPIVI=0
  86. C On boucle sur les lignes
  87. DO 1 IN=1,N
  88. IDEB=ROWPTR(IN)
  89. IFIN=ROWPTR(IN+1)-1
  90. ICOL=IDEB
  91. IF (IDEB.LE.IFIN) THEN
  92. C On cherche le terme Aii
  93. 11 CONTINUE
  94. IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
  95. ICOL=ICOL+1
  96. GOTO 11
  97. ENDIF
  98. C On ne l'a pas trouvé
  99. IF (COLIND(ICOL).NE.IN) THEN
  100. WRITE(IOIMP,*) 'diag.',IN,'inexistante'
  101. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  102. WRITE(IOIMP,*) 'est impossible.'
  103. IRET=-1
  104. GOTO 9999
  105. ELSE
  106. C On l'a trouvé
  107. C Dii <- Dii+Aii
  108. C Dii <- 1/Dii
  109. VALDIA=VAL(ICOL)
  110. VALPIV=DINV(IN)+VALDIA
  111. IF (VALPIV.EQ.ZERO) THEN
  112. WRITE(IOIMP,*) 'Pivot',IN,'nul'
  113. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  114. WRITE(IOIMP,*) 'est impossible.'
  115. IRET=-2
  116. GOTO 9999
  117. ELSE
  118. IF (ABS(VALPIV).LT.PETIT) THEN
  119. NBPIVP=NBPIVP+1
  120. ENDIF
  121. IF (VALPIV.LT.ZERO) THEN
  122. NBPIVI=NBPIVI+1
  123. ENDIF
  124. ENDIF
  125. VALPIV=ONE/VALPIV
  126. DINV(IN)=VALPIV
  127. JCOL=ICOL+1
  128. JFIN=IFIN
  129. IF (JCOL.LE.JFIN) THEN
  130. C On parcourt les valeurs Aij (j>i).
  131. C Pour chaque valeur, on regarde ligne j s'il y a une
  132. C valeur Aji correspondante.
  133. C Si oui, Djj <- Djj - AjiDiiAij
  134. C Sinon, rien !
  135. 12 CONTINUE
  136. JN=COLIND(JCOL)
  137. C Recherche de Aji
  138. KDEB=ROWPTR(JN)
  139. KFIN=ROWPTR(JN+1)-1
  140. KCOL=KDEB
  141. 121 CONTINUE
  142. IF (COLIND(KCOL).LT.IN.AND.KCOL.LT.KFIN) THEN
  143. KCOL=KCOL+1
  144. GOTO 121
  145. ENDIF
  146. C On a trouve Aji !
  147. IF (COLIND(KCOL).EQ.IN) THEN
  148. DINV(JN)=DINV(JN)-
  149. $ (VAL(KCOL)*VALPIV*VAL(JCOL))
  150. ENDIF
  151. C On boucle...
  152. JCOL=JCOL+1
  153. IF (JCOL.LE.JFIN) GOTO 12
  154. ENDIF
  155. ENDIF
  156. ELSE
  157. WRITE(IOIMP,*) 'Ligne',IN,'vide'
  158. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  159. WRITE(IOIMP,*) 'est impossible.'
  160. IRET=-3
  161. GOTO 9999
  162. ENDIF
  163. 1 CONTINUE
  164. *
  165. * Warning(s)
  166. *
  167. IF (NBPIVP.GT.0) THEN
  168. WRITE(IOIMP,*) 'WARNING !'
  169. WRITE(IOIMP,*) 'mdilun :',NBPIVP,' |pivot(s)|<',PETIT
  170. ENDIF
  171. *! IF (NBPIVI.GT.0) THEN
  172. *! WRITE(IOIMP,*) 'WARNING !'
  173. *! WRITE(IOIMP,*) 'mdilun :',NBPIVI,' pivot(s) < 0'
  174. *! ENDIF
  175. *
  176. * Normal termination
  177. *
  178. RETURN
  179. *
  180. * Error handling
  181. *
  182. 9999 CONTINUE
  183. WRITE(IOIMP,*) 'An error was detected in mdilun.eso'
  184. RETURN
  185. *
  186. * End of MDILUN
  187. *
  188. END
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.  
  196.  
  197.  
  198.  

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