Télécharger la plaquette Cast3M


fr en

Télécharger meiltp.eso

Retour à la liste

Numérotation des lignes : 

meiltp

  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.  
  64. -INC PPARAM
  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 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.       REAL*8 XNTLIG
  99. *      REAL*8 XPOND1,XPOND2
  100. *
  101. * Executable statements
  102. *
  103.       IF (XLFIL.LT.0.D0) THEN
  104.          WRITE(IOIMP,*) 'XLFIL.LT.0.D0 non acceptable'
  105.          GOTO 9999
  106.       ENDIF
  107.       PETIT=XZPREC**(0.75D0)
  108.       NBPIVP=0
  109.       NBPIVI=0
  110. c-----------------------------------------------------------------------
  111. c     initialize ju0 (points to next element to be added to alu,jlu)
  112. c     and pointer array.
  113. c-----------------------------------------------------------------------
  114.       JU0=NTT+2
  115.       JLU(1)=JU0
  116. *
  117. * Initialize permutation arrays
  118. *
  119.       DO ITT=1,NTT
  120.          JWORK(ITT)=0
  121.          IPERM(ITT)=ITT
  122.          JPERM(ITT)=ITT
  123.       ENDDO
  124. *
  125. * On boucle sur les ddl
  126. *
  127.       DO 1 ITT=1,NTT
  128. * Calcul de la norme de la ittème ligne
  129.          ISTRT=IA(ITT)
  130.          ISTOP=IA(ITT+1)-1
  131.          TNORM=ZERO
  132.          DO 11 IJA=ISTRT,ISTOP
  133.             TNORM=TNORM+ABS(A(IJA))
  134.  11      CONTINUE
  135.          TNORM=TNORM/DBLE(ISTOP-ISTRT+1)
  136.          IF (TNORM.LE.XPETIT) THEN
  137.             WRITE(IOIMP,*) 'La ligne ',ITT,' est nulle...'
  138.             GOTO 9999
  139.          ENDIF
  140. * On stocke la ittème ligne de A dans IWORK,RWORK en distinguant les
  141. * parties triangulaires inférieures et supérieurs.
  142. *
  143. * On peut le programmer autrement compte tenu du fait que les
  144. * colonnes sont ordonnées
  145. *
  146.          LENU=1
  147.          LENL=0
  148.          IWORK(ITT)=ITT
  149.          RWORK(ITT)=ZERO
  150.          JWORK(ITT)=ITT
  151.  
  152.          DO 12 IJA=ISTRT,ISTOP
  153. *!            COL=JA(IJA)
  154.             COL=JPERM(JA(IJA))
  155.             VAL=A(IJA)
  156.             IF (COL.LT.ITT) THEN
  157.                LENL=LENL+1
  158.                IWORK(LENL)=COL
  159.                RWORK(LENL)=VAL
  160.                JWORK(COL) =LENL
  161.             ELSEIF (COL.EQ.ITT) THEN
  162.                RWORK(ITT)=VAL
  163.             ELSE
  164.                LENU=LENU+1
  165.                JPOS=ITT+LENU-1
  166.                IWORK(JPOS)=COL
  167.                RWORK(JPOS)=VAL
  168.                JWORK(COL) =JPOS
  169.             ENDIF
  170.  12      CONTINUE
  171. * On a enlevé les stratégies bizarres
  172.          XNTLIG=DBLE(NNZA)/DBLE(NTT)
  173.          ILFIL=INT((XNTLIG-1.D0)*XLFIL/2.D0)
  174.          ILFILL=ILFIL
  175.          ILFILU=ILFIL
  176. * Elimination avec les lignes précédentes
  177.          JJL=0
  178.          LEN=0
  179.  13      CONTINUE
  180.          JJL=JJL+1
  181.          IF (JJL.LE.LENL) THEN
  182. c-----------------------------------------------------------------------
  183. c     in order to do the elimination in the correct order we must select
  184. c     the smallest column index among jw(k), k=jj+1, ..., lenl.
  185. c-----------------------------------------------------------------------
  186.             JROW=IWORK(JJL)
  187.             KJL=JJL
  188. c     determine smallest column index
  189.             DO 132 ILENL=JJL+1,LENL
  190.                IF (IWORK(ILENL).LT.JROW) THEN
  191.                   JROW=IWORK(ILENL)
  192.                   KJL=ILENL
  193.                ENDIF
  194.  132        CONTINUE
  195.             IF (KJL.NE.JJL) THEN
  196.                ITMP=IWORK(JJL)
  197.                IWORK(JJL)=IWORK(KJL)
  198.                IWORK(KJL)=ITMP
  199.                JWORK(JROW)=JJL
  200.                JWORK(ITMP)=KJL
  201.                RTMP=RWORK(JJL)
  202.                RWORK(JJL)=RWORK(KJL)
  203.                RWORK(KJL)=RTMP
  204.             ENDIF
  205. c     zero out element in row by setting jw(n+jrow) to zero.
  206.             JWORK(JROW)=0
  207. c     get the multiplier for row to be eliminated (jrow).
  208.             FACT=RWORK(JJL)*ALU(JROW)
  209.             IF (ABS(FACT).LE.XDTOL) THEN
  210.                GOTO 13
  211.             ENDIF
  212. c
  213. c     combine current row and row jrow
  214. c
  215.             DO 134 KLU=JU(JROW),JLU(JROW+1)-1
  216.                TERMLU=FACT*ALU(KLU)
  217. *!               COLLU =JLU(KLU)
  218.                COLLU =JPERM(JLU(KLU))
  219.                JPOS  =JWORK(COLLU)
  220.                IF (COLLU.GE.ITT) THEN
  221. c     dealing with upper part (including diagonal)
  222.                   IF (JPOS.EQ.0) THEN
  223. c     this is a fill-in element
  224.                      LENU=LENU+1
  225.                      IF (LENU.GT.NTT-ITT+1) THEN
  226.                         WRITE(IOIMP,*) 'pb. matrice U'
  227.                         GOTO 9999
  228.                      ENDIF
  229.                      IWORK(ITT+LENU-1)=COLLU
  230.                      JWORK(COLLU)=ITT+LENU-1
  231.                      RWORK(ITT+LENU-1)=-TERMLU
  232.                   ELSE
  233. c     this is not a fill-in element
  234.                      RWORK(JPOS)=RWORK(JPOS)-TERMLU
  235.                   ENDIF
  236.                ELSE
  237. c     dealing with lower part.
  238.                   IF (JPOS.EQ.0) THEN
  239. c     this is a fill-in element
  240.                      LENL=LENL+1
  241.                      IF (LENL.GT.ITT-1) THEN
  242.                         WRITE(IOIMP,*) 'pb. matrice L'
  243.                         WRITE(IOIMP,*) 'LENL=',LENL
  244.                         WRITE(IOIMP,*) 'ITT=',ITT
  245.                         GOTO 9999
  246.                      ENDIF
  247.                      IWORK(LENL) =COLLU
  248.                      JWORK(COLLU)=LENL
  249.                      RWORK(LENL) =-TERMLU
  250.                   ELSE
  251. c     this is not a fill-in element
  252.                      RWORK(JPOS)=RWORK(JPOS)-TERMLU
  253.                   ENDIF
  254.                ENDIF
  255.  134        CONTINUE
  256. c     store this pivot element -- (from left to right -- no danger of
  257. c     overlap with the working elements in L (pivots).
  258.             LEN=LEN+1
  259.             RWORK(LEN)=FACT
  260.             IWORK(LEN)=JROW
  261.             GOTO 13
  262.          ENDIF
  263. c     reset double-pointer to zero (U-part including diagonal)
  264.          DO 14 ILENU=1,LENU
  265.             JWORK(IWORK(ITT+ILENU-1))=0
  266.  14      CONTINUE
  267. c     update L-matrix
  268.          LENL=LEN
  269.          LEN=MIN(LENL,ILFILL)
  270. c     sort by quick-split
  271.          CALL QSPLIT(RWORK,IWORK,LENL,LEN,
  272.      $        IMPR,IRET)
  273. c     store L-part -- in original coordinates
  274.          DO 15 ILEN=1,LEN
  275. *!            JLU(JU0)=IWORK(ILEN)
  276.             ALU(JU0)=RWORK(ILEN)
  277.             JLU(JU0)=IPERM(IWORK(ILEN))
  278.             JU0=JU0+1
  279.  15      CONTINUE
  280. c     save pointer to beginning of row ii of U
  281.          JU(ITT)=JU0
  282. c     update U-matrix -- first apply dropping strategy
  283.          LEN=0
  284.          DO 16 ILEN=1,LENU-1
  285.             IF (ABS(RWORK(ITT+ILEN)).GT.XDTOL*TNORM) THEN
  286.                LEN=LEN+1
  287.                RWORK(ITT+LEN)=RWORK(ITT+ILEN)
  288.                IWORK(ITT+LEN)=IWORK(ITT+ILEN)
  289.             ENDIF
  290.  16      CONTINUE
  291.          LENU=LEN+1
  292. *!!!!!!!!
  293.          LEN=MIN(LENU,ILFILU+1)
  294. *!!!         LEN=MIN(LENU,ILFILU)
  295. c     sort by quick-split
  296.          CALL QSPLIT(RWORK(ITT+1),IWORK(ITT+1),LENU-1,LEN-1,
  297.      $        IMPR,IRET)
  298. *!
  299. *!    Determine next pivot
  300. *!
  301.          IMAX=ITT
  302.          XMAX=ABS(RWORK(IMAX))
  303.          XMAX0=XMAX
  304. * Je ne comprends pas bien cette ligne, je ne la mets pas
  305. *         ICUT=ITT-1+IRPIV-MOD(ITT-1,IRPIV)
  306. *         ICUT=ITT+IRPIV
  307.          DO ILEN=ITT+1,ITT+LEN-1
  308. *!!!         DO ILEN=ITT+1,ITT+LEN
  309.             XCOU=ABS(RWORK(ILEN))
  310.             IF (XCOU.GT.XMAX.AND.(XCOU*XSPIV).GT.XMAX0) THEN
  311.                IMAX=ILEN
  312.                XMAX=XCOU
  313.             ENDIF
  314.          ENDDO
  315. *
  316. * Echange des valeurs
  317. *
  318.          XTMP=RWORK(ITT)
  319.          RWORK(ITT)=RWORK(IMAX)
  320.          RWORK(IMAX)=XTMP
  321. *
  322. * Update de la numérotation ancienne
  323. *
  324.          JMAX=IWORK(IMAX)
  325.          ITMP=IPERM(ITT)
  326.          IPERM(ITT)=IPERM(JMAX)
  327.          IPERM(JMAX)=ITMP
  328. *
  329. * Update de la numérotation nouvelle
  330. *
  331.          JPERM(IPERM(ITT))=ITT
  332.          JPERM(IPERM(JMAX))=JMAX
  333. *
  334. c     copy U-Part in original coordinates
  335. *
  336.          DO 17 ILEN=ITT+1,ITT+LEN-1
  337. *!!!         DO 17 ILEN=ITT+1,ITT+LEN
  338.             JLU(JU0)=IPERM(IWORK(ILEN))
  339.             ALU(JU0)=RWORK(ILEN)
  340.             JU0=JU0+1
  341.  17      CONTINUE
  342. c     store inverse of diagonal element of u
  343.          VALPIV=RWORK(ITT)
  344.          IF (VALPIV.EQ.ZERO) THEN
  345.             WRITE(IOIMP,*) 'Pivot',ITT,' nul : ',
  346.      $           'le préconditionnement ILUTP est impossible.'
  347.             IRET=-2
  348.             GOTO 9999
  349.          ELSE
  350.             IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  351.                NBPIVP=NBPIVP+1
  352.             ENDIF
  353.             IF (VALPIV.LT.ZERO) THEN
  354.                NBPIVI=NBPIVI+1
  355.             ENDIF
  356.          ENDIF
  357.          ALU(ITT)=ONE/RWORK(ITT)
  358. c     update pointer to beginning of next row of U.
  359.          JLU(ITT+1)=JU0
  360.  1    CONTINUE
  361.       NNZLU=JLU(NTT+1)-1
  362. C
  363. C Permutation des colonnes de LU
  364. C
  365.       DO KLU=JLU(1),JLU(NTT+1)-1
  366.          JLU(KLU)=JPERM(JLU(KLU))
  367.       ENDDO
  368. C
  369. C Permutation des colonnes de A
  370. C
  371.       DO IJA=IA(1),IA(NTT+1)-1
  372.          JA(IJA)=JPERM(JA(IJA))
  373.       ENDDO
  374. *
  375. * Warning(s)
  376. *
  377.       IF (IMPR.GT.0) THEN
  378.          IF (NBPIVP.GT.0) THEN
  379.             WRITE(IOIMP,*) 'WARNING !',
  380.      $           ' meiltp :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  381.          ENDIF
  382. *!         IF (NBPIVI.GT.0) THEN
  383. *!            WRITE(IOIMP,*) 'WARNING !',
  384. *!     $           ' meiltp :',NBPIVI,' pivot(s) < 0'
  385. *!         ENDIF
  386.       ENDIF
  387. *
  388. * Normal termination
  389. *
  390.       IRET=0
  391.       RETURN
  392. *
  393. * Error handling
  394. *
  395.  9999 CONTINUE
  396.       IRET=1
  397.       WRITE(IOIMP,*) 'An error was detected in meiltp.eso'
  398.       RETURN
  399. *
  400. *     End of MEILTP
  401. *
  402.       END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  

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