Télécharger meilu0.eso

Retour à la liste

Numérotation des lignes :

meilu0
  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.  
  57. -INC PPARAM
  58. -INC CCOPTIO
  59. -INC CCREEL
  60. * .. Scalar Arguments ..
  61. INTEGER NTT,NNZ
  62. INTEGER IMPR,IRET
  63. * ..
  64. * .. Array Arguments ..
  65. REAL*8 A(NNZ),ALU(NNZ+1)
  66. INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
  67. INTEGER IWORK(NTT)
  68. *
  69. * .. Variables locales
  70. * .. Parameters
  71. REAL*8 ZERO ,ONE
  72. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  73. * ..
  74. C Nombre de pivots petits
  75. INTEGER NBPIVP
  76. C Nombre de pivots inférieurs à 0
  77. INTEGER NBPIVI
  78. INTEGER ITT,JNZ,KNZ,JCOL,JROW
  79. INTEGER JSTRT,JMID,JSTOP,JU0,JW
  80. INTEGER ISTRT,ISTOP
  81. REAL*8 TL,VALPIV,TNORM,PETIT
  82.  
  83. *
  84. * Executable statements
  85. *
  86. PETIT=XZPREC**(0.75D0)
  87. NBPIVP=0
  88. NBPIVI=0
  89. JU0=NTT+2
  90. JLU(1)=JU0
  91. *
  92. * On boucle sur les ddl
  93. *
  94. DO 1 ITT=1,NTT
  95. JSTRT=JU0
  96. *
  97. * Création de la ligne ITT pour L et U
  98. *
  99. ISTRT=IA(ITT)
  100. ISTOP=IA(ITT+1)-1
  101. TNORM=ZERO
  102. DO 12 JNZ=ISTRT,ISTOP
  103. *
  104. * Copie de la ligne ITT du format CSR au format MSR
  105. * et calcul de la norme de la ligne
  106. *
  107. JCOL=JA(JNZ)
  108. TNORM=TNORM+ABS(A(JNZ))
  109. IF (JCOL.EQ.ITT) THEN
  110. ALU(ITT) =A(JNZ)
  111. IWORK(JCOL)=ITT
  112. JU(ITT) =JU0
  113. ELSE
  114. ALU(JU0) =A(JNZ)
  115. JLU(JU0) =JA(JNZ)
  116. IWORK(JCOL)=JU0
  117. JU0 = JU0+1
  118. ENDIF
  119. 12 CONTINUE
  120. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  121. IF (TNORM.LE.XPETIT) THEN
  122. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  123. GOTO 9999
  124. ENDIF
  125. JLU(ITT+1)=JU0
  126. JSTOP=JU0-1
  127. JMID=JU(ITT)-1
  128. DO 14 JNZ=JSTRT,JMID
  129. JROW=JLU(JNZ)
  130. TL=ALU(JNZ)*ALU(JROW)
  131. ALU(JNZ)=TL
  132. *
  133. * On effectue la combinaison linéaire
  134. *
  135. DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
  136. JW=IWORK(JLU(KNZ))
  137. IF (JW.NE.0) THEN
  138. ALU(JW)=ALU(JW)-TL*ALU(KNZ)
  139. ENDIF
  140. 142 CONTINUE
  141. 14 CONTINUE
  142. *
  143. * On s'occupe de l'élément diagonal
  144. *
  145. VALPIV=ALU(ITT)
  146. IF (VALPIV.EQ.ZERO) THEN
  147. WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  148. $ 'le préconditionnement ILU(0) est impossible.'
  149. IRET=-2
  150. GOTO 9999
  151. ELSE
  152. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  153. NBPIVP=NBPIVP+1
  154. ENDIF
  155. IF (VALPIV.LT.ZERO) THEN
  156. NBPIVI=NBPIVI+1
  157. ENDIF
  158. ENDIF
  159. ALU(ITT)=ONE/ALU(ITT)
  160. *
  161. * On remet le tableau de travail à zéro
  162. *
  163. IWORK(ITT)=0
  164. DO 16 JNZ=JSTRT,JSTOP
  165. IWORK(JLU(JNZ))=0
  166. 16 CONTINUE
  167. 1 CONTINUE
  168. *
  169. * Warning(s)
  170. *
  171. IF (IMPR.GT.0) THEN
  172. IF (NBPIVP.GT.0) THEN
  173. WRITE(IOIMP,*) 'WARNING !',
  174. $ ' meilu0 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  175. ENDIF
  176. *! IF (NBPIVI.GT.0) THEN
  177. *! WRITE(IOIMP,*) 'WARNING !',
  178. *! $ ' meilu0 :',NBPIVI,' pivot(s) < 0'
  179. *! ENDIF
  180. ENDIF
  181. *
  182. * Normal termination
  183. *
  184. IRET=0
  185. RETURN
  186. *
  187. * Error handling
  188. *
  189. 9999 CONTINUE
  190. IRET=1
  191. WRITE(IOIMP,*) 'An error was detected in meilu0.eso'
  192. RETURN
  193. *
  194. * End of MEILU0
  195. *
  196. END
  197.  
  198.  
  199.  
  200.  
  201.  
  202.  
  203.  

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