Télécharger gmr2.eso

Retour à la liste

Numérotation des lignes :

  1. C GMR2 SOURCE PV 16/11/17 21:59:31 9180
  2. SUBROUTINE GMR2(KMORS,KISA,B,X,
  3. $ KPREC,INVDIA,ILUM,ILUI,
  4. $ LRES,LNMV,
  5. $ R,W,AV,
  6. $ Y,GCOS,GSIN,S,
  7. $ V, H,
  8. $ ITER,INMV,RESTRT,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  9. IMPLICIT INTEGER(I-N)
  10. IMPLICIT REAL*8 (A-H,O-Z)
  11. C***********************************************************************
  12. C NOM : GMR2
  13. C DESCRIPTION :
  14. C DESCRIPTION :
  15. C Résolution d'un système linéaire Ax=b
  16. C par une méthode GMRES(m) (restarted Generalized Minimal
  17. C Residual) préconditionné par une factorisation ILU(0)
  18. C (Crout incomplet) ou MILU(0) (Crout inc. modifié relaxé) de A.
  19. C
  20. C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  21. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  22. C mél : gounand@semt2.smts.cea.fr
  23. C REFERENCE (bibtex-like) :
  24. C @BOOK{templates,
  25. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  26. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  27. C H. Van der Vorst},
  28. C TITLE={Templates for the Solution of Linear Systems :
  29. C Building Blocks for Iterative Methods},
  30. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  31. C -> URL : http://www.netlib.org/templates/Templates.html
  32. C***********************************************************************
  33. C APPELES (BLAS 1) : DCOPY, DAXPY (subroutines)
  34. C DNRM2, DDOT, DSCAL (functions)
  35. C APPELES (BLAS 2) : DTRSV (subroutines)
  36. C APPELES (Calcul) : DMOMV : produit Matrice Morse - Vecteur
  37. C PSILU0 : Preconditionner solve
  38. C (Montée-Descente adaptée à ILU(0))
  39. C***********************************************************************
  40. C ENTREES : N, NNZ,
  41. C ROWPTR, COLIND, VAL,
  42. C B, JU, JLU, ALU, RESTRT, IMPR
  43. C ENTREES/SORTIES : X,
  44. C R,W,AV,
  45. C Y,GCOS,GSIN,S,
  46. C V,H,
  47. C ITER, RESID
  48. C SORTIES : LRES, IRET
  49. C CODE RETOUR (IRET) : 0 si ok
  50. C <0 si problème (cf. BRTOL)
  51. C >0 si l'algorithme n'a pas convergé
  52. C N : nombre de degrés de liberté
  53. C NNZ : nombre de valeurs non nulles de la matrice Morse
  54. C ROWPTR, COLIND, VAL : pointeur de ligne, index de colonne
  55. C et valeurs de la matrice Morse
  56. C B : vecteur second membre
  57. C ILUVAL : valeurs de la matrice de la factorisation
  58. C ILU(0) ou MILU(0) de A.
  59. C Cette matrice est stockée en Morse, elle a le
  60. C meme profil que A.
  61. C RESTRT : INTEGER
  62. C paramètre de redémarrage (restart) (i.e. m) pour
  63. C la méthode GMRES(m). (cf. prgmr.eso)
  64. C IMPR : niveau d'impression (0..4)
  65. C X : vecteur des inconnues.
  66. C En entrée, contient l'estimation x(0).
  67. C En sortie, la solution trouvée (si la méthode
  68. C a convergé).
  69. C R,W,AV : vecteurs de travail initialisés et détruits
  70. C Y,GCOS,GSIN,S dans prgmr.
  71. C V,H : tableaux de travail...
  72. C ITER : type INTEGER.
  73. C En entrée, contient le nombre maximum
  74. C d'itérations à effectuer.
  75. C En sortie, contient le nombre d'itérations
  76. C réellement effectuées.
  77. C RESID : type REAL*8.
  78. C En entrée, la valeur maximale du résidu à
  79. C convergence de l'algorithme mesurée comme suit :
  80. C norme[L2](b - A*x) / norme[L2]( b )
  81. C En sortie, la valeur effective de cette mesure.
  82. C LRES : pointeur sur segment MLREEL (include SMLREEL)
  83. C contient l'historique de la valeur de RESID en
  84. C fonction du nombre d'itérations effectuées.
  85. C
  86. C***********************************************************************
  87. C VERSION : 20/12/99
  88. C HISTORIQUE : v1, 16/06/98, création
  89. C HISTORIQUE : 20/12/99
  90. C Interfaçage avec le nouveau psilu0 (factorisations incomplètes
  91. C stockées au format MSR (Modified Sparse Row).
  92. C HISTORIQUE :
  93. C HISTORIQUE :
  94. C***********************************************************************
  95. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  96. C en cas de modification de ce sous-programme afin de faciliter
  97. C la maintenance !
  98. C***********************************************************************
  99. *
  100. * .. Includes .. a commenter si mise au point ok
  101. -INC CCREEL
  102. -INC CCOPTIO
  103. -INC SMLREEL
  104. POINTEUR LRES.MLREEL
  105. -INC SMLENTI
  106. POINTEUR LNMV.MLENTI
  107. POINTEUR KMORS.PMORS
  108. POINTEUR KISA.IZA
  109. POINTEUR INVDIA.IZA
  110. POINTEUR ILUM.PMORS
  111. POINTEUR ILUI.IZA
  112. POINTEUR B.IZA
  113. POINTEUR X.IZA
  114. * .. Parameters
  115. * This is a restart parameter for GMRES
  116. SEGMENT SPACE
  117. REAL*8 WRK(LDW,NLDW)
  118. ENDSEGMENT
  119. SEGMENT SPACE2
  120. POINTEUR WRK2(NIZA).IZA
  121. ENDSEGMENT
  122. * .. Work Vectors for GMRES..
  123. POINTEUR VTMP.IZA
  124. POINTEUR V.SPACE2
  125. POINTEUR H.SPACE
  126. POINTEUR R.IZA,W.IZA,AV.IZA
  127. POINTEUR Y.IZA,GCOS.IZA,GSIN.IZA
  128. POINTEUR S.IZA
  129.  
  130. INTEGER IMPR
  131. * .. Scalar Arguments ..
  132. INTEGER N,NNZ
  133. INTEGER ITER,RESTRT,IRET
  134. REAL*8 RESID
  135. * .. Tableaux de travail à supprimer
  136. C REAL*8 R(N),W(N),AV(N)
  137. C REAL*8 Y(RESTRT),GCOS(RESTRT),GSIN(RESTRT)
  138. C REAL*8 S(RESTRT+1)
  139. C REAL*8 V(N,RESTRT+1)
  140. C REAL*8 H(RESTRT+1,RESTRT)
  141. *
  142. *
  143. * .. Variables locales
  144. * .. Parameters
  145. REAL*8 ZERO,ONE
  146. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  147. * ..
  148. INTEGER I,MAXIT,J,K
  149. REAL*8 BNRM2,TOL,GDOT,GNRM2
  150. REAL*8 H1,H2,RAY
  151. * ..
  152. * .. External subroutines and functions..
  153. * ..
  154. * .. Executable Statements ..
  155. *
  156. *
  157. IRET = 0
  158. INMV = 0
  159. MAXIT = ITER
  160. TOL = RESID
  161. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans gmr2.eso'
  162. *
  163. * Set initial residual.
  164. *
  165. C r=b
  166. CALL GCOPY(B,R)
  167. * IF (GNRM2(X).NE.ZERO) THEN
  168. C r=b-Ax
  169. * CALL GMOMV ('N',ONE,KMORS,KISA,X,ONE,R)
  170. * ENDIF
  171. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  172. INMV=INMV+1
  173. *
  174. ITER = 0
  175. BNRM2 = GNRM2(B)
  176. IF(IMPR.GT.6) THEN
  177. WRITE(IOIMP,*) '||B||=',BNRM2
  178. ENDIF
  179. RNRM2 = GNRM2(R)
  180. IF (ICALRS.EQ.1) BNRM2=RNRM2
  181. *
  182. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  183. RESID = RNRM2 / BNRM2
  184. IF(IMPR.GT.6) THEN
  185. WRITE(IOIMP,*) 'Résidu initial=',RESID
  186. ENDIF
  187. LRES.PROG(ITER+1)=RESID
  188. LNMV.LECT(ITER+1)=INMV
  189. IF (RESID.LE.TOL) GOTO 30
  190. *
  191. 10 CONTINUE
  192. *
  193. * Construct the first column of V, and initialize S to the
  194. * elementary vector E1 scaled by RNORM.
  195. *
  196. C Preconditionner solve (ILU(0)) : Mv(1)=r
  197. VTMP=V.WRK2(1)
  198. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  199. $ ,VTMP,R)
  200.  
  201. C s=||v(1)|| e(1)
  202. C v(1)=v(1) / ||v(1)||
  203. S.A(1) = GNRM2(VTMP)
  204. CALL GSCAL(ONE/S.A(1),VTMP)
  205. *
  206. DO 50 I = 1, RESTRT
  207. ITER = ITER + 1
  208. C av=A v(i)
  209. VTMP=V.WRK2(I)
  210. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,VTMP,ZERO,AV)
  211. INMV=INMV+1
  212. C Preconditionner solve (ILU(0)) : Mw=av
  213. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  214. $ ,W,AV)
  215. *
  216. * Construct I-th column of H so that it is orthonormal to
  217. * the previous I-1 columns.
  218. *
  219. DO 60 K = 1, I
  220. VTMP=V.WRK2(K)
  221. H.WRK(K,I) = GDOT(W,VTMP,IDDOT)
  222. CALL GAXPY(-ONE*H.WRK(K,I),VTMP,W)
  223. 60 CONTINUE
  224. H.WRK(I+1,I)=GNRM2(W)
  225. VTMP=V.WRK2(I+1)
  226. CALL GCOPY(W,VTMP)
  227. CALL GSCAL(ONE/H.WRK(I+1,I),VTMP)
  228. *
  229. C IF (I.GT.1)
  230. *
  231. * Apply Givens rotations to the I-th column of H. This
  232. * effectively reduces the Hessenberg matrix to upper
  233. * triangular form during the RESTRT iterations.
  234. *
  235. DO 70 K = 1, I-1
  236. H1= (GCOS.A(K)*H.WRK(K,I)) + (GSIN.A(K)*H.WRK(K+1,I))
  237. H2=-(GSIN.A(K)*H.WRK(K,I)) + (GCOS.A(K)*H.WRK(K+1,I))
  238. H.WRK(K,I) =H1
  239. H.WRK(K+1,I)=H2
  240. 70 CONTINUE
  241.  
  242. *
  243. * Approximate residual norm. Check tolerance. If okay, compute
  244. * final approximation vector X and quit.
  245. *
  246. RAY=SQRT((H.WRK(I,I)*H.WRK(I,I)) + (H.WRK(I+1,I)*H.WRK(I+1,I)))
  247. GCOS.A(I) = H.WRK(I,I) / RAY
  248. GSIN.A(I) = H.WRK(I+1,I) / RAY
  249. H.WRK(I,I)=RAY
  250. H.WRK(I+1,I)=ZERO
  251. S.A(I+1)=-GSIN.A(I)*S.A(I)
  252. S.A(I) = GCOS.A(I)*S.A(I)
  253. RESID=ABS(S.A(I+1)) / BNRM2
  254. LRES.PROG(ITER+1)=RESID
  255. LNMV.LECT(ITER+1)=INMV
  256. * IF (RESID.LE.TOL.OR.ITER.EQ.MAXIT) THEN
  257. IF (RESID.LE.TOL.OR.INMV.GE.MAXIT) THEN
  258. * Solve the system Hy=S
  259. CALL DCOPY(I,S.A,1,Y.A,1)
  260. CALL DTRSV('UPPER','NOTRANS','NONUNIT',I,H.WRK,RESTRT+1,
  261. $ Y.A,1)
  262. * Updating solution
  263. C CALL UPDATE(I,N,RESTRT,X,WORK2(1,H),WORK(1,Y),
  264. C $ WORK(1,S),WORK(1,V))
  265. DO 503 J=I,1,-1
  266. VTMP=V.WRK2(J)
  267. CALL GAXPY(Y.A(J),VTMP,X)
  268. 503 CONTINUE
  269. IF (RESID.LE.TOL) THEN
  270. GOTO 30
  271. ENDIF
  272. * IF (ITER.EQ.MAXIT) THEN
  273. IF (INMV.GE.MAXIT) THEN
  274. GOTO 40
  275. ENDIF
  276. ENDIF
  277. 50 CONTINUE
  278. *
  279. * Compute current solution vector X.
  280. *
  281. * Solve the system Hy=S
  282. CALL DCOPY(RESTRT,S.A,1,Y.A,1)
  283. CALL DTRSV('UPPER','NOTRANS','NONUNIT',
  284. $ RESTRT,H.WRK,RESTRT+1,Y.A,1)
  285. * Updating solution
  286. DO 303 J=RESTRT,1,-1
  287. VTMP=V.WRK2(J)
  288. CALL GAXPY(Y.A(J),VTMP,X)
  289. 303 CONTINUE
  290. *
  291. * Compute residual vector R, find norm,
  292. * then check for tolerance.
  293. *
  294. CALL GCOPY(B,R)
  295. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  296. INMV=INMV+1
  297. S.A(RESTRT+1) = GNRM2(R)
  298. RESID = S.A(RESTRT+1) / BNRM2
  299. LRES.PROG(ITER+1)=RESID
  300. LNMV.LECT(ITER+1)=INMV
  301. IF (RESID.LE.TOL) GOTO 30
  302. *
  303. * Iteration fails.
  304. *
  305. * IF (ITER.EQ.MAXIT) THEN
  306. IF (INMV.GE.MAXIT) THEN
  307. GOTO 40
  308. ENDIF
  309. *
  310. * Do some more
  311. *
  312. GOTO 10
  313. *
  314. * Iteration successful
  315. *
  316. 30 CONTINUE
  317. IRET = 0
  318. RETURN
  319. *
  320. * Iteration fails
  321. *
  322. 40 CONTINUE
  323. IRET = 1
  324. RETURN
  325. *
  326. * End of GMR2.
  327. *
  328. END
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  

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