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.  
  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. *
  158. * Perform Preconditionned Conjugate Gradient Squared iteration.
  159. *
  160. ITER = ITER + 1
  161. IF(IMPR.GT.7) THEN
  162. WRITE(IOIMP,*) 'ITER=',ITER
  163. ENDIF
  164. *
  165. C rho(i-1)=r[tilde]'r(i-1)
  166. RHO = GDOT(RTLD,R,IDDOT)
  167. * IF ((ITER.LE.10).OR.((ITER/100)*100.EQ.ITER)) THEN
  168. * WRITE(IOIMP,*) 'ITER=',ITER,' RHO=',RHO
  169. * ENDIF
  170. C if rho(i-1)=0 method fails
  171. IF (ABS(RHO).LT.RHOTOL) THEN
  172. RHO=SIGN(XZPREC,RHO)
  173. IWARN=IWARN+1
  174. ENDIF
  175. *
  176. * Compute vector P.
  177. *
  178. IF (ITER.GT.1) THEN
  179. C beta(i-1)=(rho(i-1)/rho(i-2))
  180. BETA = RHO / RHO1
  181. C u(i)=r(i-1)+beta(i-1)*q(i-1)
  182. CALL GCOPY(R,U)
  183. CALL GAXPY(BETA,Q,U)
  184. C p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
  185. CALL GSCAL(BETA**2,P)
  186. CALL GAXPY(BETA,Q,P)
  187. CALL GAXPY(ONE,U,P)
  188. ELSE
  189. C u(1)=r(0)
  190. C p(1)=u(1)
  191. CALL GCOPY(R,U)
  192. CALL GCOPY(U,P)
  193. ENDIF
  194. *
  195. * Compute direction adjusting vector PHAT and scalar ALPHA.
  196. *
  197. C
  198. C Preconditionner solve
  199. C Mp[hat]=p(i)
  200. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  201. $ ,PHAT,P)
  202. C vhat=Ap[hat]
  203. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
  204. INMV=INMV+1
  205. C alpha(i)=rho(i-1)/(r[tilde]'vhat)
  206. TAU = GDOT(RTLD,VHAT,IDDOT)
  207. IF (ABS(TAU).LT.TAUTOL) THEN
  208. TAU=SIGN(XZPREC,TAU)
  209. IWARN=IWARN+1
  210. ENDIF
  211. ALPHA = RHO / TAU
  212. C q(i)=u(i)-alpha(i)vhat
  213. CALL GCOPY(U,Q)
  214. CALL GAXPY(-ALPHA,VHAT,Q)
  215. *
  216. * Compute direction adjusting vector UHAT.
  217. * PHAT is being used as temporary storage here.
  218. *
  219. C M uhat=u(i)+q(i)
  220. CALL GCOPY(Q,PHAT)
  221. CALL GAXPY(ONE,U,PHAT)
  222. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  223. $ ,UHAT,PHAT)
  224. *
  225. * Compute new solution approximation vector X.
  226. *
  227. * x'(i)=x'(i-1)+alpha(i)uhat
  228. CALL GAXPY(ALPHA,UHAT,XP)
  229. *
  230. * Compute residual R and check for tolerance
  231. * using Neumaier's strategy
  232. *
  233. C r(i)=b'-Ax'(i)
  234. CALL GCOPY(BP,R)
  235. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
  236. INMV=INMV+1
  237. MU=GNRM2(R)
  238. RESID = MU / BNRM2
  239. CC qhat= A uhat
  240. C CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
  241. C INMV=INMV+1
  242. CC r(i)=r(i-1)-alpha(i) qhat
  243. C CALL GAXPY(-ALPHA,QHAT,R)
  244. C RESID = GNRM2(R) / BNRM2
  245.  
  246.  
  247. IF(IMPR.GT.6) THEN
  248. WRITE(IOIMP,*) 'Résidu =',RESID
  249. ENDIF
  250. LRES.PROG(ITER+1)=RESID
  251. LNMV.LECT(ITER+1)=INMV
  252. IF (RESID.LE.TOL) THEN
  253. * x(i)=x(i)+x'(i)
  254. CALL GAXPY(ONE,XP,X)
  255. GOTO 30
  256. ENDIF
  257. * Neumaier's strategy
  258. IF (MU.LT.MUP) THEN
  259. * x(i)=x(i)+x'(i)
  260. CALL GAXPY(ONE,XP,X)
  261. * x'(i)=0
  262. CALL GSCAL(ZERO,XP)
  263. * b'=r(i)
  264. CALL GCOPY(R,BP)
  265. * mu'=mu
  266. MUP = MU
  267. * WRITE(IOIMP,*) 'ITER=',ITER
  268. * WRITE(IOIMP,*) 'MUP=',MUP
  269. ENDIF
  270. *
  271. * Iteration fails.
  272. *
  273. * IF (ITER.EQ.MAXIT) THEN
  274. IF (INMV.GE.MAXIT) THEN
  275. IRET = 1
  276. IF (IMPR.GT.0) THEN
  277. C WRITE(IOIMP,*)
  278. C $ 'cgsn.eso : Convergence to tolerance not achieved'
  279. C WRITE(IOIMP,*) 'INMV=',INMV,
  280. C $ ' TOL=',TOL,
  281. C $ ' RESID=',RESID
  282. IF (IWARN.GT.0) THEN
  283. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
  284. $ ,' times.'
  285. ENDIF
  286. ENDIF
  287. RETURN
  288. ENDIF
  289. *
  290. * Do some more
  291. *
  292. RHO1 = RHO
  293. GOTO 10
  294. *
  295. * Iteration successful
  296. *
  297. 30 CONTINUE
  298. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  299. WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
  300. ENDIF
  301. IRET = 0
  302. RETURN
  303. *
  304. * End of CGSN.
  305. *
  306. END
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  

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