Télécharger prcg2.eso

Retour à la liste

Numérotation des lignes :

prcg2
  1. C PRCG2 SOURCE GOUNAND 22/08/25 21:15:09 11434
  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.  
  114. -INC PPARAM
  115. -INC CCOPTIO
  116. -INC SMLREEL
  117. INTEGER JG
  118. POINTEUR LRES.MLREEL
  119. -INC SMLENTI
  120. POINTEUR LNMV.MLENTI
  121. POINTEUR MAPREC.MATRIK
  122. POINTEUR KMORS.PMORS
  123. POINTEUR KISA.IZA
  124. POINTEUR KS2B.IZA
  125. POINTEUR INCX.IZA
  126. POINTEUR INVDIA.IZA
  127. POINTEUR ILUM.PMORS
  128. POINTEUR ILUI.IZA
  129. * .. Work Vectors for CG
  130. POINTEUR IR.IZA,IP.IZA,IQ.IZA,IZ.IZA
  131. * .. Scalar Arguments ..
  132. INTEGER ITER, KPREC, IMPR, IRET
  133. REAL*8 RESID
  134. INTEGER NBVA,NJA,NTT,NTTPRE
  135. * .. Executable Statements ..
  136. *
  137. IRET = 0
  138. *
  139. * On récupère les paramètres
  140. *
  141. SEGACT MATRIK
  142. SEGACT MAPREC
  143. IF (KSYM.NE.0) THEN
  144. IF (IMPR.GT.0) THEN
  145. WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
  146. $ 'use Bi-CGSTAB instead !'
  147. * ENDIF
  148. * IRET=-1
  149. * GOTO 9999
  150. ENDIF
  151. ENDIF
  152. * Pour le préconditionneur
  153. ILUM =MAPREC.KIDMAT(6)
  154. ILUI =MAPREC.KIDMAT(7)
  155. IDMAT=MAPREC.KIDMAT(1)
  156. SEGACT IDMAT
  157. INVDIA=IDIAG
  158. SEGDES IDMAT
  159. SEGACT KMORS
  160. NTT =KMORS.IA(/1)-1
  161. * NJA =KMORS.JA(/1)
  162. SEGACT KISA
  163. SEGACT KS2B
  164. SEGACT INCX*MOD
  165. C Paramètres des préconditionnements diagonaux et D-ILU
  166. IF (KPREC.NE.0) THEN
  167. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  168. C Est-il compatible avec la matrice ?
  169. SEGACT INVDIA
  170. NTTPRE=INVDIA.A(/1)
  171. IF (NTTPRE.NE.NTT) THEN
  172. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  173. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  174. WRITE(IOIMP,*) 'NTT=',NTT
  175. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  176. IRET=-2
  177. GOTO 9999
  178. ENDIF
  179. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
  180. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  181. SEGACT ILUM
  182. SEGACT ILUI
  183. NTTPRE=ILUM.IA(/1)
  184. IF (NTTPRE.NE.NTT) THEN
  185. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  186. $ 'ne sont pas compatibles...'
  187. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  188. IRET=-2
  189. GOTO 9999
  190. ENDIF
  191. ELSE
  192. WRITE(IOIMP,*) 'Erreur de programmation'
  193. GOTO 9999
  194. ENDIF
  195. ENDIF
  196. C
  197. C Initialisation de l'historique de convergence
  198. C
  199. JG=ITER+1
  200. SEGINI LNMV
  201. SEGINI LRES
  202. C
  203. C
  204. C Initialisation des vecteurs de travail pour le gradient conjugué
  205. C
  206. NBVA=NTT
  207. SEGINI IR,IP,IQ,IZ
  208. C
  209. C Appel du gradient conjugué
  210. C
  211. CALL CG2(KMORS,KISA,KS2B,INCX,
  212. $ KPREC,INVDIA,ILUM,ILUI,
  213. $ LRES,LNMV,
  214. $ IR,IP,IQ,IZ,
  215. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  216. * Gestion du CTRL-C
  217. if (ierr.NE.0) return
  218. C
  219. C Désactivation
  220. C
  221. SEGSUP IR,IP,IQ,IZ
  222. JG=ITER+1
  223. SEGADJ LRES
  224. SEGDES LRES
  225. SEGADJ LNMV
  226. SEGDES LNMV
  227. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  228. SEGDES INVDIA
  229. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  230. SEGDES ILUI
  231. SEGDES ILUM
  232. ENDIF
  233. SEGDES INCX
  234. SEGDES KS2B
  235. SEGDES KISA
  236. SEGDES KMORS
  237. SEGDES MAPREC
  238. SEGDES MATRIK
  239. *
  240. * Normal termination
  241. *
  242. RETURN
  243. *
  244. * Format handling
  245. *
  246. *
  247. * Error handling
  248. *
  249. 9999 CONTINUE
  250. WRITE(IOIMP,*) 'An error was detected in prcg2.eso'
  251. RETURN
  252. *
  253. * End of PRCG2.
  254. *
  255. END
  256.  
  257.  

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