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.  
  131. -INC PPARAM
  132. -INC CCOPTIO
  133. -INC SMLREEL
  134. INTEGER JG
  135. POINTEUR LRES.MLREEL
  136. -INC SMLENTI
  137. POINTEUR LNMV.MLENTI
  138. POINTEUR MAPREC.MATRIK
  139. POINTEUR KMORS.PMORS
  140. POINTEUR KISA.IZA
  141. POINTEUR KS2B.IZA
  142. POINTEUR INCX.IZA
  143. POINTEUR INVDIA.IZA
  144. POINTEUR ILUM.PMORS
  145. POINTEUR ILUI.IZA
  146. * .. Parameters
  147. * This is a restart parameter for GMRES
  148. INTEGER RESTRT
  149. INTEGER LDW,NLDW
  150. SEGMENT SPACE
  151. REAL*8 WRK(LDW,NLDW)
  152. ENDSEGMENT
  153. SEGMENT SPACE2
  154. POINTEUR WRK2(NIZA).IZA
  155. ENDSEGMENT
  156. * .. Work Vectors for GMRES..
  157. POINTEUR ITMP.IZA
  158. POINTEUR IV.SPACE2
  159. POINTEUR IH.SPACE
  160. POINTEUR IR.IZA,IW.IZA,IAV.IZA
  161. POINTEUR IY.IZA,IGCOS.IZA,IGSIN.IZA
  162. POINTEUR IS.IZA
  163. * .. Scalar Arguments ..
  164. INTEGER ITER, KPREC, IMPR, IRET
  165. REAL*8 RESID
  166. INTEGER NBVA,NJA,NTT,NTTPRE
  167. * .. Executable Statements ..
  168. *
  169. IRET = 0
  170. *
  171. * On récupère les paramètres
  172. *
  173. SEGACT MATRIK
  174. SEGACT MAPREC
  175. IF (KSYM.EQ.0) THEN
  176. IF (IMPR.GT.0) THEN
  177. WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
  178. $ 'use a Conjugate Gradient instead !'
  179. ENDIF
  180. C IRET=-2
  181. C GOTO 9999
  182. ENDIF
  183. * Pour le préconditionneur
  184. ILUM =MAPREC.KIDMAT(6)
  185. ILUI =MAPREC.KIDMAT(7)
  186. IDMAT=MAPREC.KIDMAT(1)
  187. SEGACT IDMAT
  188. INVDIA=IDIAG
  189. SEGDES IDMAT
  190.  
  191. SEGACT KMORS
  192. NTT =KMORS.IA(/1)-1
  193. * NJA =KMORS.JA(/1)
  194. SEGACT KISA
  195. SEGACT KS2B
  196. SEGACT INCX*MOD
  197. C Paramètres des préconditionnements diagonaux et D-ILU
  198. IF (KPREC.NE.0) THEN
  199. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  200. C Est-il compatible avec la matrice ?
  201. SEGACT INVDIA
  202. NTTPRE=INVDIA.A(/1)
  203. IF (NTTPRE.NE.NTT) THEN
  204. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  205. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  206. WRITE(IOIMP,*) 'NTT=',NTT
  207. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  208. IRET=-2
  209. GOTO 9999
  210. ENDIF
  211. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
  212. C ilutp, ilutpg, ilutpg2
  213. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  214. SEGACT ILUM
  215. SEGACT ILUI
  216. NTTPRE=ILUM.IA(/1)
  217. IF (NTTPRE.NE.NTT) THEN
  218. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  219. $ 'ne sont pas compatibles...'
  220. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  221. IRET=-2
  222. GOTO 9999
  223. ENDIF
  224. ENDIF
  225. ENDIF
  226. C
  227. C Initialisation de l'historique de convergence
  228. C
  229. JG=ITER+1
  230. SEGINI LNMV
  231. SEGINI LRES
  232. C
  233. C
  234. C Initialisation des vecteurs de travail pour GMRES
  235. C
  236. NBVA=NTT
  237. SEGINI IR,IW,IAV
  238. NBVA=RESTRT
  239. SEGINI IY,IGCOS,IGSIN
  240. NBVA=RESTRT+1
  241. SEGINI IS
  242. NBVA=NTT
  243. NIZA=RESTRT+1
  244. SEGINI IV
  245. DO I=1,NIZA
  246. SEGINI ITMP
  247. IV.WRK2(I)=ITMP
  248. ENDDO
  249. LDW=RESTRT+1
  250. NLDW=RESTRT
  251. SEGINI IH
  252. C
  253. CALL GMR2(KMORS,KISA,KS2B,INCX,
  254. $ KPREC,INVDIA,ILUM,ILUI,
  255. $ LRES,LNMV,
  256. $ IR,IW,IAV,
  257. $ IY,IGCOS,IGSIN,IS,
  258. $ IV, IH,
  259. $ ITER,INMV,RESTRT,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  260. C
  261. C Désactivation-suppression
  262. C
  263. SEGSUP IH
  264. DO I=1,NIZA
  265. ITMP=IV.WRK2(I)
  266. SEGSUP ITMP
  267. ENDDO
  268. SEGSUP IV
  269. SEGSUP IS
  270. SEGSUP IY,IGCOS,IGSIN
  271. SEGSUP IR,IW,IAV
  272. JG=ITER+1
  273. SEGADJ LRES
  274. SEGDES LRES
  275. SEGADJ LNMV
  276. SEGDES LNMV
  277. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  278. SEGDES INVDIA
  279. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  280. SEGDES ILUI
  281. SEGDES ILUM
  282. ENDIF
  283. SEGDES INCX
  284. SEGDES KS2B
  285. SEGDES KISA
  286. SEGDES KMORS
  287. SEGDES MAPREC
  288. SEGDES MATRIK
  289. C
  290. C A problem has occured in the GMRES method
  291. C
  292. IF (IRET.LT.0) GOTO 9999
  293. *
  294. * Normal termination
  295. *
  296. RETURN
  297. *
  298. * Format handling
  299. *
  300. *
  301. * Error handling
  302. *
  303. 9999 CONTINUE
  304. WRITE(IOIMP,*) 'An error was detected in prgmr2.eso'
  305. RETURN
  306. *
  307. * End of prgmr2
  308. *
  309. END
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  

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