Télécharger cg2.eso

Retour à la liste

Numérotation des lignes :

  1. C CG2 SOURCE PV 16/11/17 21:58:21 9180
  2. SUBROUTINE CG2(KMORS,KISA,B,X,
  3. $ KPREC,INVDIA,ILUM,ILUI,
  4. $ LRES,LNMV,
  5. $ R,P,Q,Z,
  6. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. C***********************************************************************
  10. C NOM : CG2
  11. C DESCRIPTION :
  12. C Résolution d'un système linéaire Ax=b
  13. C par la méthode du gradient conjugué préconditionné.
  14. C
  15. C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  16. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  17. C mél : gounand@semt2.smts.cea.fr
  18. C REFERENCE (bibtex-like) :
  19. C @BOOK{templates,
  20. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  21. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  22. C H. Van der Vorst},
  23. C TITLE={Templates for the Solution of Linear Systems :
  24. C Building Blocks for Iterative Methods},
  25. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  26. C -> URL : http://www.netlib.org/templates/Templates.html
  27. C***********************************************************************
  28. C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
  29. C GNRM2, GDOT (functions)
  30. C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
  31. C PSOLVE : Preconditionner solve
  32. C***********************************************************************
  33. C ENTREES :
  34. C ENTREES/SORTIES :
  35. C SORTIES :
  36. C CODE RETOUR (IRET) : 0 si ok
  37. C <0 si problème (non utilisé !)
  38. C >0 si l'algorithme n'a pas convergé
  39. C***********************************************************************
  40. C VERSION : v1, 09/02/06
  41. C HISTORIQUE : v1, 09/02/06, création
  42. C HISTORIQUE :
  43. C HISTORIQUE :
  44. C***********************************************************************
  45. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  46. C en cas de modification de ce sous-programme afin de faciliter
  47. C la maintenance !
  48. C***********************************************************************
  49. *
  50. * .. Includes .. a commenter si mise au point ok
  51. -INC CCOPTIO
  52. -INC CCREEL
  53. -INC SMLREEL
  54. POINTEUR LRES.MLREEL
  55. -INC SMLENTI
  56. POINTEUR LNMV.MLENTI
  57. POINTEUR KMORS.PMORS
  58. POINTEUR KISA.IZA
  59. POINTEUR INVDIA.IZA
  60. POINTEUR ILUM.PMORS
  61. POINTEUR ILUI.IZA
  62. POINTEUR B.IZA
  63. POINTEUR X.IZA
  64. * .. Work Vectors for CG
  65. POINTEUR R.IZA,P.IZA,Q.IZA,Z.IZA
  66. * .. Scalar Arguments ..
  67. INTEGER ITER, KPREC, IMPR, IRET
  68. REAL*8 RESID
  69. *
  70. *
  71. * .. Variables locales
  72. * .. Parameters
  73. REAL*8 ZERO , ONE
  74. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  75. * ..
  76. INTEGER MAXIT
  77. REAL*8 ALPHA,BETA,RHO,RHO1,TOL,GNRM2,OMEGA
  78. REAL*8 GDOT,BNRM2
  79. * ..
  80. * .. External subroutines and functions..
  81. EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
  82. * ..
  83. * .. Executable Statements ..
  84. *
  85.  
  86. IRET = 0
  87. INMV = 0
  88. IWARN = 0
  89. IGRAN = 0
  90. MGRAN = 5
  91. MAXIT = ITER
  92. TOL = RESID
  93. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cg2.eso'
  94. *
  95. * Set parameter tolerances.
  96. *
  97. RHOTOL = BRTOL
  98. OMETOL = BRTOL
  99. * RRMAX = 10.D0
  100. RRMAX = 1.D0/(XZPREC**0.25D0)
  101. *
  102. * Calcul du résidu initial
  103. *
  104. C r(0)=b
  105. CALL GCOPY(B,R)
  106.  
  107. C r(0)=b-Ax
  108. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  109. INMV=INMV+1
  110. *
  111. ITER = 0
  112. BNRM2 = GNRM2(B)
  113. IF(IMPR.GT.5) THEN
  114. WRITE(IOIMP,*) '||B||=',BNRM2
  115. ENDIF
  116. RNRM2 = GNRM2(R)
  117. IF (ICALRS.EQ.1) BNRM2=RNRM2
  118. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  119. RESID = RNRM2 / BNRM2
  120. IF(IMPR.GT.5) THEN
  121. WRITE(IOIMP,*) 'Initial residual=',RESID
  122. ENDIF
  123. RESMAX = RESID * RRMAX
  124. LRES.PROG(ITER+1)=RESID
  125. LNMV.LECT(ITER+1)=INMV
  126. *
  127. IF (RESID.LE.TOL) THEN
  128. IRET = 0
  129. GOTO 30
  130. ENDIF
  131.  
  132. *
  133. 10 CONTINUE
  134. *
  135. * Perform (Preconditioned) Conjugate Gradient Iteration.
  136. *
  137. ITER = ITER +1
  138. IF(IMPR.GT.7) THEN
  139. WRITE(IOIMP,*) 'ITER=',ITER
  140. ENDIF
  141. *
  142. * Preconditioner Solve
  143. C Mz(i-1)=r(i-1)
  144. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  145. $ ,Z,R)
  146. C rho(i-1)=r(i-1)'z(i-1)
  147. RHO = GDOT(R,Z,IDDOT)
  148. * WRITE(IOIMP,*) ' RHO=',RHO
  149. IF (ABS(RHO).LT.RHOTOL) THEN
  150. RHO=SIGN(XZPREC,RHO)
  151. IWARN=IWARN+1
  152. ENDIF
  153. *
  154. * Compute direction vector P
  155. *
  156. IF (ITER.GT.1) THEN
  157. C beta(i-1)=rho(i-1)/rho(i-2)
  158. BETA = RHO / RHO1
  159. C p(i)=z(i-1)+beta(i-1)p(i-1)
  160. CALL GAXPY(BETA,P,Z)
  161. CALL GCOPY(Z,P)
  162. ELSE
  163. C p(1)=z(0)
  164. CALL GCOPY(Z,P)
  165. ENDIF
  166. *
  167. * Compute scalar ALPHA (save A*p to q)
  168. *
  169. C q(i)=Ap(i)
  170. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,P,ZERO,Q)
  171. INMV=INMV+1
  172. C alpha(i)=rho(i-1)/(p(i)'q(i))
  173. OMEGA = GDOT(P,Q,IDDOT)
  174. * WRITE(IOIMP,*) ' OMEGA=',OMEGA
  175. IF (ABS(OMEGA).LT.OMETOL) THEN
  176. OMEGA=SIGN(XZPREC,OMEGA)
  177. IWARN=IWARN+1
  178. ENDIF
  179. ALPHA = RHO / OMEGA
  180. *
  181. * Compute current solution vector x
  182. *
  183. C x(i)=x(i-1)+alpha(i)p(i)
  184. CALL GAXPY(ALPHA,P,X)
  185. *
  186. * Compute residual vector R, find norm,
  187. * then check for tolerance
  188. *
  189. C r(i)=r(i-1)-alpha(i)r(i)
  190. CALL GAXPY(-ALPHA,Q,R)
  191. *
  192. RESID = GNRM2(R) / BNRM2
  193. IF(IMPR.GT.7) THEN
  194. WRITE(IOIMP,*) 'Résidu=',RESID
  195. ENDIF
  196. LRES.PROG(ITER+1)=RESID
  197. LNMV.LECT(ITER+1)=INMV
  198. IF (RESID.GT.RESMAX) THEN
  199. IGRAN=IGRAN+1
  200. ELSE
  201. IGRAN=0
  202. ENDIF
  203. *
  204. * Iteration successful
  205. *
  206. IF (RESID.LE.TOL) THEN
  207. IRET = 0
  208. *
  209. * Iteration fails
  210. *
  211. * ELSE IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  212. ELSE IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  213. IRET = 1
  214. *
  215. * Do some more
  216. *
  217. ELSE
  218. RHO1 = RHO
  219. GOTO 10
  220. ENDIF
  221. *
  222. * Normal termination
  223. *
  224. 30 CONTINUE
  225. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  226. WRITE(IOIMP,*) 'CG nearly failed ',IWARN
  227. $ ,' times.'
  228. ENDIF
  229. RETURN
  230. *
  231. * End of CG2
  232. *
  233. END
  234.  
  235.  
  236.  
  237.  
  238.  
  239.  
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  

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