Télécharger meiltp.eso

Retour à la liste

Numérotation des lignes :

  1. C MEILTP SOURCE GOUNAND 05/06/06 21:15:36 5078
  2. SUBROUTINE MEILTP(NTT,NNZA,A,JA,IA,
  3. $ XLFIL,XDTOL,XSPIV,NNZMLU,
  4. $ ALU,JLU,JU,NNZLU,
  5. $ IWORK,JWORK,RWORK,
  6. $ IPERM,JPERM,
  7. $ IMPR,IRET)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. IMPLICIT INTEGER (I-N)
  10. C***********************************************************************
  11. C NOM : MEILTP
  12. C DESCRIPTION :
  13. C Calcul du préconditionneur ILUTP d'une matrice Morse
  14. C ILUT : Incomplete LU factorization with dual truncation strategy
  15. C La matrice Morse est au format CSR (Column Sparse Row)
  16. C Le préconditionneur est au format MSR (Modified Sparse Row,
  17. C stockage de l'inverse de la diagonale) (ALU, JLU).
  18. C JU est un tableau supplémentaire pour le repérage des diagonales.
  19. C (voir la doc. de Sparskit version 2)
  20. C
  21. C LANGAGE : FORTRAN 77 (sauf E/S)
  22. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
  23. C mél : gounand@semt2.smts.cea.fr
  24. C REFERENCE (bibtex-like) :
  25. C @BOOK{templates,
  26. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  27. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  28. C H. Van der Vorst},
  29. C TITLE={Templates for the Solution of Linear Systems :
  30. C Building Blocks for Iterative Methods},
  31. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  32. C -> URL : http://www.netlib.org/templates/Templates.html
  33. C Sparskit : a basic tool kit for sparse matrix computations
  34. C Version 2 (Youcef Saad)
  35. C -> URL : http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  36. C REMARQUES :
  37. C Quelques parties pourraient être améliorés :
  38. C - l'algorithme de quicksplit (cf.qsplit.eso)
  39. C - la partie où l'on doit effectuer l'élimination dans le bon ordre
  40. C pourrait utiliser un algorithme de tri type quicksort
  41. C
  42. C On a tenté de corriger les bugs de la version Saad
  43. C signalés par *!
  44. C
  45. C***********************************************************************
  46. C APPELES : QSPLIT
  47. C APPELE PAR : PRILUT
  48. C***********************************************************************
  49. C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
  50. C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
  51. C SORTIES : ALU,JLU,JU,NNZLU
  52. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  53. C***********************************************************************
  54. C VERSION : v1, 08/04/04, version initiale
  55. C HISTORIQUE : v1, 08/04/04, 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,NNZA,NNZMLU,NNZLU
  67. INTEGER IMPR,IRET
  68. REAL*8 XDTOL,XLFIL
  69. * ..
  70. * .. Array Arguments ..
  71. REAL*8 A(NNZA),ALU(NNZMLU+1)
  72. INTEGER JA(NNZA),IA(NTT+1),JU(NTT),JLU(NNZMLU+1)
  73. INTEGER IWORK(NTT)
  74. INTEGER JWORK(NTT)
  75. REAL*8 RWORK(NTT)
  76. INTEGER IPERM(NTT)
  77. INTEGER JPERM(NTT)
  78. *
  79. * .. Variables locales
  80. * .. Parameters
  81. REAL*8 ZERO ,ONE
  82. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  83. * ..
  84. C Nombre de pivots petits
  85. INTEGER NBPIVP
  86. C Nombre de pivots inférieurs à 0
  87. INTEGER NBPIVI
  88. *
  89. INTEGER ILFIL,ILFILL,ILFILU
  90. INTEGER COL,COLLU
  91. INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
  92. INTEGER JJL,JPOS,JROW,KJL,KLU
  93. INTEGER LEN,LENL,LENU
  94. INTEGER JU0,ITMP
  95. REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
  96. REAL*8 XNTLIG
  97. * REAL*8 XPOND1,XPOND2
  98. *
  99. * Executable statements
  100. *
  101. IF (XLFIL.LT.0.D0) THEN
  102. WRITE(IOIMP,*) 'XLFIL.LT.0.D0 non acceptable'
  103. GOTO 9999
  104. ENDIF
  105. PETIT=XZPREC**(0.75D0)
  106. NBPIVP=0
  107. NBPIVI=0
  108. c-----------------------------------------------------------------------
  109. c initialize ju0 (points to next element to be added to alu,jlu)
  110. c and pointer array.
  111. c-----------------------------------------------------------------------
  112. JU0=NTT+2
  113. JLU(1)=JU0
  114. *
  115. * Initialize permutation arrays
  116. *
  117. DO ITT=1,NTT
  118. JWORK(ITT)=0
  119. IPERM(ITT)=ITT
  120. JPERM(ITT)=ITT
  121. ENDDO
  122. *
  123. * On boucle sur les ddl
  124. *
  125. DO 1 ITT=1,NTT
  126. * Calcul de la norme de la ittème ligne
  127. ISTRT=IA(ITT)
  128. ISTOP=IA(ITT+1)-1
  129. TNORM=ZERO
  130. DO 11 IJA=ISTRT,ISTOP
  131. TNORM=TNORM+ABS(A(IJA))
  132. 11 CONTINUE
  133. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  134. IF (TNORM.LE.XPETIT) THEN
  135. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  136. GOTO 9999
  137. ENDIF
  138. * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
  139. * parties triangulaires inférieures et supérieurs.
  140. *
  141. * On peut le programmer autrement compte tenu du fait que les
  142. * colonnes sont ordonnées
  143. *
  144. LENU=1
  145. LENL=0
  146. IWORK(ITT)=ITT
  147. RWORK(ITT)=ZERO
  148. JWORK(ITT)=ITT
  149.  
  150. DO 12 IJA=ISTRT,ISTOP
  151. *! COL=JA(IJA)
  152. COL=JPERM(JA(IJA))
  153. VAL=A(IJA)
  154. IF (COL.LT.ITT) THEN
  155. LENL=LENL+1
  156. IWORK(LENL)=COL
  157. RWORK(LENL)=VAL
  158. JWORK(COL) =LENL
  159. ELSEIF (COL.EQ.ITT) THEN
  160. RWORK(ITT)=VAL
  161. ELSE
  162. LENU=LENU+1
  163. JPOS=ITT+LENU-1
  164. IWORK(JPOS)=COL
  165. RWORK(JPOS)=VAL
  166. JWORK(COL) =JPOS
  167. ENDIF
  168. 12 CONTINUE
  169. * On a enlevé les stratégies bizarres
  170. XNTLIG=DBLE(NNZA)/DBLE(NTT)
  171. ILFIL=INT((XNTLIG-1.D0)*XLFIL/2.D0)
  172. ILFILL=ILFIL
  173. ILFILU=ILFIL
  174. * Elimination avec les lignes précédentes
  175. JJL=0
  176. LEN=0
  177. 13 CONTINUE
  178. JJL=JJL+1
  179. IF (JJL.LE.LENL) THEN
  180. c-----------------------------------------------------------------------
  181. c in order to do the elimination in the correct order we must select
  182. c the smallest column index among jw(k), k=jj+1, ..., lenl.
  183. c-----------------------------------------------------------------------
  184. JROW=IWORK(JJL)
  185. KJL=JJL
  186. c determine smallest column index
  187. DO 132 ILENL=JJL+1,LENL
  188. IF (IWORK(ILENL).LT.JROW) THEN
  189. JROW=IWORK(ILENL)
  190. KJL=ILENL
  191. ENDIF
  192. 132 CONTINUE
  193. IF (KJL.NE.JJL) THEN
  194. ITMP=IWORK(JJL)
  195. IWORK(JJL)=IWORK(KJL)
  196. IWORK(KJL)=ITMP
  197. JWORK(JROW)=JJL
  198. JWORK(ITMP)=KJL
  199. RTMP=RWORK(JJL)
  200. RWORK(JJL)=RWORK(KJL)
  201. RWORK(KJL)=RTMP
  202. ENDIF
  203. c zero out element in row by setting jw(n+jrow) to zero.
  204. JWORK(JROW)=0
  205. c get the multiplier for row to be eliminated (jrow).
  206. FACT=RWORK(JJL)*ALU(JROW)
  207. IF (ABS(FACT).LE.XDTOL) THEN
  208. GOTO 13
  209. ENDIF
  210. c
  211. c combine current row and row jrow
  212. c
  213. DO 134 KLU=JU(JROW),JLU(JROW+1)-1
  214. TERMLU=FACT*ALU(KLU)
  215. *! COLLU =JLU(KLU)
  216. COLLU =JPERM(JLU(KLU))
  217. JPOS =JWORK(COLLU)
  218. IF (COLLU.GE.ITT) THEN
  219. c dealing with upper part (including diagonal)
  220. IF (JPOS.EQ.0) THEN
  221. c this is a fill-in element
  222. LENU=LENU+1
  223. IF (LENU.GT.NTT-ITT+1) THEN
  224. WRITE(IOIMP,*) 'pb. matrice U'
  225. GOTO 9999
  226. ENDIF
  227. IWORK(ITT+LENU-1)=COLLU
  228. JWORK(COLLU)=ITT+LENU-1
  229. RWORK(ITT+LENU-1)=-TERMLU
  230. ELSE
  231. c this is not a fill-in element
  232. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  233. ENDIF
  234. ELSE
  235. c dealing with lower part.
  236. IF (JPOS.EQ.0) THEN
  237. c this is a fill-in element
  238. LENL=LENL+1
  239. IF (LENL.GT.ITT-1) THEN
  240. WRITE(IOIMP,*) 'pb. matrice L'
  241. WRITE(IOIMP,*) 'LENL=',LENL
  242. WRITE(IOIMP,*) 'ITT=',ITT
  243. GOTO 9999
  244. ENDIF
  245. IWORK(LENL) =COLLU
  246. JWORK(COLLU)=LENL
  247. RWORK(LENL) =-TERMLU
  248. ELSE
  249. c this is not a fill-in element
  250. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  251. ENDIF
  252. ENDIF
  253. 134 CONTINUE
  254. c store this pivot element -- (from left to right -- no danger of
  255. c overlap with the working elements in L (pivots).
  256. LEN=LEN+1
  257. RWORK(LEN)=FACT
  258. IWORK(LEN)=JROW
  259. GOTO 13
  260. ENDIF
  261. c reset double-pointer to zero (U-part including diagonal)
  262. DO 14 ILENU=1,LENU
  263. JWORK(IWORK(ITT+ILENU-1))=0
  264. 14 CONTINUE
  265. c update L-matrix
  266. LENL=LEN
  267. LEN=MIN(LENL,ILFILL)
  268. c sort by quick-split
  269. CALL QSPLIT(RWORK,IWORK,LENL,LEN,
  270. $ IMPR,IRET)
  271. c store L-part -- in original coordinates
  272. DO 15 ILEN=1,LEN
  273. *! JLU(JU0)=IWORK(ILEN)
  274. ALU(JU0)=RWORK(ILEN)
  275. JLU(JU0)=IPERM(IWORK(ILEN))
  276. JU0=JU0+1
  277. 15 CONTINUE
  278. c save pointer to beginning of row ii of U
  279. JU(ITT)=JU0
  280. c update U-matrix -- first apply dropping strategy
  281. LEN=0
  282. DO 16 ILEN=1,LENU-1
  283. IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
  284. LEN=LEN+1
  285. RWORK(ITT+LEN)=RWORK(ITT+ILEN)
  286. IWORK(ITT+LEN)=IWORK(ITT+ILEN)
  287. ENDIF
  288. 16 CONTINUE
  289. LENU=LEN+1
  290. *!!!!!!!!
  291. LEN=MIN(LENU,ILFILU+1)
  292. *!!! LEN=MIN(LENU,ILFILU)
  293. c sort by quick-split
  294. CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU-1,LEN-1,
  295. $ IMPR,IRET)
  296. *!
  297. *! Determine next pivot
  298. *!
  299. IMAX=ITT
  300. XMAX=ABS(RWORK(IMAX))
  301. XMAX0=XMAX
  302. * Je ne comprends pas bien cette ligne, je ne la mets pas
  303. * ICUT=ITT-1+IRPIV-MOD(ITT-1,IRPIV)
  304. * ICUT=ITT+IRPIV
  305. DO ILEN=ITT+1,ITT+LEN-1
  306. *!!! DO ILEN=ITT+1,ITT+LEN
  307. XCOU=ABS(RWORK(ILEN))
  308. IF (XCOU.GT.XMAX.AND.(XCOU*XSPIV).GT.XMAX0) THEN
  309. IMAX=ILEN
  310. XMAX=XCOU
  311. ENDIF
  312. ENDDO
  313. *
  314. * Echange des valeurs
  315. *
  316. XTMP=RWORK(ITT)
  317. RWORK(ITT)=RWORK(IMAX)
  318. RWORK(IMAX)=XTMP
  319. *
  320. * Update de la numérotation ancienne
  321. *
  322. JMAX=IWORK(IMAX)
  323. ITMP=IPERM(ITT)
  324. IPERM(ITT)=IPERM(JMAX)
  325. IPERM(JMAX)=ITMP
  326. *
  327. * Update de la numérotation nouvelle
  328. *
  329. JPERM(IPERM(ITT))=ITT
  330. JPERM(IPERM(JMAX))=JMAX
  331. *
  332. c copy U-Part in original coordinates
  333. *
  334. DO 17 ILEN=ITT+1,ITT+LEN-1
  335. *!!! DO 17 ILEN=ITT+1,ITT+LEN
  336. JLU(JU0)=IPERM(IWORK(ILEN))
  337. ALU(JU0)=RWORK(ILEN)
  338. JU0=JU0+1
  339. 17 CONTINUE
  340. c store inverse of diagonal element of u
  341. VALPIV=RWORK(ITT)
  342. IF (VALPIV.EQ.ZERO) THEN
  343. WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
  344. $ 'le préconditionnement ILUTP est impossible.'
  345. IRET=-2
  346. GOTO 9999
  347. ELSE
  348. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  349. NBPIVP=NBPIVP+1
  350. ENDIF
  351. IF (VALPIV.LT.ZERO) THEN
  352. NBPIVI=NBPIVI+1
  353. ENDIF
  354. ENDIF
  355. ALU(ITT)=ONE/RWORK(ITT)
  356. c update pointer to beginning of next row of U.
  357. JLU(ITT+1)=JU0
  358. 1 CONTINUE
  359. NNZLU=JLU(NTT+1)-1
  360. C
  361. C Permutation des colonnes de LU
  362. C
  363. DO KLU=JLU(1),JLU(NTT+1)-1
  364. JLU(KLU)=JPERM(JLU(KLU))
  365. ENDDO
  366. C
  367. C Permutation des colonnes de A
  368. C
  369. DO IJA=IA(1),IA(NTT+1)-1
  370. JA(IJA)=JPERM(JA(IJA))
  371. ENDDO
  372. *
  373. * Warning(s)
  374. *
  375. IF (IMPR.GT.0) THEN
  376. IF (NBPIVP.GT.0) THEN
  377. WRITE(IOIMP,*) 'WARNING !',
  378. $ ' meiltp :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  379. ENDIF
  380. *! IF (NBPIVI.GT.0) THEN
  381. *! WRITE(IOIMP,*) 'WARNING !',
  382. *! $ ' meiltp :',NBPIVI,' pivot(s) < 0'
  383. *! ENDIF
  384. ENDIF
  385. *
  386. * Normal termination
  387. *
  388. IRET=0
  389. RETURN
  390. *
  391. * Error handling
  392. *
  393. 9999 CONTINUE
  394. IRET=1
  395. WRITE(IOIMP,*) 'An error was detected in meiltp.eso'
  396. RETURN
  397. *
  398. * End of MEILTP
  399. *
  400. END
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  

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