Télécharger meilu0.eso

Retour à la liste

Numérotation des lignes :

  1. C MEILU0 SOURCE GOUNAND 05/02/16 21:16:06 5029
  2. SUBROUTINE MEILU0(NTT,NNZ,A,JA,IA,
  3. $ ALU,JLU,JU,
  4. $ IWORK,
  5. $ IMPR,IRET)
  6. IMPLICIT REAL*8 (A-H,O-Z)
  7. IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM : MEILU0
  10. C DESCRIPTION :
  11. C Calcul du préconditionneur ILU(0) d'une matrice Morse.
  12. C ILU(0) : Incomplete LU factorization of level 0
  13. C appelée aussi Choleski ou Crout incomplet
  14. C La matrice Morse est au format CSR (Column Sparse Row)
  15. C Le préconditionneur est au format MSR (Modified Sparse Row,
  16. C stockage de l'inverse de la diagonale) (ALU, JLU).
  17. C JU est un tableau supplémentaire pour le repérage des diagonales.
  18. C (voir la doc. de Sparskit version 2)
  19. C
  20. C LANGAGE : FORTRAN 77 (sauf E/S)
  21. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  22. C mél : gounand@semt2.smts.cea.fr
  23. C REFERENCE (bibtex-like) :
  24. C @BOOK{templates,
  25. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  26. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  27. C H. Van der Vorst},
  28. C TITLE={Templates for the Solution of Linear Systems :
  29. C Building Blocks for Iterative Methods},
  30. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  31. C -> URL : http://www.netlib.org/templates/Templates.html
  32. C Sparskit : a basic tool kit for sparse matrix computations
  33. C Version 2 (Youcef Saad)
  34. C -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  35. C
  36. C***********************************************************************
  37. C APPELES : -
  38. C APPELE PAR : PRILU0
  39. C***********************************************************************
  40. C ENTREES : NTT, NNZ, IA, JA, A
  41. C ENTREES/SORTIES : IWORK (tableau de travail)
  42. C SORTIES : JU, JLU, ALU
  43. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  44. C***********************************************************************
  45. C VERSION : v1.1, 21/03/00
  46. C HISTORIQUE : v1.1, 22/03/00, gounand
  47. C Amélioration de la détection des pivots nulles.
  48. C HISTORIQUE : v1, 20/12/99, création
  49. C HISTORIQUE :
  50. C HISTORIQUE :
  51. C***********************************************************************
  52. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  53. C en cas de modification de ce sous-programme afin de faciliter
  54. C la maintenance !
  55. C***********************************************************************
  56. -INC CCOPTIO
  57. -INC CCREEL
  58. * .. Scalar Arguments ..
  59. INTEGER NTT,NNZ
  60. INTEGER IMPR,IRET
  61. * ..
  62. * .. Array Arguments ..
  63. REAL*8 A(NNZ),ALU(NNZ+1)
  64. INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
  65. INTEGER IWORK(NTT)
  66. *
  67. * .. Variables locales
  68. * .. Parameters
  69. REAL*8 ZERO ,ONE
  70. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  71. * ..
  72. C Nombre de pivots petits
  73. INTEGER NBPIVP
  74. C Nombre de pivots inférieurs à 0
  75. INTEGER NBPIVI
  76. INTEGER ITT,JNZ,KNZ,JCOL,JROW
  77. INTEGER JSTRT,JMID,JSTOP,JU0,JW
  78. INTEGER ISTRT,ISTOP
  79. REAL*8 TL,VALPIV,TNORM,PETIT
  80.  
  81. *
  82. * Executable statements
  83. *
  84. PETIT=XZPREC**(0.75D0)
  85. NBPIVP=0
  86. NBPIVI=0
  87. JU0=NTT+2
  88. JLU(1)=JU0
  89. *
  90. * On boucle sur les ddl
  91. *
  92. DO 1 ITT=1,NTT
  93. JSTRT=JU0
  94. *
  95. * Création de la ligne ITT pour L et U
  96. *
  97. ISTRT=IA(ITT)
  98. ISTOP=IA(ITT+1)-1
  99. TNORM=ZERO
  100. DO 12 JNZ=ISTRT,ISTOP
  101. *
  102. * Copie de la ligne ITT du format CSR au format MSR
  103. * et calcul de la norme de la ligne
  104. *
  105. JCOL=JA(JNZ)
  106. TNORM=TNORM+ABS(A(JNZ))
  107. IF (JCOL.EQ.ITT) THEN
  108. ALU(ITT) =A(JNZ)
  109. IWORK(JCOL)=ITT
  110. JU(ITT) =JU0
  111. ELSE
  112. ALU(JU0) =A(JNZ)
  113. JLU(JU0) =JA(JNZ)
  114. IWORK(JCOL)=JU0
  115. JU0 = JU0+1
  116. ENDIF
  117. 12 CONTINUE
  118. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  119. IF (TNORM.LE.XPETIT) THEN
  120. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  121. GOTO 9999
  122. ENDIF
  123. JLU(ITT+1)=JU0
  124. JSTOP=JU0-1
  125. JMID=JU(ITT)-1
  126. DO 14 JNZ=JSTRT,JMID
  127. JROW=JLU(JNZ)
  128. TL=ALU(JNZ)*ALU(JROW)
  129. ALU(JNZ)=TL
  130. *
  131. * On effectue la combinaison linéaire
  132. *
  133. DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
  134. JW=IWORK(JLU(KNZ))
  135. IF (JW.NE.0) THEN
  136. ALU(JW)=ALU(JW)-TL*ALU(KNZ)
  137. ENDIF
  138. 142 CONTINUE
  139. 14 CONTINUE
  140. *
  141. * On s'occupe de l'élément diagonal
  142. *
  143. VALPIV=ALU(ITT)
  144. IF (VALPIV.EQ.ZERO) THEN
  145. WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  146. $ 'le préconditionnement ILU(0) est impossible.'
  147. IRET=-2
  148. GOTO 9999
  149. ELSE
  150. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  151. NBPIVP=NBPIVP+1
  152. ENDIF
  153. IF (VALPIV.LT.ZERO) THEN
  154. NBPIVI=NBPIVI+1
  155. ENDIF
  156. ENDIF
  157. ALU(ITT)=ONE/ALU(ITT)
  158. *
  159. * On remet le tableau de travail à zéro
  160. *
  161. IWORK(ITT)=0
  162. DO 16 JNZ=JSTRT,JSTOP
  163. IWORK(JLU(JNZ))=0
  164. 16 CONTINUE
  165. 1 CONTINUE
  166. *
  167. * Warning(s)
  168. *
  169. IF (IMPR.GT.0) THEN
  170. IF (NBPIVP.GT.0) THEN
  171. WRITE(IOIMP,*) 'WARNING !',
  172. $ ' meilu0 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  173. ENDIF
  174. *! IF (NBPIVI.GT.0) THEN
  175. *! WRITE(IOIMP,*) 'WARNING !',
  176. *! $ ' meilu0 :',NBPIVI,' pivot(s) < 0'
  177. *! ENDIF
  178. ENDIF
  179. *
  180. * Normal termination
  181. *
  182. IRET=0
  183. RETURN
  184. *
  185. * Error handling
  186. *
  187. 9999 CONTINUE
  188. IRET=1
  189. WRITE(IOIMP,*) 'An error was detected in meilu0.eso'
  190. RETURN
  191. *
  192. * End of MEILU0
  193. *
  194. END
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  

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