Télécharger miltp3.eso

Retour à la liste

Numérotation des lignes :

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

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