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.  
  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. 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. ENDIF
  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. C pour j>i : Djj <- Djj-AjiDiiAij
  126. JCOL=ICOL+1
  127. JFIN=IFIN
  128. IF (JCOL.LE.JFIN) THEN
  129. 12 CONTINUE
  130. JN=COLIND(JCOL)
  131. AIJ=VAL(JCOL)
  132. DINV(JN)=DINV(JN)-(VALPIV*AIJ*AIJ)
  133. JCOL=JCOL+1
  134. IF (JCOL.LE.JFIN) GOTO 12
  135. ENDIF
  136. ELSE
  137. WRITE(IOIMP,*) 'Ligne',IN,'vide'
  138. WRITE(IOIMP,*) 'le préconditionnement D-ILU'
  139. WRITE(IOIMP,*) 'est impossible.'
  140. IRET=-3
  141. GOTO 9999
  142. ENDIF
  143. 1 CONTINUE
  144. *
  145. * Warning(s)
  146. *
  147. IF (NBPIVP.GT.0) THEN
  148. WRITE(IOIMP,*) 'WARNING !'
  149. WRITE(IOIMP,*) 'mdilus :',NBPIVP,' |pivot(s)|<',PETIT
  150. ENDIF
  151. IF (NBPIVI.GT.0) THEN
  152. WRITE(IOIMP,*) 'WARNING !'
  153. WRITE(IOIMP,*) 'mdilus :',NBPIVI,' pivot(s) < 0'
  154. ENDIF
  155. *
  156. * Normal termination
  157. *
  158. RETURN
  159. *
  160. * Format handling
  161. *
  162. *
  163. * Error handling
  164. *
  165. 9999 CONTINUE
  166. WRITE(IOIMP,*) 'An error was detected in mdilus.eso'
  167. RETURN
  168. *
  169. * End of MDILUS
  170. *
  171. END
  172.  
  173.  
  174.  
  175.  
  176.  
  177.  
  178.  
  179.  
  180.  

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