Télécharger cgsn.eso

Retour à la liste

Numérotation des lignes :

cgsn
  1. C CGSN SOURCE GOUNAND 22/08/25 21:15:03 11434
  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.  
  56. -INC PPARAM
  57. -INC CCOPTIO
  58. -INC SMLREEL
  59. POINTEUR LRES.MLREEL
  60. -INC SMLENTI
  61. POINTEUR LNMV.MLENTI
  62. POINTEUR KMORS.PMORS
  63. POINTEUR KISA.IZA
  64. POINTEUR INVDIA.IZA
  65. POINTEUR ILUM.PMORS
  66. POINTEUR ILUI.IZA
  67. POINTEUR B.IZA
  68. POINTEUR X.IZA
  69. * .. Work Vectors for CGSN
  70. POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
  71. POINTEUR Q.IZA,QHAT.IZA,UHAT.IZA
  72. POINTEUR U.IZA,VHAT.IZA
  73. * Note QHAT and U, UHAT,VHAT share the same memory space
  74. POINTEUR XP.IZA,BP.IZA
  75. *
  76. INTEGER IMPR
  77. * .. Scalar Arguments ..
  78. INTEGER N,NNZ
  79. INTEGER ITER,IRET
  80. REAL*8 BRTOL,RESID
  81. * ..
  82. *
  83. *
  84. * .. Variables locales
  85. * .. Parameters
  86. REAL*8 ZERO,ONE
  87. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  88. * ..
  89. INTEGER MAXIT,III
  90. REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
  91. REAL*8 RHOTOL
  92. REAL*8 GDOT,GNRM2,RNRM2
  93. REAL*8 MU,MUP
  94. * ..
  95. * .. External subroutines and functions..
  96. * ..
  97. * .. Executable Statements ..
  98. *
  99. *
  100. IRET = 0
  101. INMV = 0
  102. IWARN = 0
  103. MAXIT = ITER
  104. TOL = RESID
  105. *
  106. * Same memory space
  107. U=QHAT
  108. VHAT=UHAT
  109. *
  110. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgsn.eso'
  111. *
  112. * Set breakdown parameter tolerances.
  113. *
  114. RHOTOL = BRTOL
  115. TAUTOL = BRTOL
  116. *
  117. * Set initial residual.
  118. *
  119. C r(0)=b
  120. CALL GCOPY(B,R)
  121. C r(0)=b-Ax
  122. * IF (GNRM2(X).NE.ZERO) THEN
  123. C r(0)=b-Ax
  124. * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
  125. * ENDIF
  126. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  127. INMV=INMV+1
  128. *
  129. ITER = 0
  130. BNRM2 = GNRM2(B)
  131. IF(IMPR.GT.6) THEN
  132. WRITE(IOIMP,*) '||B||=',BNRM2
  133. ENDIF
  134. RNRM2 = GNRM2(R)
  135. IF (ICALRS.EQ.1) BNRM2=RNRM2
  136. *
  137. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  138. RESID = RNRM2 / BNRM2
  139. IF(IMPR.GT.6) THEN
  140. WRITE(IOIMP,*) 'Résidu initial=',RESID
  141. ENDIF
  142. LRES.PROG(ITER+1)=RESID
  143. LNMV.LECT(ITER+1)=INMV
  144. IF (RESID.LE.TOL) GOTO 30
  145. C r[tilde]=r(0)
  146. CALL GCOPY(R,RTLD)
  147. C Neumaier's Strategy
  148. C x=x(0) x'=0
  149. C Normalement, c'est déjà le cas
  150. C b'=r(0) mu' = ||b'||
  151. CALL GCOPY(R,BP)
  152. MUP=RNRM2
  153. *! WRITE(IOIMP,*) 'MUP0=',MUP
  154. *! WRITE(IOIMP,*) '||XP||',GNRM2(XP)
  155. *
  156. 10 CONTINUE
  157. * Gestion du CTRL-C
  158. if (ierr.NE.0) return
  159. *
  160. * Perform Preconditionned Conjugate Gradient Squared iteration.
  161. *
  162. ITER = ITER + 1
  163. IF(IMPR.GT.7) THEN
  164. WRITE(IOIMP,*) 'ITER=',ITER
  165. ENDIF
  166. *
  167. C rho(i-1)=r[tilde]'r(i-1)
  168. RHO = GDOT(RTLD,R,IDDOT)
  169. * IF ((ITER.LE.10).OR.((ITER/100)*100.EQ.ITER)) THEN
  170. * WRITE(IOIMP,*) 'ITER=',ITER,' RHO=',RHO
  171. * ENDIF
  172. C if rho(i-1)=0 method fails
  173. IF (ABS(RHO).LT.RHOTOL) THEN
  174. RHO=SIGN(XZPREC,RHO)
  175. IWARN=IWARN+1
  176. ENDIF
  177. *
  178. * Compute vector P.
  179. *
  180. IF (ITER.GT.1) THEN
  181. C beta(i-1)=(rho(i-1)/rho(i-2))
  182. BETA = RHO / RHO1
  183. C u(i)=r(i-1)+beta(i-1)*q(i-1)
  184. CALL GCOPY(R,U)
  185. CALL GAXPY(BETA,Q,U)
  186. C p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
  187. CALL GSCAL(BETA**2,P)
  188. CALL GAXPY(BETA,Q,P)
  189. CALL GAXPY(ONE,U,P)
  190. ELSE
  191. C u(1)=r(0)
  192. C p(1)=u(1)
  193. CALL GCOPY(R,U)
  194. CALL GCOPY(U,P)
  195. ENDIF
  196. *
  197. * Compute direction adjusting vector PHAT and scalar ALPHA.
  198. *
  199. C
  200. C Preconditionner solve
  201. C Mp[hat]=p(i)
  202. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  203. $ ,PHAT,P)
  204. C vhat=Ap[hat]
  205. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
  206. INMV=INMV+1
  207. C alpha(i)=rho(i-1)/(r[tilde]'vhat)
  208. TAU = GDOT(RTLD,VHAT,IDDOT)
  209. IF (ABS(TAU).LT.TAUTOL) THEN
  210. TAU=SIGN(XZPREC,TAU)
  211. IWARN=IWARN+1
  212. ENDIF
  213. ALPHA = RHO / TAU
  214. C q(i)=u(i)-alpha(i)vhat
  215. CALL GCOPY(U,Q)
  216. CALL GAXPY(-ALPHA,VHAT,Q)
  217. *
  218. * Compute direction adjusting vector UHAT.
  219. * PHAT is being used as temporary storage here.
  220. *
  221. C M uhat=u(i)+q(i)
  222. CALL GCOPY(Q,PHAT)
  223. CALL GAXPY(ONE,U,PHAT)
  224. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  225. $ ,UHAT,PHAT)
  226. *
  227. * Compute new solution approximation vector X.
  228. *
  229. * x'(i)=x'(i-1)+alpha(i)uhat
  230. CALL GAXPY(ALPHA,UHAT,XP)
  231. *
  232. * Compute residual R and check for tolerance
  233. * using Neumaier's strategy
  234. *
  235. C r(i)=b'-Ax'(i)
  236. CALL GCOPY(BP,R)
  237. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
  238. INMV=INMV+1
  239. MU=GNRM2(R)
  240. RESID = MU / BNRM2
  241. CC qhat= A uhat
  242. C CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
  243. C INMV=INMV+1
  244. CC r(i)=r(i-1)-alpha(i) qhat
  245. C CALL GAXPY(-ALPHA,QHAT,R)
  246. C RESID = GNRM2(R) / BNRM2
  247.  
  248.  
  249. IF(IMPR.GT.6) THEN
  250. WRITE(IOIMP,*) 'Résidu =',RESID
  251. ENDIF
  252. LRES.PROG(ITER+1)=RESID
  253. LNMV.LECT(ITER+1)=INMV
  254. IF (RESID.LE.TOL) THEN
  255. * x(i)=x(i)+x'(i)
  256. CALL GAXPY(ONE,XP,X)
  257. GOTO 30
  258. ENDIF
  259. * Neumaier's strategy
  260. IF (MU.LT.MUP) THEN
  261. * x(i)=x(i)+x'(i)
  262. CALL GAXPY(ONE,XP,X)
  263. * x'(i)=0
  264. CALL GSCAL(ZERO,XP)
  265. * b'=r(i)
  266. CALL GCOPY(R,BP)
  267. * mu'=mu
  268. MUP = MU
  269. * WRITE(IOIMP,*) 'ITER=',ITER
  270. * WRITE(IOIMP,*) 'MUP=',MUP
  271. ENDIF
  272. *
  273. * Iteration fails.
  274. *
  275. * IF (ITER.EQ.MAXIT) THEN
  276. IF (INMV.GE.MAXIT) THEN
  277. IRET = 1
  278. IF (IMPR.GT.0) THEN
  279. C WRITE(IOIMP,*)
  280. C $ 'cgsn.eso : Convergence to tolerance not achieved'
  281. C WRITE(IOIMP,*) 'INMV=',INMV,
  282. C $ ' TOL=',TOL,
  283. C $ ' RESID=',RESID
  284. IF (IWARN.GT.0) THEN
  285. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
  286. $ ,' times.'
  287. ENDIF
  288. ENDIF
  289. RETURN
  290. ENDIF
  291. *
  292. * Do some more
  293. *
  294. RHO1 = RHO
  295. GOTO 10
  296. *
  297. * Iteration successful
  298. *
  299. 30 CONTINUE
  300. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  301. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
  302. ENDIF
  303. IRET = 0
  304. RETURN
  305. *
  306. * End of CGSN.
  307. *
  308. END
  309.  
  310.  

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