Télécharger la plaquette Cast3M


fr en

Télécharger meilup.eso

Retour à la liste

Numérotation des lignes : 

meilup

  1. C MEILUP    SOURCE    GOUNAND   11/07/21    21:15:25     7046
  2.       SUBROUTINE MEILUP(NTT,NNZ,A,JA,IA,
  3.      $     ALU,JLU,JU,
  4.      $     IWORK,
  5.      $     RXILUP,IFILTR,
  6.      $     IMPR,IRET)
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8.       IMPLICIT INTEGER (I-N)
  9. C***********************************************************************
  10. C NOM         : MEILUP
  11. C DESCRIPTION :
  12. C     Calcul du préconditionneur ILU(0) d'une matrice Morse.
  13. C     modifié PV
  14. C     ILU(0) : Incomplete LU factorization of level 0
  15. C              appelée aussi Choleski ou Crout incomplet
  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
  38. C***********************************************************************
  39. C APPELES          : -
  40. C APPELE PAR       : PRILU0
  41. C***********************************************************************
  42. C ENTREES            : NTT, NNZ, IA, JA, A
  43. C ENTREES/SORTIES    : IWORK (tableau de travail)
  44. C SORTIES            : JU, JLU, ALU
  45. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  46. C***********************************************************************
  47. C VERSION    : v1.1, 21/03/00
  48. C HISTORIQUE : v1.1, 22/03/00, gounand
  49. C Amélioration de la détection des pivots nulles.
  50. C HISTORIQUE : v1,   20/12/99, création
  51. C HISTORIQUE :
  52. C HISTORIQUE :
  53. C***********************************************************************
  54. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  55. C en cas de modification de ce sous-programme afin de faciliter
  56. C la maintenance !
  57. C***********************************************************************
  58.  
  59. -INC PPARAM
  60. -INC CCOPTIO
  61. -INC CCREEL
  62. *     .. Scalar Arguments ..
  63.       INTEGER  NTT,NNZ
  64.       INTEGER  IMPR,IRET
  65. *     ..
  66. *     .. Array Arguments ..
  67.       REAL*8  A(NNZ),ALU(NNZ+1)
  68.       INTEGER JA(NNZ),IA(NTT+1),JU(NTT),JLU(NNZ+1)
  69.       INTEGER IWORK(NTT)
  70. *
  71. *     .. Variables locales
  72. *     .. Parameters
  73.       REAL*8     ZERO      ,ONE
  74.       PARAMETER (ZERO=0.0D0,ONE=1.0D0)
  75. *     ..
  76. C     Nombre de pivots petits
  77.       INTEGER NBPIVP
  78. C     Nombre de pivots inférieurs à 0
  79.       INTEGER NBPIVI
  80.       INTEGER ITT,JNZ,KNZ,JCOL,JROW
  81.       INTEGER JSTRT,JMID,JSTOP,JU0,JW
  82.       INTEGER ISTRT,ISTOP
  83.       REAL*8 TL,VALPIV,TNORM,PETIT,XOPIV,XNPIV
  84.  
  85. *
  86. * Executable statements
  87. *
  88.       PETIT=XZPREC**(0.75D0)
  89.       NBPIVP=0
  90.       NBPIVI=0
  91.       JU0=NTT+2
  92.       JLU(1)=JU0
  93. *
  94. * On boucle sur les ddl
  95. *
  96.       DO 1 ITT=1,NTT
  97.          JSTRT=JU0
  98. *
  99. * Création de la ligne ITT pour L et U
  100. *
  101.          ISTRT=IA(ITT)
  102.          ISTOP=IA(ITT+1)-1
  103.          TNORM=ZERO
  104.          XOPIV=ZERO
  105.          DO 12 JNZ=ISTRT,ISTOP
  106. *
  107. * Copie de la ligne ITT du format CSR au format MSR
  108. * et calcul de la norme de la ligne
  109. *
  110.             JCOL=JA(JNZ)
  111.             TNORM=TNORM+ABS(A(JNZ))
  112.             IF (JCOL.EQ.ITT) THEN
  113.                XOPIV      =A(JNZ)
  114.                ALU(ITT)   =A(JNZ)
  115.                IWORK(JCOL)=ITT
  116.                JU(ITT)    =JU0
  117.             ELSE
  118.                ALU(JU0)   =A(JNZ)
  119.                JLU(JU0)   =JA(JNZ)
  120.                IWORK(JCOL)=JU0
  121.                JU0 = JU0+1
  122.             ENDIF
  123.  12      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.          JLU(ITT+1)=JU0
  130.          JSTOP=JU0-1
  131.          JMID=JU(ITT)-1
  132.          DO 14 JNZ=JSTRT,JMID
  133.             JROW=JLU(JNZ)
  134.             TL=ALU(JNZ)*ALU(JROW)
  135.             ALU(JNZ)=TL
  136. *
  137. * On effectue la combinaison linéaire
  138. *
  139.             DO 142 KNZ=JU(JROW),JLU(JROW+1)-1
  140.                JW=IWORK(JLU(KNZ))
  141.                IF (JW.NE.0) THEN
  142.                   ALU(JW)=ALU(JW)-TL*ALU(KNZ)
  143.                ENDIF
  144.  142        CONTINUE
  145.  14      CONTINUE
  146. *
  147. *  On s'occupe de l'élément diagonal
  148. *
  149. *ILU(0)         VALPIV=ALU(ITT)
  150.          IF (IFILTR.EQ.1) THEN
  151. *New deb
  152.             XNPIV=ALU(ITT)
  153.             VALPIV=XNPIV
  154.             IF (ABS(XOPIV).GT.XPETIT) THEN
  155.                IF (XNPIV.LT.(ABS(XOPIV)*RXILUP)) THEN
  156.                   VALPIV=XOPIV
  157.                ENDIF
  158.             ENDIF
  159.          ENDIF
  160. *New fin
  161.          IF (IFILTR.EQ.0) THEN
  162.             VALPIV=XOPIV
  163.          ENDIF
  164. *
  165.          IF (VALPIV.EQ.ZERO) THEN
  166.             WRITE(IOIMP,*) 'Pivot',ITT,'nul : ',
  167.      $           'le préconditionnement ILU(0)-PV est impossible.'
  168.             IRET=-2
  169.             GOTO 9999
  170.          ELSE
  171.             IF (ABS(VALPIV).LT.PETIT*TNORM) THEN
  172.                NBPIVP=NBPIVP+1
  173.             ENDIF
  174.             IF (VALPIV.LT.ZERO) THEN
  175.                NBPIVI=NBPIVI+1
  176.             ENDIF
  177.          ENDIF
  178.          ALU(ITT)=ONE/VALPIV
  179. *
  180. * On remet le tableau de travail à zéro
  181. *
  182.          IWORK(ITT)=0
  183.          DO 16 JNZ=JSTRT,JSTOP
  184.             IWORK(JLU(JNZ))=0
  185.  16      CONTINUE
  186.  1    CONTINUE
  187. *
  188. * Warning(s)
  189. *
  190.       IF (IMPR.GT.0) THEN
  191.          IF (NBPIVP.GT.0) THEN
  192.             WRITE(IOIMP,*) 'WARNING !',
  193.      $           ' meilup :',NBPIVP,' |pivot(s)|/|ligne|<',PETIT
  194.          ENDIF
  195. *!         IF (NBPIVI.GT.0) THEN
  196. *!            WRITE(IOIMP,*) 'WARNING !',
  197. *!     $           ' meilup :',NBPIVI,' pivot(s) < 0'
  198. *!         ENDIF
  199.       ENDIF
  200. *
  201. * Normal termination
  202. *
  203.       IRET=0
  204.       RETURN
  205. *
  206. * Error handling
  207. *
  208.  9999 CONTINUE
  209.       IRET=1
  210.       WRITE(IOIMP,*) 'An error was detected in meilup.eso'
  211.       RETURN
  212. *
  213. *     End of MEILUP
  214. *
  215.       END
  216.  
  217.  
  218.  
  219.  
  220.  
  221.  
  222.  
  223.  

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