Télécharger prcg2.eso

Retour à la liste

Numérotation des lignes :

  1. C PRCG2 SOURCE PV 16/11/17 22:01:03 9180
  2. SUBROUTINE PRCG2(KMORS,KISA,KS2B,MATRIK,MAPREC,
  3. $ LRES,LNMV,INCX,ITER,INMV,RESID,KPREC,BRTOL,ICALRS,IDDOT,
  4. $ IMVEC,
  5. $ IMPR,IRET)
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. C***********************************************************************
  9. C NOM : PRCG2
  10. C DESCRIPTION :
  11. C Préparation de la résolution d'un système linéaire Ax=b
  12. C par la méthode du gradient conjugué (préconditionnée ou non)
  13. C
  14. C A doit etre symétrique.
  15. C
  16. C Les 4 préconditionnements disponibles pour cette méthode
  17. C sont :
  18. C * Jacobi (diagonal)
  19. C * D-ILU (Diagonal Incomplete LU factorization)
  20. C * ILU(0) (Incomplete LU factorization
  21. C of level zero ie Choleski incomplet)
  22. C * MILU(0) (Modified ILU(0)) relaxé
  23. C 5 : préconditionnement ILUT (dual truncation)
  24. C 6 : préconditionnement ILUT2 (une variante du
  25. C précédent qui remplit mieux la mémoire et
  26. C fonctionne mieux quelquefois)
  27. C
  28. C Ces préconditionneurs sont supposés déjà construits (ie
  29. C factorisés), voir les subroutines MEIDIA, PRDILU, PRILU0,
  30. C PRMILU et PRILUT.
  31. C
  32. C Ce sous-programme est en fait une interface aux sous-programmes :
  33. C cg, cgd, cgdilu, cgilu0
  34. C qui sont en Fortran presque pur (pour raison de rapidité)
  35. C et effectuent la résolution du système linéaire
  36. C proprement dite.
  37. C
  38. C LANGAGE : ESOPE
  39. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  40. C mél : gounand@semt2.smts.cea.fr
  41. C REFERENCE (bibtex-like) :
  42. C @BOOK{templates,
  43. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  44. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  45. C H. Van der Vorst},
  46. C TITLE={Templates for the Solution of Linear Systems :
  47. C Building Blocks for Iterative Methods},
  48. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  49. C -> URL : http://www.netlib.org/templates/Templates.html
  50. C***********************************************************************
  51. C APPELES (Calcul) : CG, CGD, CGDILU, CGILU0
  52. C***********************************************************************
  53. C ENTREES : MATRIK, MAPREC, KPREC, IMPR
  54. C ENTREES/SORTIES : INCX, ITER, RESID
  55. C SORTIES : LRES, IRET
  56. C CODE RETOUR (IRET) : 0 si ok
  57. C <0 si problème
  58. C >0 si l'algorithme appelé n'a pas convergé
  59. C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
  60. C on pioche dedans les informations nécessaires
  61. C (différents pointeurs, nb. de ddl...)
  62. C MAPREC : pointeur sur segment MATRIK de l'include SMMATRIK
  63. C on va l'utiliser comme préconditionneur
  64. C KPREC : type de préconditionnement à effectuer
  65. C 0 : pas de préconditionnement
  66. C 1 : préconditionnement Jacobi (diagonal)
  67. C 2 : préconditionnement D-ILU
  68. C 3 : préconditionnement ILU(0) (Choleski incomplet)
  69. C 4 : préconditionnement MILU(0)
  70. C (Choleski incomplet modifié)
  71. C 5 : préconditionnement ILUT (dual truncation)
  72. C 6 : préconditionnement ILUT2 (une variante du
  73. C précédent qui remplit mieux la mémoire et
  74. C fonctionne mieux quelquefois)
  75. C IMPR : niveau d'impression (0..4)
  76. C INCX : pointeur sur segment IZA de l'include SMMATRIK
  77. C C'est le vecteur des inconnues.
  78. C En entrée, contient l'estimation x(0).
  79. C En sortie, la solution trouvée (si la méthode
  80. C a convergé).
  81. C ITER : type INTEGER.
  82. C En entrée, contient le nombre maximum
  83. C d'itérations à effectuer.
  84. C En sortie, contient le nombre d'itérations
  85. C réellement effectuées.
  86. C RESID : type REAL*8.
  87. C En entrée, la valeur maximale du résidu à
  88. C convergence de l'algorithme mesurée comme suit :
  89. C norme[L2](b - A*x) / norme[L2]( b )
  90. C En sortie, la valeur effective de cette mesure.
  91. C LRES : pointeur sur segment MLREEL (include SMLREEL)
  92. C contient l'historique de la valeur de RESID en
  93. C fonction du nombre d'itérations effectuées.
  94. C
  95. C***********************************************************************
  96. C VERSION : v1, 02/04/98, version initiale
  97. C HISTORIQUE : v1, 02/04/98, création
  98. C HISTORIQUE : 09/02/99 : ajout préconditionneur <> MATRIK
  99. C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau cgilu0
  100. C (factorisations incomplètes stockées au format MSR (Modified Sparse
  101. C Row) (voir la doc de Sparskit version 2+)
  102. C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  103. C HISTORIQUE : 22/03/00 : ajout des préconditionneurs ILUT
  104. C HISTORIQUE : 10/02/06 : Nettoyage
  105. C HISTORIQUE :
  106. C***********************************************************************
  107. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  108. C en cas de modification de ce sous-programme afin de faciliter
  109. C la maintenance !
  110. C***********************************************************************
  111. *
  112. * .. Includes et pointeurs associés ..
  113. -INC CCOPTIO
  114. -INC SMLREEL
  115. INTEGER JG
  116. POINTEUR LRES.MLREEL
  117. -INC SMLENTI
  118. POINTEUR LNMV.MLENTI
  119. POINTEUR MAPREC.MATRIK
  120. POINTEUR KMORS.PMORS
  121. POINTEUR KISA.IZA
  122. POINTEUR KS2B.IZA
  123. POINTEUR INCX.IZA
  124. POINTEUR INVDIA.IZA
  125. POINTEUR ILUM.PMORS
  126. POINTEUR ILUI.IZA
  127. * .. Work Vectors for CG
  128. POINTEUR IR.IZA,IP.IZA,IQ.IZA,IZ.IZA
  129. * .. Scalar Arguments ..
  130. INTEGER ITER, KPREC, IMPR, IRET
  131. REAL*8 RESID
  132. INTEGER NBVA,NJA,NTT,NTTPRE
  133. * .. Executable Statements ..
  134. *
  135. IRET = 0
  136. *
  137. * On récupère les paramètres
  138. *
  139. SEGACT MATRIK
  140. SEGACT MAPREC
  141. IF (KSYM.NE.0) THEN
  142. IF (IMPR.GT.0) THEN
  143. WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
  144. $ 'use Bi-CGSTAB instead !'
  145. * ENDIF
  146. * IRET=-1
  147. * GOTO 9999
  148. ENDIF
  149. ENDIF
  150. * Pour le préconditionneur
  151. ILUM =MAPREC.KIDMAT(6)
  152. ILUI =MAPREC.KIDMAT(7)
  153. IDMAT=MAPREC.KIDMAT(1)
  154. SEGACT IDMAT
  155. INVDIA=IDIAG
  156. SEGDES IDMAT
  157. SEGACT KMORS
  158. NTT =KMORS.IA(/1)-1
  159. * NJA =KMORS.JA(/1)
  160. SEGACT KISA
  161. SEGACT KS2B
  162. SEGACT INCX*MOD
  163. C Paramètres des préconditionnements diagonaux et D-ILU
  164. IF (KPREC.NE.0) THEN
  165. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  166. C Est-il compatible avec la matrice ?
  167. SEGACT INVDIA
  168. NTTPRE=INVDIA.A(/1)
  169. IF (NTTPRE.NE.NTT) THEN
  170. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  171. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  172. WRITE(IOIMP,*) 'NTT=',NTT
  173. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  174. IRET=-2
  175. GOTO 9999
  176. ENDIF
  177. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
  178. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  179. SEGACT ILUM
  180. SEGACT ILUI
  181. NTTPRE=ILUM.IA(/1)
  182. IF (NTTPRE.NE.NTT) THEN
  183. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  184. $ 'ne sont pas compatibles...'
  185. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  186. IRET=-2
  187. GOTO 9999
  188. ENDIF
  189. ELSE
  190. WRITE(IOIMP,*) 'Erreur de programmation'
  191. GOTO 9999
  192. ENDIF
  193. ENDIF
  194. C
  195. C Initialisation de l'historique de convergence
  196. C
  197. JG=ITER+1
  198. SEGINI LNMV
  199. SEGINI LRES
  200. C
  201. C
  202. C Initialisation des vecteurs de travail pour le gradient conjugué
  203. C
  204. NBVA=NTT
  205. SEGINI IR,IP,IQ,IZ
  206. C
  207. C Appel du gradient conjugué
  208. C
  209. CALL CG2(KMORS,KISA,KS2B,INCX,
  210. $ KPREC,INVDIA,ILUM,ILUI,
  211. $ LRES,LNMV,
  212. $ IR,IP,IQ,IZ,
  213. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  214. C
  215. C Désactivation
  216. C
  217. SEGSUP IR,IP,IQ,IZ
  218. JG=ITER+1
  219. SEGADJ LRES
  220. SEGDES LRES
  221. SEGADJ LNMV
  222. SEGDES LNMV
  223. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  224. SEGDES INVDIA
  225. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  226. SEGDES ILUI
  227. SEGDES ILUM
  228. ENDIF
  229. SEGDES INCX
  230. SEGDES KS2B
  231. SEGDES KISA
  232. SEGDES KMORS
  233. SEGDES MAPREC
  234. SEGDES MATRIK
  235. *
  236. * Normal termination
  237. *
  238. RETURN
  239. *
  240. * Format handling
  241. *
  242. *
  243. * Error handling
  244. *
  245. 9999 CONTINUE
  246. WRITE(IOIMP,*) 'An error was detected in prcg2.eso'
  247. RETURN
  248. *
  249. * End of PRCG2.
  250. *
  251. END
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  

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