Télécharger mdilus.eso

Retour à la liste

Numérotation des lignes :

  1. C MDILUS SOURCE CHAT 05/01/13 01:40:17 5004
  2. SUBROUTINE MDILUS(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 : MDILUS
  10. C DESCRIPTION :
  11. C Calcul du préconditionneur D-ILU d'une matrice Morse
  12. C SYMETRIQUE.
  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 + chouhia 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. REAL*8 DINV(N)
  67. *
  68. * .. Variables locales
  69. * .. Parameters
  70. REAL*8 ZERO ,PETIT ,ONE
  71. PARAMETER (ZERO=0.0D0,PETIT=1.0D-14,ONE=1.0D0)
  72. * ..
  73. C Nombre de pivots petits
  74. INTEGER NBPIVP
  75. C Nombre de pivots inférieurs à 0
  76. INTEGER NBPIVI
  77. REAL*8 VALPIV,VALDIA
  78. C***
  79. IRET=0
  80. NBPIVP=0
  81. NBPIVI=0
  82. C On boucle sur les lignes
  83. DO 1 IN=1,N
  84. IDEB=ROWPTR(IN)
  85. IFIN=ROWPTR(IN+1)-1
  86. ICOL=IDEB
  87. IF (IDEB.LE.IFIN) THEN
  88. C On cherche le terme Aii
  89. 11 CONTINUE
  90. IF (COLIND(ICOL).LT.IN.AND.ICOL.LT.IFIN) THEN
  91. ICOL=ICOL+1
  92. GOTO 11
  93. ENDIF
  94. C On ne l'a pas trouvé
  95. IF (COLIND(ICOL).NE.IN) THEN
  96. WRITE(IOIMP,*) 'diag.',IN,'inexistante'
  97. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  98. WRITE(IOIMP,*) 'est impossible.'
  99. IRET=-1
  100. GOTO 9999
  101. ENDIF
  102. C On l'a trouvé
  103. C Dii <- Dii+Aii
  104. C Dii <- 1/Dii
  105. VALDIA=VAL(ICOL)
  106. VALPIV=DINV(IN)+VALDIA
  107. IF (VALPIV.EQ.ZERO) THEN
  108. WRITE(IOIMP,*) 'Pivot',IN,'nul'
  109. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  110. WRITE(IOIMP,*) 'est impossible.'
  111. IRET=-2
  112. GOTO 9999
  113. ELSE
  114. IF (ABS(VALPIV).LT.PETIT) THEN
  115. NBPIVP=NBPIVP+1
  116. ENDIF
  117. IF (VALPIV.LT.ZERO) THEN
  118. NBPIVI=NBPIVI+1
  119. ENDIF
  120. ENDIF
  121. VALPIV=ONE/VALPIV
  122. DINV(IN)=VALPIV
  123. C pour j>i : Djj <- Djj-AjiDiiAij
  124. JCOL=ICOL+1
  125. JFIN=IFIN
  126. IF (JCOL.LE.JFIN) THEN
  127. 12 CONTINUE
  128. JN=COLIND(JCOL)
  129. AIJ=VAL(JCOL)
  130. DINV(JN)=DINV(JN)-(VALPIV*AIJ*AIJ)
  131. JCOL=JCOL+1
  132. IF (JCOL.LE.JFIN) GOTO 12
  133. ENDIF
  134. ELSE
  135. WRITE(IOIMP,*) 'Ligne',IN,'vide'
  136. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  137. WRITE(IOIMP,*) 'est impossible.'
  138. IRET=-3
  139. GOTO 9999
  140. ENDIF
  141. 1 CONTINUE
  142. *
  143. * Warning(s)
  144. *
  145. IF (NBPIVP.GT.0) THEN
  146. WRITE(IOIMP,*) 'WARNING !'
  147. WRITE(IOIMP,*) 'mdilus :',NBPIVP,' |pivot(s)|<',PETIT
  148. ENDIF
  149. IF (NBPIVI.GT.0) THEN
  150. WRITE(IOIMP,*) 'WARNING !'
  151. WRITE(IOIMP,*) 'mdilus :',NBPIVI,' pivot(s) < 0'
  152. ENDIF
  153. *
  154. * Normal termination
  155. *
  156. RETURN
  157. *
  158. * Format handling
  159. *
  160. *
  161. * Error handling
  162. *
  163. 9999 CONTINUE
  164. WRITE(IOIMP,*) 'An error was detected in mdilus.eso'
  165. RETURN
  166. *
  167. * End of MDILUS
  168. *
  169. END
  170.  
  171.  
  172.  
  173.  
  174.  
  175.  
  176.  
  177.  
  178.  

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