Télécharger memilu.eso

Retour à la liste

Numérotation des lignes :

  1. C MEMILU SOURCE GOUNAND 06/03/06 21:18:28 5319
  2. SUBROUTINE MEMILU(NTT,NNZ,A,JA,IA,
  3. $ ALU,JLU,JU,
  4. $ IWORK,
  5. $ RXMILU,
  6. $ IMPR,IRET)
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. C***********************************************************************
  10. C NOM : MEMILU
  11. C DESCRIPTION :
  12. C Calcul du préconditionneur MILU(0) d'une matrice Morse
  13. C MILU(0) : Modified Incomplete LU factorization of level 0
  14. C appelée aussi Choleski ou Crout incomplet modifié !
  15. C Ici, c'est meme un MILU(0) relaxé (avec RXMILU)...
  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, RXMILU
  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 RXMILU : type REAL*8.
  47. C facteur de relaxation (normalement compris
  48. C entre 0.D0 et 1.D0)
  49. C 0.D0 : on se ramène à Choleski incomplet
  50. C 1.D0 : on se ramène à MILU(0) non relaxé.
  51. C***********************************************************************
  52. C VERSION : v1.1, 21/03/00
  53. C HISTORIQUE : v1.1, 22/03/00, gounand
  54. C Amélioration de la détection des pivots nulles.
  55. C HISTORIQUE : v1, 20/12/99, création
  56. C HISTORIQUE :
  57. C HISTORIQUE :
  58. C***********************************************************************
  59. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  60. C en cas de modification de ce sous-programme afin de faciliter
  61. C la maintenance !
  62. C***********************************************************************
  63. -INC CCOPTIO
  64. -INC CCREEL
  65. * .. Scalar Arguments ..
  66. INTEGER NTT,NNZ
  67. REAL*8 RXMILU
  68. INTEGER IMPR,IRET
  69. * ..
  70. * .. Array Arguments ..
  71. REAL*8 A(NNZ),ALU(NNZ+1)
  72. INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
  73. INTEGER IWORK(NNZ)
  74. *
  75. * .. Variables locales
  76. * .. Parameters
  77. REAL*8 ZERO ,ONE
  78. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  79. * ..
  80. C Nombre de pivots petits
  81. INTEGER NBPIVP
  82. C Nombre de pivots inférieurs à 0
  83. INTEGER NBPIVI
  84. INTEGER ITT,JNZ,KNZ,JCOL,JROW
  85. INTEGER JSTRT,JMID,JSTOP,JU0,JW
  86. INTEGER ISTRT,ISTOP
  87. REAL*8 TL,VALPIV,TNORM,PETIT
  88. REAL*8 FILVAL
  89. *
  90. * Executable statements
  91. *
  92. PETIT=XZPREC**(0.75D0)
  93. NBPIVP=0
  94. NBPIVI=0
  95. JU0=NTT+2
  96. JLU(1)=JU0
  97. *
  98. * On boucle sur les ddl
  99. *
  100. DO 1 ITT=1,NTT
  101. JSTRT=JU0
  102. *
  103. * Création de la ligne ITT pour L et U
  104. *
  105. ISTRT=IA(ITT)
  106. ISTOP=IA(ITT+1)-1
  107. TNORM=ZERO
  108. DO 12 JNZ=ISTRT,ISTOP
  109. *
  110. * Copie de la ligne ITT du format CSR au format MSR
  111. *
  112. JCOL=JA(JNZ)
  113. TNORM=TNORM+ABS(A(JNZ))
  114. IF (JCOL.EQ.ITT) THEN
  115. ALU(ITT) =A(JNZ)
  116. IWORK(JCOL)=ITT
  117. JU(ITT) =JU0
  118. ELSE
  119. ALU(JU0) =A(JNZ)
  120. JLU(JU0) =JA(JNZ)
  121. IWORK(JCOL)=JU0
  122. JU0 = JU0+1
  123. ENDIF
  124. 12 CONTINUE
  125. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  126. IF (TNORM.LE.XPETIT) THEN
  127. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  128. GOTO 9999
  129. ENDIF
  130. JLU(ITT+1)=JU0
  131. JSTOP=JU0-1
  132. JMID=JU(ITT)-1
  133. FILVAL=ZERO
  134. DO 14 JNZ=JSTRT,JMID
  135. JROW=JLU(JNZ)
  136. TL=ALU(JNZ)*ALU(JROW)
  137. ALU(JNZ)=TL
  138. *
  139. * On effectue la combinaison linéaire
  140. *
  141. DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
  142. JW=IWORK(JLU(KNZ))
  143. IF (JW.NE.0) THEN
  144. ALU(JW)=ALU(JW)-TL*ALU(KNZ)
  145. ELSE
  146. FILVAL=FILVAL+TL*ALU(KNZ)
  147. ENDIF
  148. 142 CONTINUE
  149. 14 CONTINUE
  150. *
  151. * On s'occupe de l'élément diagonal
  152. *
  153. ALU(ITT)=ALU(ITT)-(FILVAL*RXMILU)
  154. VALPIV=ALU(ITT)
  155. IF (VALPIV.EQ.ZERO) THEN
  156. WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  157. $ 'le préconditionnement MILU(0) est impossible.'
  158. IRET=-2
  159. GOTO 9999
  160. ELSE
  161. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  162. NBPIVP=NBPIVP+1
  163. ENDIF
  164. IF (VALPIV.LT.ZERO) THEN
  165. NBPIVI=NBPIVI+1
  166. ENDIF
  167. ENDIF
  168. ALU(ITT)=ONE/ALU(ITT)
  169. *
  170. * On remet le tableau de travail à zéro
  171. *
  172. IWORK(ITT)=0
  173. DO 16 JNZ=JSTRT,JSTOP
  174. IWORK(JLU(JNZ))=0
  175. 16 CONTINUE
  176. 1 CONTINUE
  177. *
  178. * Warning(s)
  179. *
  180. IF (IMPR.GT.0) THEN
  181. IF (NBPIVP.GT.0) THEN
  182. WRITE(IOIMP,*) 'WARNING !',
  183. $ ' memilu :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  184. ENDIF
  185. * IF (NBPIVI.GT.0) THEN
  186. * WRITE(IOIMP,*) 'WARNING !',
  187. * $ ' memilu :',NBPIVI,' pivot(s) < 0'
  188. * ENDIF
  189. ENDIF
  190. *
  191. * Normal termination
  192. *
  193. IRET=0
  194. RETURN
  195. *
  196. * Error handling
  197. *
  198. 9999 CONTINUE
  199. WRITE(IOIMP,*) 'An error was detected in memilu.eso'
  200. RETURN
  201. *
  202. * End of MEMILU
  203. *
  204. END
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.  
  212.  

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