Télécharger la plaquette Cast3M


fr en

Télécharger miltp2.eso

Retour à la liste

Numérotation des lignes : 

miltp2

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

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