Télécharger cgsn.eso

Retour à la liste

Numérotation des lignes :

  1. C CGSN SOURCE PV 16/11/17 21:58:21 9180
  2. SUBROUTINE CGSN(KMORS,KISA,B,X,
  3. $ KPREC,INVDIA,ILUM,ILUI,
  4. $ LRES,LNMV,
  5. $ R,RTLD,P,PHAT,
  6. $ Q,QHAT,UHAT,
  7. $ XP,BP,
  8. $ ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  9. IMPLICIT REAL*8 (A-H,O-Z)
  10. IMPLICIT INTEGER(I-N)
  11. C***********************************************************************
  12. C NOM : CGSN
  13. C DESCRIPTION :
  14. C Résolution d'un système linéaire Ax=b
  15. C par une méthode de gradient conjugué au carré (CGS) avec la
  16. C stratégie de stabilisation de Neumaier préconditionnée ou non.
  17. C
  18. C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  19. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  20. C mél : gounand@semt2.smts.cea.fr
  21. C REFERENCE (bibtex-like) :
  22. C @BOOK{templates,
  23. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  24. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  25. C H. Van der Vorst},
  26. C TITLE={Templates for the Solution of Linear Systems :
  27. C Building Blocks for Iterative Methods},
  28. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  29. C -> URL : http://www.netlib.org/templates/Templates.html
  30. C***********************************************************************
  31. C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
  32. C GNRM2, GDOT (functions)
  33. C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
  34. C PSOLVE : Preconditionner solve
  35. C***********************************************************************
  36. C ENTREES :
  37. C ENTREES/SORTIES :
  38. C SORTIES :
  39. C CODE RETOUR (IRET) : 0 si ok
  40. C <0 si problème (non utilisé !)
  41. C >0 si l'algorithme n'a pas convergé
  42. C***********************************************************************
  43. C VERSION : v1, 10/02/06
  44. C HISTORIQUE : v1, 10/02/06, création
  45. C HISTORIQUE :
  46. C HISTORIQUE :
  47. C***********************************************************************
  48. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  49. C en cas de modification de ce sous-programme afin de faciliter
  50. C la maintenance !
  51. C***********************************************************************
  52. *
  53. * .. Includes .. a commenter si mise au point ok
  54. -INC CCREEL
  55. -INC CCOPTIO
  56. -INC SMLREEL
  57. POINTEUR LRES.MLREEL
  58. -INC SMLENTI
  59. POINTEUR LNMV.MLENTI
  60. POINTEUR KMORS.PMORS
  61. POINTEUR KISA.IZA
  62. POINTEUR INVDIA.IZA
  63. POINTEUR ILUM.PMORS
  64. POINTEUR ILUI.IZA
  65. POINTEUR B.IZA
  66. POINTEUR X.IZA
  67. * .. Work Vectors for CGSN
  68. POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
  69. POINTEUR Q.IZA,QHAT.IZA,UHAT.IZA
  70. POINTEUR U.IZA,VHAT.IZA
  71. * Note QHAT and U, UHAT,VHAT share the same memory space
  72. POINTEUR XP.IZA,BP.IZA
  73. *
  74. INTEGER IMPR
  75. * .. Scalar Arguments ..
  76. INTEGER N,NNZ
  77. INTEGER ITER,IRET
  78. REAL*8 BRTOL,RESID
  79. * ..
  80. *
  81. *
  82. * .. Variables locales
  83. * .. Parameters
  84. REAL*8 ZERO,ONE
  85. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  86. * ..
  87. INTEGER MAXIT,III
  88. REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
  89. REAL*8 RHOTOL
  90. REAL*8 GDOT,GNRM2,RNRM2
  91. REAL*8 MU,MUP
  92. * ..
  93. * .. External subroutines and functions..
  94. * ..
  95. * .. Executable Statements ..
  96. *
  97. *
  98. IRET = 0
  99. INMV = 0
  100. IWARN = 0
  101. MAXIT = ITER
  102. TOL = RESID
  103. *
  104. * Same memory space
  105. U=QHAT
  106. VHAT=UHAT
  107. *
  108. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgsn.eso'
  109. *
  110. * Set breakdown parameter tolerances.
  111. *
  112. RHOTOL = BRTOL
  113. TAUTOL = BRTOL
  114. *
  115. * Set initial residual.
  116. *
  117. C r(0)=b
  118. CALL GCOPY(B,R)
  119. C r(0)=b-Ax
  120. * IF (GNRM2(X).NE.ZERO) THEN
  121. C r(0)=b-Ax
  122. * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
  123. * ENDIF
  124. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  125. INMV=INMV+1
  126. *
  127. ITER = 0
  128. BNRM2 = GNRM2(B)
  129. IF(IMPR.GT.6) THEN
  130. WRITE(IOIMP,*) '||B||=',BNRM2
  131. ENDIF
  132. RNRM2 = GNRM2(R)
  133. IF (ICALRS.EQ.1) BNRM2=RNRM2
  134. *
  135. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  136. RESID = RNRM2 / BNRM2
  137. IF(IMPR.GT.6) THEN
  138. WRITE(IOIMP,*) 'Résidu initial=',RESID
  139. ENDIF
  140. LRES.PROG(ITER+1)=RESID
  141. LNMV.LECT(ITER+1)=INMV
  142. IF (RESID.LE.TOL) GOTO 30
  143. C r[tilde]=r(0)
  144. CALL GCOPY(R,RTLD)
  145. C Neumaier's Strategy
  146. C x=x(0) x'=0
  147. C Normalement, c'est déjà le cas
  148. C b'=r(0) mu' = ||b'||
  149. CALL GCOPY(R,BP)
  150. MUP=RNRM2
  151. *! WRITE(IOIMP,*) 'MUP0=',MUP
  152. *! WRITE(IOIMP,*) '||XP||',GNRM2(XP)
  153. *
  154. 10 CONTINUE
  155. *
  156. * Perform Preconditionned Conjugate Gradient Squared iteration.
  157. *
  158. ITER = ITER + 1
  159. IF(IMPR.GT.7) THEN
  160. WRITE(IOIMP,*) 'ITER=',ITER
  161. ENDIF
  162. *
  163. C rho(i-1)=r[tilde]'r(i-1)
  164. RHO = GDOT(RTLD,R,IDDOT)
  165. * IF ((ITER.LE.10).OR.((ITER/100)*100.EQ.ITER)) THEN
  166. * WRITE(IOIMP,*) 'ITER=',ITER,' RHO=',RHO
  167. * ENDIF
  168. C if rho(i-1)=0 method fails
  169. IF (ABS(RHO).LT.RHOTOL) THEN
  170. RHO=SIGN(XZPREC,RHO)
  171. IWARN=IWARN+1
  172. ENDIF
  173. *
  174. * Compute vector P.
  175. *
  176. IF (ITER.GT.1) THEN
  177. C beta(i-1)=(rho(i-1)/rho(i-2))
  178. BETA = RHO / RHO1
  179. C u(i)=r(i-1)+beta(i-1)*q(i-1)
  180. CALL GCOPY(R,U)
  181. CALL GAXPY(BETA,Q,U)
  182. C p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
  183. CALL GSCAL(BETA**2,P)
  184. CALL GAXPY(BETA,Q,P)
  185. CALL GAXPY(ONE,U,P)
  186. ELSE
  187. C u(1)=r(0)
  188. C p(1)=u(1)
  189. CALL GCOPY(R,U)
  190. CALL GCOPY(U,P)
  191. ENDIF
  192. *
  193. * Compute direction adjusting vector PHAT and scalar ALPHA.
  194. *
  195. C
  196. C Preconditionner solve
  197. C Mp[hat]=p(i)
  198. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  199. $ ,PHAT,P)
  200. C vhat=Ap[hat]
  201. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
  202. INMV=INMV+1
  203. C alpha(i)=rho(i-1)/(r[tilde]'vhat)
  204. TAU = GDOT(RTLD,VHAT,IDDOT)
  205. IF (ABS(TAU).LT.TAUTOL) THEN
  206. TAU=SIGN(XZPREC,TAU)
  207. IWARN=IWARN+1
  208. ENDIF
  209. ALPHA = RHO / TAU
  210. C q(i)=u(i)-alpha(i)vhat
  211. CALL GCOPY(U,Q)
  212. CALL GAXPY(-ALPHA,VHAT,Q)
  213. *
  214. * Compute direction adjusting vector UHAT.
  215. * PHAT is being used as temporary storage here.
  216. *
  217. C M uhat=u(i)+q(i)
  218. CALL GCOPY(Q,PHAT)
  219. CALL GAXPY(ONE,U,PHAT)
  220. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  221. $ ,UHAT,PHAT)
  222. *
  223. * Compute new solution approximation vector X.
  224. *
  225. * x'(i)=x'(i-1)+alpha(i)uhat
  226. CALL GAXPY(ALPHA,UHAT,XP)
  227. *
  228. * Compute residual R and check for tolerance
  229. * using Neumaier's strategy
  230. *
  231. C r(i)=b'-Ax'(i)
  232. CALL GCOPY(BP,R)
  233. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
  234. INMV=INMV+1
  235. MU=GNRM2(R)
  236. RESID = MU / BNRM2
  237. CC qhat= A uhat
  238. C CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
  239. C INMV=INMV+1
  240. CC r(i)=r(i-1)-alpha(i) qhat
  241. C CALL GAXPY(-ALPHA,QHAT,R)
  242. C RESID = GNRM2(R) / BNRM2
  243.  
  244.  
  245. IF(IMPR.GT.6) THEN
  246. WRITE(IOIMP,*) 'Résidu =',RESID
  247. ENDIF
  248. LRES.PROG(ITER+1)=RESID
  249. LNMV.LECT(ITER+1)=INMV
  250. IF (RESID.LE.TOL) THEN
  251. * x(i)=x(i)+x'(i)
  252. CALL GAXPY(ONE,XP,X)
  253. GOTO 30
  254. ENDIF
  255. * Neumaier's strategy
  256. IF (MU.LT.MUP) THEN
  257. * x(i)=x(i)+x'(i)
  258. CALL GAXPY(ONE,XP,X)
  259. * x'(i)=0
  260. CALL GSCAL(ZERO,XP)
  261. * b'=r(i)
  262. CALL GCOPY(R,BP)
  263. * mu'=mu
  264. MUP = MU
  265. * WRITE(IOIMP,*) 'ITER=',ITER
  266. * WRITE(IOIMP,*) 'MUP=',MUP
  267. ENDIF
  268. *
  269. * Iteration fails.
  270. *
  271. * IF (ITER.EQ.MAXIT) THEN
  272. IF (INMV.GE.MAXIT) THEN
  273. IRET = 1
  274. IF (IMPR.GT.0) THEN
  275. C WRITE(IOIMP,*)
  276. C $ 'cgsn.eso : Convergence to tolerance not achieved'
  277. C WRITE(IOIMP,*) 'INMV=',INMV,
  278. C $ ' TOL=',TOL,
  279. C $ ' RESID=',RESID
  280. IF (IWARN.GT.0) THEN
  281. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
  282. $ ,' times.'
  283. ENDIF
  284. ENDIF
  285. RETURN
  286. ENDIF
  287. *
  288. * Do some more
  289. *
  290. RHO1 = RHO
  291. GOTO 10
  292. *
  293. * Iteration successful
  294. *
  295. 30 CONTINUE
  296. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  297. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
  298. ENDIF
  299. IRET = 0
  300. RETURN
  301. *
  302. * End of CGSN.
  303. *
  304. END
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  

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