Télécharger meilup.eso

Retour à la liste

Numérotation des lignes :

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

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