Télécharger miltp2.eso

Retour à la liste

Numérotation des lignes :

  1. C MILTP2 SOURCE GOUNAND 05/02/16 21:16:16 5029
  2. SUBROUTINE MILTP2(NTT,NNZA,A,JA,IA,
  3. $ XLFIL,XDTOL,XSPIV,NNZMLU,
  4. $ ALU,JLU,JU,NNZLU,
  5. $ IWORK,JWORK,RWORK,LWORK,
  6. $ IPERM,JPERM,
  7. $ IMPR,IRET)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. IMPLICIT INTEGER (I-N)
  10. C***********************************************************************
  11. C NOM : MILTP2
  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 applique une stratégie de non-dropping lorsqu'on regarde les
  43. C termes qui seraient dans ILU(0)
  44. C signalé par *!!!!!!!!!!!!
  45. C
  46. C***********************************************************************
  47. C APPELES : QSPLIT
  48. C APPELE PAR : PRILUT
  49. C***********************************************************************
  50. C ENTREES : NTT,NNZA,A,JA,IA,XLFIL,XDTOL,NNZMLU,IVARI
  51. C ENTREES/SORTIES : IWORK,JWORK,RWORK (tableaux de travail)
  52. C SORTIES : ALU,JLU,JU,NNZLU
  53. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  54. C***********************************************************************
  55. C VERSION : v1, 08/04/04, version initiale
  56. C HISTORIQUE : v1, 08/04/04, création
  57. C HISTORIQUE :
  58. C HISTORIQUE :
  59. C***********************************************************************
  60. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  61. C en cas de modification de ce sous-programme afin de faciliter
  62. C la maintenance !
  63. C***********************************************************************
  64. -INC CCOPTIO
  65. -INC CCREEL
  66. * .. Scalar Arguments ..
  67. INTEGER NTT,NNZA,NNZMLU,NNZLU
  68. INTEGER 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. LOGICAL LWORK(NTT)
  78. INTEGER IPERM(NTT)
  79. INTEGER JPERM(NTT)
  80. *
  81. * .. Variables locales
  82. * .. Parameters
  83. REAL*8 ZERO ,ONE
  84. PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  85. * ..
  86. C Nombre de pivots petits
  87. INTEGER NBPIVP
  88. C Nombre de pivots inférieurs à 0
  89. INTEGER NBPIVI
  90. *
  91. INTEGER ILFIL,ILFILL,ILFILU
  92. INTEGER COL,COLLU
  93. INTEGER IJA,ILEN,ILENL,ILENU,ITT,ISTRT,ISTOP
  94. INTEGER JJL,JPOS,JROW,KJL,KLU
  95. INTEGER LEN,LENL,LENU
  96. INTEGER JU0,ITMP
  97. REAL*8 FACT,RTMP,TERMLU,TNORM,VAL,VALPIV,PETIT
  98. LOGICAL LTMP,LFACT
  99. REAL*8 XNTLIG
  100. * REAL*8 XPOND1,XPOND2
  101. *
  102. * Executable statements
  103. *
  104. IF (XLFIL.LT.1.D0) THEN
  105. WRITE(IOIMP,*) 'XLFIL.LT.1.D0 non acceptable'
  106. GOTO 9999
  107. ENDIF
  108. PETIT=XZPREC**(0.75D0)
  109. NBPIVP=0
  110. NBPIVI=0
  111. c-----------------------------------------------------------------------
  112. c initialize ju0 (points to next element to be added to alu,jlu)
  113. c and pointer array.
  114. c-----------------------------------------------------------------------
  115. JU0=NTT+2
  116. JLU(1)=JU0
  117. *
  118. * Initialize permutation arrays and lwork
  119. *
  120. DO ITT=1,NTT
  121. JWORK(ITT)=0
  122. IPERM(ITT)=ITT
  123. JPERM(ITT)=ITT
  124. *!!!!!!!!!!!!!!
  125. LWORK(ITT)=.FALSE.
  126. ENDDO
  127. *
  128. * On boucle sur les ddl
  129. *
  130. DO 1 ITT=1,NTT
  131. * Calcul de la norme de la ittème ligne
  132. ISTRT=IA(ITT)
  133. ISTOP=IA(ITT+1)-1
  134. TNORM=ZERO
  135. DO 11 IJA=ISTRT,ISTOP
  136. TNORM=TNORM+ABS(A(IJA))
  137. 11 CONTINUE
  138. TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  139. IF (TNORM.LE.XPETIT) THEN
  140. WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  141. GOTO 9999
  142. ENDIF
  143. * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
  144. * parties triangulaires inférieures et supérieurs.
  145. *
  146. * On peut le programmer autrement compte tenu du fait que les
  147. * colonnes sont ordonnées
  148. *
  149. LENU=1
  150. LENL=0
  151. IWORK(ITT)=ITT
  152. RWORK(ITT)=ZERO
  153. JWORK(ITT)=ITT
  154. LWORK(ITT)=.TRUE.
  155.  
  156. DO 12 IJA=ISTRT,ISTOP
  157. *! COL=JA(IJA)
  158. COL=JPERM(JA(IJA))
  159. VAL=A(IJA)
  160. IF (COL.LT.ITT) THEN
  161. LENL=LENL+1
  162. IWORK(LENL)=COL
  163. RWORK(LENL)=VAL
  164. *!!!!!!!!!!!!!!
  165. LWORK(LENL)=.TRUE.
  166. JWORK(COL) =LENL
  167. ELSEIF (COL.EQ.ITT) THEN
  168. RWORK(ITT)=VAL
  169. ELSE
  170. LENU=LENU+1
  171. JPOS=ITT+LENU-1
  172. IWORK(JPOS)=COL
  173. RWORK(JPOS)=VAL
  174. *!!!!!!!!!!!!!!
  175. LWORK(JPOS)=.TRUE.
  176. JWORK(COL) =JPOS
  177. ENDIF
  178. 12 CONTINUE
  179. * On a enlevé les stratégies bizarres
  180. XNTLIG=DBLE(NNZA)/DBLE(NTT)
  181. *!! ILFIL=INT(XNTLIG*XLFIL/2.D0)
  182. ILFIL=INT(XNTLIG*(XLFIL-1.D0)/2.D0)
  183. ILFILL=ILFIL
  184. ILFILU=ILFIL
  185. * Elimination avec les lignes précédentes
  186. JJL=0
  187. LEN=0
  188. 13 CONTINUE
  189. JJL=JJL+1
  190. IF (JJL.LE.LENL) THEN
  191. c-----------------------------------------------------------------------
  192. c in order to do the elimination in the correct order we must select
  193. c the smallest column index among jw(k), k=jj+1, ..., lenl.
  194. c-----------------------------------------------------------------------
  195. JROW=IWORK(JJL)
  196. KJL=JJL
  197. c determine smallest column index
  198. DO 132 ILENL=JJL+1,LENL
  199. IF (IWORK(ILENL).LT.JROW) THEN
  200. JROW=IWORK(ILENL)
  201. KJL=ILENL
  202. ENDIF
  203. 132 CONTINUE
  204. IF (KJL.NE.JJL) THEN
  205. ITMP=IWORK(JJL)
  206. IWORK(JJL)=IWORK(KJL)
  207. IWORK(KJL)=ITMP
  208. JWORK(JROW)=JJL
  209. JWORK(ITMP)=KJL
  210. RTMP=RWORK(JJL)
  211. RWORK(JJL)=RWORK(KJL)
  212. RWORK(KJL)=RTMP
  213. *!!!!!!!!!!!!!!
  214. LTMP=LWORK(JJL)
  215. LWORK(JJL)=LWORK(KJL)
  216. LWORK(KJL)=LTMP
  217. ENDIF
  218. c zero out element in row by setting jw(n+jrow) to zero.
  219. JWORK(JROW)=0
  220. c get the multiplier for row to be eliminated (jrow).
  221. FACT=RWORK(JJL)*ALU(JROW)
  222. LFACT=LWORK(JJL)
  223. *!!!!!!!!!!!!!!
  224. IF (ABS(FACT).LE.XDTOL.AND.(.NOT.LFACT)) THEN
  225. * IF (ABS(FACT).LE.XDTOL) THEN
  226. GOTO 13
  227. ENDIF
  228. c
  229. c combine current row and row jrow
  230. c
  231. DO 134 KLU=JU(JROW),JLU(JROW+1)-1
  232. TERMLU=FACT*ALU(KLU)
  233. *! COLLU =JLU(KLU)
  234. COLLU =JPERM(JLU(KLU))
  235. JPOS =JWORK(COLLU)
  236. IF (COLLU.GE.ITT) THEN
  237. c dealing with upper part (including diagonal)
  238. IF (JPOS.EQ.0) THEN
  239. c this is a fill-in element
  240. LENU=LENU+1
  241. IF (LENU.GT.NTT-ITT+1) THEN
  242. WRITE(IOIMP,*) 'pb. matrice U'
  243. GOTO 9999
  244. ENDIF
  245. IWORK(ITT+LENU-1)=COLLU
  246. JWORK(COLLU)=ITT+LENU-1
  247. RWORK(ITT+LENU-1)=-TERMLU
  248. LWORK(ITT+LENU-1)=.FALSE.
  249. ELSE
  250. c this is not a fill-in element
  251. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  252. ENDIF
  253. ELSE
  254. c dealing with lower part.
  255. IF (JPOS.EQ.0) THEN
  256. c this is a fill-in element
  257. LENL=LENL+1
  258. IF (LENL.GT.ITT-1) THEN
  259. WRITE(IOIMP,*) 'pb. matrice L'
  260. WRITE(IOIMP,*) 'LENL=',LENL
  261. WRITE(IOIMP,*) 'ITT=',ITT
  262. GOTO 9999
  263. ENDIF
  264. IWORK(LENL) =COLLU
  265. JWORK(COLLU)=LENL
  266. RWORK(LENL) =-TERMLU
  267. LWORK(LENL) =.FALSE.
  268. ELSE
  269. c this is not a fill-in element
  270. RWORK(JPOS)=RWORK(JPOS)-TERMLU
  271. ENDIF
  272. ENDIF
  273. 134 CONTINUE
  274. c store this pivot element -- (from left to right -- no danger of
  275. c overlap with the working elements in L (pivots).
  276. LEN=LEN+1
  277. * ?????????
  278. LWORK(LEN)=LFACT
  279. RWORK(LEN)=FACT
  280. IWORK(LEN)=JROW
  281. GOTO 13
  282. ENDIF
  283. c reset double-pointer to zero (U-part including diagonal)
  284. DO 14 ILENU=1,LENU
  285. JWORK(IWORK(ITT+ILENU-1))=0
  286. 14 CONTINUE
  287. c update L-matrix
  288. LENL=LEN
  289. LEN=MIN(LENL,ILFILL)
  290. c sort by quick-split
  291. CALL QSPLIT(RWORK,IWORK,LENL,LEN,
  292. $ IMPR,IRET)
  293. c store L-part -- in original coordinates
  294. DO 15 ILEN=1,LEN
  295. *! JLU(JU0)=IWORK(ILEN)
  296. ALU(JU0)=RWORK(ILEN)
  297. JLU(JU0)=IPERM(IWORK(ILEN))
  298. JU0=JU0+1
  299. *!!!!!!!!!!!!!!
  300. *! LWORK(ILEN)=.FALSE.
  301. 15 CONTINUE
  302. *!!!!!!!!!!!!!!
  303. DO ILEN=LEN+1,LENL
  304. IF (LWORK(ILEN)) THEN
  305. ALU(JU0)=RWORK(ILEN)
  306. JLU(JU0)=IPERM(IWORK(ILEN))
  307. JU0=JU0+1
  308. *!! LWORK(ILEN)=.FALSE.
  309. ENDIF
  310. ENDDO
  311.  
  312. c save pointer to beginning of row ii of U
  313. JU(ITT)=JU0
  314. c update U-matrix -- first apply dropping strategy
  315. LEN=0
  316. DO 16 ILEN=1,LENU-1
  317. IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
  318. LEN=LEN+1
  319. RWORK(ITT+LEN)=RWORK(ITT+ILEN)
  320. IWORK(ITT+LEN)=IWORK(ITT+ILEN)
  321. *!!!!!!!!!!!!!!
  322. LWORK(ITT+LEN)=LWORK(ITT+ILEN)
  323. *! LWORK(ITT+ILEN)=.FALSE.
  324. ENDIF
  325. 16 CONTINUE
  326. LENU=LEN+1
  327. *!!!!!!!!
  328. LEN=MIN(LENU,ILFILU+1)
  329. *!!! LEN=MIN(LENU,ILFILU)
  330. c sort by quick-split
  331. CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU-1,LEN-1,
  332. $ IMPR,IRET)
  333. *!
  334. *! Determine next pivot
  335. *!
  336. IMAX=ITT
  337. XMAX=ABS(RWORK(IMAX))
  338. XMAX0=XMAX
  339. * Je ne comprends pas bien cette ligne, je ne la mets pas
  340. * ICUT=ITT-1+IRPIV-MOD(ITT-1,IRPIV)
  341. * ICUT=ITT+IRPIV
  342. DO ILEN=ITT+1,ITT+LEN-1
  343. *!!! DO ILEN=ITT+1,ITT+LEN
  344. XCOU=ABS(RWORK(ILEN))
  345. IF (XCOU.GT.XMAX.AND.(XCOU*XSPIV).GT.XMAX0) THEN
  346. IMAX=ILEN
  347. XMAX=XCOU
  348. ENDIF
  349. ENDDO
  350. *
  351. * Echange des valeurs
  352. *
  353. XTMP=RWORK(ITT)
  354. RWORK(ITT)=RWORK(IMAX)
  355. RWORK(IMAX)=XTMP
  356. *!!!!!!!!!!!!!!
  357. LTMP=LWORK(ITT)
  358. LWORK(ITT)=LWORK(IMAX)
  359. LWORK(IMAX)=LTMP
  360. *
  361. * Update de la numérotation ancienne
  362. *
  363. JMAX=IWORK(IMAX)
  364. ITMP=IPERM(ITT)
  365. IPERM(ITT)=IPERM(JMAX)
  366. IPERM(JMAX)=ITMP
  367. *
  368. * Update de la numérotation nouvelle
  369. *
  370. JPERM(IPERM(ITT))=ITT
  371. JPERM(IPERM(JMAX))=JMAX
  372. *
  373. c copy U-Part in original coordinates
  374. *
  375. DO 17 ILEN=ITT+1,ITT+LEN-1
  376. *!!! DO 17 ILEN=ITT+1,ITT+LEN
  377. JLU(JU0)=IPERM(IWORK(ILEN))
  378. ALU(JU0)=RWORK(ILEN)
  379. JU0=JU0+1
  380. *!!!!!!!!!!!!!!
  381. *! LWORK(ILEN)=.FALSE.
  382. 17 CONTINUE
  383. *!!!!!!!!!!!!!!
  384. DO ILEN=ITT+LEN,ITT+LENU-1
  385. IF (LWORK(ILEN)) THEN
  386. ALU(JU0)=RWORK(ILEN)
  387. JLU(JU0)=IPERM(IWORK(ILEN))
  388. JU0=JU0+1
  389. *! LWORK(ILEN)=.FALSE.
  390. ENDIF
  391. ENDDO
  392. c store inverse of diagonal element of u
  393. VALPIV=RWORK(ITT)
  394. IF (VALPIV.EQ.ZERO) THEN
  395. WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
  396. $ 'le préconditionnement ILUTP est impossible.'
  397. IRET=-2
  398. GOTO 9999
  399. ELSE
  400. IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  401. NBPIVP=NBPIVP+1
  402. ENDIF
  403. IF (VALPIV.LT.ZERO) THEN
  404. NBPIVI=NBPIVI+1
  405. ENDIF
  406. ENDIF
  407. ALU(ITT)=ONE/RWORK(ITT)
  408. c update pointer to beginning of next row of U.
  409. JLU(ITT+1)=JU0
  410. 1 CONTINUE
  411. NNZLU=JLU(NTT+1)-1
  412. C
  413. C Permutation des colonnes de LU
  414. C
  415. DO KLU=JLU(1),JLU(NTT+1)-1
  416. JLU(KLU)=JPERM(JLU(KLU))
  417. ENDDO
  418. C
  419. C Permutation des colonnes de A
  420. C
  421. DO IJA=IA(1),IA(NTT+1)-1
  422. JA(IJA)=JPERM(JA(IJA))
  423. ENDDO
  424. *
  425. * Warning(s)
  426. *
  427. IF (IMPR.GT.0) THEN
  428. IF (NBPIVP.GT.0) THEN
  429. WRITE(IOIMP,*) 'WARNING !',
  430. $ ' miltp2 :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  431. ENDIF
  432. *! IF (NBPIVI.GT.0) THEN
  433. *! WRITE(IOIMP,*) 'WARNING !',
  434. *! $ ' miltp2 :',NBPIVI,' pivot(s) < 0'
  435. *! ENDIF
  436. ENDIF
  437. *
  438. * Normal termination
  439. *
  440. IRET=0
  441. RETURN
  442. *
  443. * Error handling
  444. *
  445. 9999 CONTINUE
  446. IRET=1
  447. WRITE(IOIMP,*) 'An error was detected in miltp2.eso'
  448. RETURN
  449. *
  450. * End of MILTP2
  451. *
  452. END
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  

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