Télécharger la plaquette Cast3M


fr en

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