Télécharger la plaquette Cast3M


fr en

Télécharger psdilu.eso

Retour à la liste

Numérotation des lignes : 

psdilu

  1. C PSDILU    SOURCE    CHAT      05/01/13    02:37:23     5004
  2.       SUBROUTINE PSDILU(N,NNZ,ROWPTR,COLIND,VAL,
  3.      $     DINV,
  4.      $     Z,R)
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7. C***********************************************************************
  8. C NOM         : PSDILU
  9. C DESCRIPTION :
  10. C    PSDILU solves the linear system Mz = r
  11. C    M is a preconditionning matrix written as :
  12. C                 M=(D+LA)D-1(D+UA)
  13. C    It is the D-ILU incomplete factorization of matrix A
  14. C    (computed in medilu.eso)
  15. C    LA and UA correspond to the strictly lower and strictly
  16. C    upper part of the matrix A.
  17. C    Ax=b is a system we wish to solve by an iterative method.
  18. C
  19. C    A is stored in Compressed Row Storage (CRS, i.e. Morse)
  20. C    In DINV, the INVERSE of the diagonal elements D
  21. C    of the D-ILU factorisation are stored.
  22. C    They are calculated in subroutine mdilus (A symmetric)
  23. C                        or subroutine mdilun (A non-symm.)
  24. C
  25. C    No tests are made
  26. C
  27. C    Currently, this subroutine is called by : cgdilu, bcgsdi
  28. C
  29. C
  30. C LANGAGE     : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  31. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  32. C               mél : gounand@semt2.smts.cea.fr
  33. C REFERENCE (bibtex-like) :
  34. C  @BOOK{templates,
  35. C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  36. C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  37. C                H. Van der Vorst},
  38. C        TITLE={Templates for the Solution of Linear Systems :
  39. C               Building Blocks for Iterative Methods},
  40. C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  41. C  -> URL : http://www.netlib.org/templates/Templates.html
  42. C***********************************************************************
  43. C APPELES : -
  44. C***********************************************************************
  45. C ENTREES            : N, NNZ, ROWPTR, COLIND, VAL,
  46. C                      DINV, R
  47. C ENTREES/SORTIES    : Z
  48. C SORTIES            : -
  49. C CODE RETOUR (IRET) : -
  50. C     N      : nombre de degrés de liberté
  51. C     NNZ    : nombre de valeurs non nulles de la matrice Morse
  52. C     ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
  53. C                           et valeurs de la matrice Morse A
  54. C     DINV   : vecteur contenant les inverses des pivots
  55. C              de la factorisation D-ILU de A
  56. C     R      : vecteur second membre du système à résoudre.
  57. C     Z      : vecteur des inconnues du système à résoudre.
  58. C              initialisé dans la subroutine appelante.
  59. C***********************************************************************
  60. C VERSION    : v1, 01/04/98, version initiale
  61. C HISTORIQUE : v1, 01/04/98, création
  62. C HISTORIQUE :
  63. C HISTORIQUE :
  64. C***********************************************************************
  65. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  66. C en cas de modification de ce sous-programme afin de faciliter
  67. C la maintenance !
  68. C***********************************************************************
  69. *
  70. *     .. Scalar Arguments ..
  71.       INTEGER  N,NNZ
  72. *     ..
  73. *     .. Array Arguments ..
  74. *     .. Matrices stockées en Morse
  75.       INTEGER ROWPTR(N+1)
  76.       INTEGER COLIND(NNZ)
  77.       REAL*8  VAL(NNZ)
  78. *     .. Vecteurs
  79.       REAL*8  DINV(N)
  80.       REAL*8  Z(N),R(N)
  81. *
  82. *
  83. *     .. Variables locales
  84. *     .. Parameters
  85.       REAL*8     ZERO
  86.       PARAMETER (ZERO=0.0D0)
  87. *     ..
  88.       INTEGER IN,JN
  89.       REAL*8 TMP
  90. *     .. Executable Statements ..
  91. *
  92. *
  93. C     On effectue une montée-descente classique
  94. C
  95. C     Montée :
  96. C     Zi=(1/Dii)(Ri- Somme(j<i) [LijZj])
  97. C
  98. C Une optimisation effectuée : sortie de la boucle sitot
  99. C qu'on atteint un indice j supérieur ou égal à i
  100. C (on suppose qu'ils sont ordonnés dans l'ordre croissant
  101. C dans COLIND)
  102. C
  103.       Z(1)=DINV(1)*R(1)
  104.       DO 100 IN=2,N
  105.          TMP=ZERO
  106.          IR=ROWPTR(IN)
  107.          IRFIN=ROWPTR(IN+1)-1
  108.  101     CONTINUE
  109.          JN=COLIND(IR)
  110.          IF (JN.LT.IN) THEN
  111.             TMP=TMP+(VAL(IR)*Z(JN))
  112.             IR=IR+1
  113.             IF (IR.LE.IRFIN) GOTO 101
  114.          ENDIF
  115.          Z(IN)=DINV(IN)*(R(IN)-TMP)
  116.  100  CONTINUE
  117. C
  118. C     Descente :
  119. C     Zi=Zi-(1/Dii)(Somme(j>i) [UijZj])
  120. C
  121. C Meme optimisation que ci-dessus
  122. C
  123.       DO 200 IN=N-1,1,-1
  124.          TMP=ZERO
  125.          IR=ROWPTR(IN+1)-1
  126.          IRFIN=ROWPTR(IN)
  127.  201     CONTINUE
  128.          JN=COLIND(IR)
  129.          IF (JN.GT.IN) THEN
  130.             TMP=TMP+(VAL(IR)*Z(JN))
  131.             IR=IR-1
  132.             IF (IR.GE.IRFIN) GOTO 201
  133.          ENDIF
  134.          Z(IN)=Z(IN)-DINV(IN)*TMP
  135.  200  CONTINUE
  136. *
  137. *     Normal termination
  138. *
  139.       RETURN
  140. *
  141. *     End of PSDILU.
  142. *
  143.       END
  144.  
  145.  
  146.  
  147.  

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