Télécharger meilut.eso

Retour à la liste

Numérotation des lignes :

meilut
  1. C MEILUT SOURCE GOUNAND 05/06/06 21:15:39 5078
  2. SUBROUTINE MEILUT(NTT,NNZA,A,JA,IA,
  3. $ XLFIL,XDTOL,NNZMLU,
  4. $ ALU,JLU,JU,NNZLU,
  5. $ IWORK,JWORK,RWORK,
  6. $ IVARI,IMPR,IRET)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. IMPLICIT INTEGER (I-N)
  9. C***********************************************************************
  10. C NOM : MEILUT
  11. C DESCRIPTION :
  12. C Calcul du préconditionneur ILUT d'une matrice Morse (IVARI=0)
  13. C ou d'une variante qui remplit mieux la mémoire et fonctionne mieux
  14. C quelquefois (IVARI=1).
  15. C ILUT : Incomplete LU factorization with dual truncation strategy
  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 REMARQUES :
  38. C Quelques parties pourraient être améliorés :
  39. C - l'algorithme de quicksplit (cf.qsplit.eso)
  40. C - la partie où l'on doit effectuer l'élimination dans le bon ordre
  41. C pourrait utiliser un algorithme de tri type quicksort
  42. C
  43. C***********************************************************************
  44. C APPELES : QSPLIT
  45. C APPELE PAR : PRILUT
  46. C***********************************************************************
  47. C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
  48. C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
  49. C SORTIES : ALU,JLU,JU,NNZLU
  50. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  51. C***********************************************************************
  52. C VERSION : v1, 20/03/00, réécriture totale
  53. C VERSION : v1, 23/02/00, version initiale
  54. C HISTORIQUE : v1, 23/02/00, création
  55. C HISTORIQUE :
  56. C HISTORIQUE :
  57. C***********************************************************************
  58. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  59. C en cas de modification de ce sous-programme afin de faciliter
  60. C la maintenance !
  61. C***********************************************************************
  62.  
  63. -INC PPARAM
  64. -INC CCOPTIO
  65. -INC CCREEL
  66. * .. Scalar Arguments ..
  67. INTEGER NTT,NNZA,NNZMLU,NNZLU
  68. INTEGER IVARI,IMPR,IRET
  69. REAL*8 XDTOL,XLFIL
  70. * ..
  71. * .. Array Arguments ..
  72. REAL*8 A(NNZA),ALU(NNZMLU+1)
  73. INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
  74. INTEGER IWORK(NTT)
  75. INTEGER JWORK(NTT)
  76. REAL*8 RWORK(NTT)
  77. *
  78. * .. Variables locales
  79. * .. Parameters
  80. REAL*8 ZERO ,ONE
  81. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  82. * ..
  83. C Nombre de pivots petits
  84. INTEGER NBPIVP
  85. C Nombre de pivots inférieurs à 0
  86. INTEGER NBPIVI
  87. *
  88. INTEGER ILFIL,ILFILL,ILFILU
  89. INTEGER COL,COLLU
  90. INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
  91. INTEGER JJL,JPOS,JROW,KJL,KLU
  92. INTEGER LEN,LENL,LENU
  93. INTEGER JU0,ITMP
  94. REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
  95. REAL*8 XNTLIG
  96. * REAL*8 XPOND1,XPOND2
  97. *
  98. * Executable statements
  99. *
  100. IF (XLFIL.LT.0.D0) THEN
  101. WRITE(IOIMP,*) 'XLFIL.LT.0.D0 non acceptable'
  102. GOTO 9999
  103. ENDIF
  104. PETIT=XZPREC**(0.75D0)
  105. NBPIVP=0
  106. NBPIVI=0
  107. c-----------------------------------------------------------------------
  108. c initialize ju0 (points to next element to be added to alu,jlu)
  109. c and pointer array.
  110. c-----------------------------------------------------------------------
  111. JU0=NTT+2
  112. JLU(1)=JU0
  113. *
  114. * On boucle sur les ddl
  115. *
  116. DO 1 ITT=1,NTT
  117. * Calcul de la norme de la ittème ligne
  118. ISTRT=IA(ITT)
  119. ISTOP=IA(ITT+1)-1
  120. TNORM=ZERO
  121. DO 11 IJA=ISTRT,ISTOP
  122. TNORM=TNORM+ABS(A(IJA))
  123. 11 CONTINUE
  124. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  125. IF (TNORM.LE.XPETIT) THEN
  126. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  127. GOTO 9999
  128. ENDIF
  129. * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
  130. * parties triangulaires inférieures et supérieurs.
  131. *
  132. * On peut le programmer autrement compte tenu du fait que les
  133. * colonnes sont ordonnées
  134. *
  135. LENL=0
  136. LENU=0
  137. DO 12 IJA=ISTRT,ISTOP
  138. COL=JA(IJA)
  139. VAL=A(IJA)
  140. IF (COL.LT.ITT) THEN
  141. LENL=LENL+1
  142. IWORK(LENL)=COL
  143. RWORK(LENL)=VAL
  144. JWORK(COL) =LENL
  145. ELSEIF (COL.EQ.ITT) THEN
  146. IWORK(ITT)=ITT
  147. RWORK(ITT)=VAL
  148. JWORK(ITT)=ITT
  149. ELSE
  150. LENU=LENU+1
  151. IWORK(ITT+LENU)=COL
  152. RWORK(ITT+LENU)=VAL
  153. JWORK(COL) =ITT+LENU
  154. ENDIF
  155. 12 CONTINUE
  156. IF (IVARI.EQ.1) THEN
  157. * Stratégie : on remplit autant que possible
  158. * normale meilug
  159. XNTLIG=DBLE(NNZMLU-JU0+1)/DBLE(NTT-ITT+1)
  160. * on retire la diagonale
  161. * ILFIL=INT(XNTLIG/2.D0)
  162. ILFIL=INT((XNTLIG-1.D0)/2.D0)
  163. * Stratégie meilug pondérée globalement
  164. * XNTLIG=(DBLE(NNZMLU-JU0+1)/DBLE(NTT-ITT+1))
  165. * $ *(1.5D0-0.5D0*(DBLE(ITT-1)/DBLE(NTT-1)))
  166. * ILFIL=INT(XNTLIG/2.D0)
  167. * Stratégie normale (idem meilut)
  168. ELSEIF (IVARI.EQ.0) THEN
  169. XNTLIG=DBLE(NNZA)/DBLE(NTT)
  170. * Pour faire comme dans meiltp
  171. * ILFIL=INT(XNTLIG*XLFIL/2.D0)
  172. ILFIL=INT((XNTLIG-1.D0)*XLFIL/2.D0)
  173. ELSE
  174. WRITE(IOIMP,*) 'Erreur de programmation IVARI=',IVARI
  175. ENDIF
  176. * On a enlevé les :
  177. * Corrections pour les extrémités de la matrice
  178. * IF (ILFIL.GT.ITT-1) THEN
  179. * ILFILL=ITT-1
  180. * ILFILU=(2*ILFIL)-ILFILL
  181. * ELSEIF (ILFIL.GT.(NTT-ITT)) THEN
  182. * ILFILU=NTT-ITT
  183. * ILFILL=(2*ILFIL)-ILFILU
  184. * ELSE
  185. * ILFILL=ILFIL
  186. * ILFILU=ILFIL
  187. * ENDIF
  188. * WRITE(IOIMP,*) 'itt,ilfil,ilfill,ilfilu=',ITT,ILFIL,ILFILL
  189. * $ ,ILFILU
  190. * Stratégie normale
  191. ILFILL=ILFIL
  192. ILFILU=ILFIL
  193. * Stratégie pondérée localement
  194. * XPOND1=DBLE(LENL)/DBLE(LENL+LENU)
  195. * XPOND2=DBLE(LENU)/DBLE(LENL+LENU)
  196. * ILFILL=INT(XPOND1*DBLE(ILFIL*2))
  197. * ILFILU=(ILFIL*2)-ILFILL
  198. * Elimination avec les lignes précédentes
  199. JJL=0
  200. LEN=0
  201. 13 CONTINUE
  202. JJL=JJL+1
  203. IF (JJL.LE.LENL) THEN
  204. c-----------------------------------------------------------------------
  205. c in order to do the elimination in the correct order we must select
  206. c the smallest column index among jw(k), k=jj+1, ..., lenl.
  207. c-----------------------------------------------------------------------
  208. JROW=IWORK(JJL)
  209. KJL=JJL
  210. c determine smallest column index
  211. DO 132 ILENL=JJL+1,LENL
  212. IF (IWORK(ILENL).LT.JROW) THEN
  213. JROW=IWORK(ILENL)
  214. KJL=ILENL
  215. ENDIF
  216. 132 CONTINUE
  217. IF (KJL.NE.JJL) THEN
  218. ITMP=IWORK(JJL)
  219. IWORK(JJL)=IWORK(KJL)
  220. IWORK(KJL)=ITMP
  221. JWORK(JROW)=JJL
  222. JWORK(ITMP)=KJL
  223. RTMP=RWORK(JJL)
  224. RWORK(JJL)=RWORK(KJL)
  225. RWORK(KJL)=RTMP
  226. ENDIF
  227. c zero out element in row by setting jw(n+jrow) to zero.
  228. JWORK(JROW)=0
  229. c get the multiplier for row to be eliminated (jrow).
  230. FACT=RWORK(JJL)*ALU(JROW)
  231. IF (ABS(FACT).LE.XDTOL) THEN
  232. GOTO 13
  233. ENDIF
  234. c
  235. c combine current row and row jrow
  236. c
  237. DO 134 KLU=JU(JROW),JLU(JROW+1)-1
  238. TERMLU=FACT*ALU(KLU)
  239. COLLU =JLU(KLU)
  240. JPOS =JWORK(COLLU)
  241. IF (COLLU.GE.ITT) THEN
  242. c dealing with upper part (including diagonal)
  243. IF (JPOS.EQ.0) THEN
  244. c this is a fill-in element
  245. LENU=LENU+1
  246. IF (LENU.GT.NTT-ITT) THEN
  247. WRITE(IOIMP,*) 'pb. matrice U'
  248. GOTO 9999
  249. ENDIF
  250. IWORK(ITT+LENU)=COLLU
  251. JWORK(COLLU)=ITT+LENU
  252. RWORK(ITT+LENU)=-TERMLU
  253. ELSE
  254. c this is not a fill-in element
  255. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  256. ENDIF
  257. ELSE
  258. c dealing with lower part.
  259. IF (JPOS.EQ.0) THEN
  260. c this is a fill-in element
  261. LENL=LENL+1
  262. IF (LENL.GT.ITT-1) THEN
  263. WRITE(IOIMP,*) 'pb. matrice L'
  264. GOTO 9999
  265. ENDIF
  266. IWORK(LENL) =COLLU
  267. JWORK(COLLU)=LENL
  268. RWORK(LENL) =-TERMLU
  269. ELSE
  270. c this is not a fill-in element
  271. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  272. ENDIF
  273. ENDIF
  274. 134 CONTINUE
  275. c store this pivot element -- (from left to right -- no danger of
  276. c overlap with the working elements in L (pivots).
  277. LEN=LEN+1
  278. IWORK(LEN)=JROW
  279. RWORK(LEN)=FACT
  280. GOTO 13
  281. ENDIF
  282. c reset double-pointer to zero (U-part including diagonal)
  283. DO 14 ILENU=0,LENU
  284. JWORK(IWORK(ITT+ILENU))=0
  285. 14 CONTINUE
  286. c update L-matrix
  287. LENL=LEN
  288. LEN=MIN(LENL,ILFILL)
  289. c sort by quick-split
  290. CALL QSPLIT(RWORK,IWORK,LENL,LEN,
  291. $ IMPR,IRET)
  292. c store L-part
  293. DO 15 ILEN=1,LEN
  294. JLU(JU0)=IWORK(ILEN)
  295. ALU(JU0)=RWORK(ILEN)
  296. JU0=JU0+1
  297. 15 CONTINUE
  298. c save pointer to beginning of row ii of U
  299. JU(ITT)=JU0
  300. c update U-matrix -- first apply dropping strategy
  301. LEN=0
  302. DO 16 ILEN=1,LENU
  303. IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
  304. LEN=LEN+1
  305. IWORK(ITT+LEN)=IWORK(ITT+ILEN)
  306. RWORK(ITT+LEN)=RWORK(ITT+ILEN)
  307. ENDIF
  308. 16 CONTINUE
  309. LENU=LEN
  310. LEN=MIN(LENU,ILFILU)
  311. c sort by quick-split
  312. CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU,LEN,
  313. $ IMPR,IRET)
  314. c copy
  315. DO 17 ILEN=ITT+1,ITT+LEN
  316. JLU(JU0)=IWORK(ILEN)
  317. ALU(JU0)=RWORK(ILEN)
  318. JU0=JU0+1
  319. 17 CONTINUE
  320. c store inverse of diagonal element of u
  321. VALPIV=RWORK(ITT)
  322. IF (VALPIV.EQ.ZERO) THEN
  323. WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  324. $ 'le préconditionnement ILUT est impossible.'
  325. IRET=-2
  326. GOTO 9999
  327. ELSE
  328. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  329. NBPIVP=NBPIVP+1
  330. ENDIF
  331. IF (VALPIV.LT.ZERO) THEN
  332. NBPIVI=NBPIVI+1
  333. ENDIF
  334. ENDIF
  335. ALU(ITT)=ONE/RWORK(ITT)
  336. c update pointer to beginning of next row of U.
  337. JLU(ITT+1)=JU0
  338. 1 CONTINUE
  339. NNZLU=JLU(NTT+1)-1
  340. *
  341. * Warning(s)
  342. *
  343. IF (IMPR.GT.0) THEN
  344. IF (NBPIVP.GT.0) THEN
  345. WRITE(IOIMP,*) 'WARNING !',
  346. $ ' meilut :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  347. ENDIF
  348. *! IF (NBPIVI.GT.0) THEN
  349. *! WRITE(IOIMP,*) 'WARNING !',
  350. *! $ ' meilut :',NBPIVI,' pivot(s) < 0'
  351. *! ENDIF
  352. ENDIF
  353. *
  354. * Normal termination
  355. *
  356. IRET=0
  357. RETURN
  358. *
  359. * Error handling
  360. *
  361. 9999 CONTINUE
  362. IRET=1
  363. WRITE(IOIMP,*) 'An error was detected in meilut.eso'
  364. RETURN
  365. *
  366. * End of MEILUT
  367. *
  368. END
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  

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