Télécharger psdilu.eso

Retour à la liste

Numérotation des lignes :

  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