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

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