Télécharger prgmr2.eso

Retour à la liste

Numérotation des lignes :

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

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