Télécharger bcgg.eso

Retour à la liste

Numérotation des lignes :

  1. C BCGG SOURCE PV 16/11/17 21:58:17 9180
  2. SUBROUTINE BCGG(KMORS,KISA,B,X,
  3. $ KPREC,INVDIA,ILUM,ILUI,
  4. $ LRES,LNMV,
  5. $ RTLD,XTLD,UHAT,R,U,Z,ZM1,Y,YP,BP,TRAV,
  6. $ ITER,INMV,BRTOL,LBCG,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. IMPLICIT INTEGER(I-N)
  9. C***********************************************************************
  10. C NOM : BCGG
  11. C DESCRIPTION :
  12. C Résolution d'un système linéaire Ax=b
  13. C par une méthode BiCGSTAB(l) préconditionnée ou non avec reliable
  14. C updates et convex combination.
  15. C
  16. C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  17. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  18. C mél : gounand@semt2.smts.cea.fr
  19. C REFERENCE (bibtex-like) :
  20. C @BOOK{templates,
  21. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  22. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  23. C H. Van der Vorst},
  24. C TITLE={Templates for the Solution of Linear Systems :
  25. C Building Blocks for Iterative Methods},
  26. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  27. C -> URL : http://www.netlib.org/templates/Templates.html
  28. C@TechReport{fokkema,
  29. C author = {DR Fokkema},
  30. C title = {Enhanced implementation of BiCGSTAB(l) for solving
  31. C linear systems of equations},
  32. C institution = {?},
  33. C year = {1996}}
  34. C***********************************************************************
  35. C ENTREES :
  36. C ENTREES/SORTIES :
  37. C SORTIES :
  38. C CODE RETOUR (IRET) : 0 si ok
  39. C <0 si problème (non utilisé !)
  40. C >0 si l'algorithme n'a pas convergé
  41. C***********************************************************************
  42. C VERSION : v1, 28/02/06, version initiale
  43. C HISTORIQUE : v1, 28/02/06, création
  44. C HISTORIQUE :
  45. C HISTORIQUE :
  46. C***********************************************************************
  47. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  48. C en cas de modification de ce sous-programme afin de faciliter
  49. C la maintenance !
  50. C***********************************************************************
  51. *
  52. * .. Includes .. a commenter si mise au point ok
  53. -INC CCREEL
  54. -INC CCOPTIO
  55. -INC SMLREEL
  56. POINTEUR LRES.MLREEL
  57. -INC SMLENTI
  58. POINTEUR LNMV.MLENTI
  59. POINTEUR TRAV.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 BiCGSTAB(l)
  68. SEGMENT SPACE
  69. REAL*8 IJ(NI,NJ)
  70. ENDSEGMENT
  71. POINTEUR Z.SPACE,ZM1.SPACE
  72. SEGMENT SPACE2
  73. POINTEUR WRK(NIZA).IZA
  74. ENDSEGMENT
  75. POINTEUR R.SPACE2,U.SPACE2
  76. * .. Work Vectors for BiCGStab
  77. POINTEUR RTLD.IZA,XTLD.IZA,BP.IZA
  78. POINTEUR RI.IZA,UI.IZA,UIP.IZA,RIP.IZA
  79. POINTEUR UHAT.IZA,RHAT.IZA,Y.IZA,YP.IZA
  80. *
  81. INTEGER IMPR
  82. * .. Scalar Arguments ..
  83. INTEGER N,NNZ
  84. INTEGER ITER,IRET
  85. REAL*8 BRTOL,RESID
  86. * ..
  87. *
  88. *
  89. * .. Variables locales
  90. * .. Parameters
  91. REAL*8 ZERO,ONE
  92. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  93. * ..
  94. INTEGER MAXIT,III
  95. REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
  96. REAL*8 RHOTOL,OMETOL
  97. REAL*8 GDOT,GNRM2,PETIT,ZMXRNR,ZMXRNX
  98. LOGICAL LCPRES,LUPDAX
  99. * ..
  100. * .. External subroutines and functions..
  101. * ..
  102. * .. Executable Statements ..
  103. *
  104. *
  105. IRET = 0
  106. INMV = 0
  107. IWARN = 0
  108. IGRAN = 0
  109. MGRAN = 3
  110. MAXIT = ITER
  111. TOL = RESID
  112. *
  113. * Same memory space
  114. RHAT=UHAT
  115. *
  116. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgg.eso'
  117. *
  118. * Set parameter tolerances.
  119. *
  120. RHOTOL = BRTOL
  121. SIGTOL = BRTOL
  122. RRMAX = 1.D0/(XZPREC**0.25D0)
  123. * XZPREC* 1.D0 plante si tous les termes de la matrice Z sont égaux
  124. * sur IBM
  125. PETIT = XZPREC*10.D0
  126. *
  127. * Set initial residual.
  128. *
  129. C r(1)=b
  130. RI=R.WRK(1)
  131. CALL GCOPY(B,RI)
  132. C r(1)=b-Ax
  133. * IF (GNRM2(X).NE.ZERO) THEN
  134. C r(1)=b-Ax
  135. * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
  136. * ENDIF
  137. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,RI)
  138. INMV=INMV+1
  139. CALL GCOPY(RI,BP)
  140. *
  141. ITER = 0
  142. BNRM2 = GNRM2(B)
  143. IF(IMPR.GT.6) THEN
  144. WRITE(IOIMP,*) '||B||=',BNRM2
  145. ENDIF
  146. RNRM2 = GNRM2(RI)
  147. IF (ICALRS.EQ.1) BNRM2=RNRM2
  148. *
  149. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  150. RESID = RNRM2 / BNRM2
  151. IF(IMPR.GT.6) THEN
  152. WRITE(IOIMP,*) 'Résidu initial=',RESID
  153. ENDIF
  154. RESMAX = RESID * RRMAX
  155. LRES.PROG(ITER+1)=RESID
  156. LNMV.LECT(ITER+1)=INMV
  157. IF (RESID.LE.TOL) GOTO 30
  158. C
  159. C Variables for reliable updating
  160. C
  161. ZMXRNR=RNRM2
  162. ZMXRNX=RNRM2
  163. LCPRES=.FALSE.
  164. LUPDAX=.FALSE.
  165. XDELTA=1.D-2
  166. C r[tilde]=r(1)
  167. CALL GCOPY(RI,RTLD)
  168. C alpha = rho0 = omega = 1
  169. ALPHA = ONE
  170. RHO0 = ONE
  171. OMEGA = ONE
  172. *
  173. 10 CONTINUE
  174. *
  175. * Perform BiConjugate Gradient part.
  176. *
  177. ITER = ITER + 1
  178. IF(IMPR.GT.7) THEN
  179. WRITE(IOIMP,*) 'ITER=',ITER
  180. ENDIF
  181. *
  182. C rho0 = -omega rho0
  183. RHO0 = -OMEGA * RHO0
  184. DO J = 1 , LBCG
  185. RI = R.WRK(J)
  186. RHO1 = GDOT(RI,RTLD,IDDOT)
  187. IF (ABS(RHO0).LT.RHOTOL) THEN
  188. RHO0=SIGN(XZPREC,RHO0)
  189. IWARN=IWARN+1
  190. ENDIF
  191. BETA = ALPHA * (RHO1 / RHO0)
  192. RHO0 = RHO1
  193. DO I = 1 , J
  194. UI = U.WRK(I)
  195. RI = R.WRK(I)
  196. CALL GSCAL(-BETA,UI)
  197. CALL GAXPY(ONE,RI,UI)
  198. ENDDO
  199. UI = U.WRK(J)
  200. UIP = U.WRK(J+1)
  201. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  202. $ ,UHAT,UI)
  203. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,UHAT,ZERO,UIP)
  204. INMV=INMV+1
  205. C
  206. SIGMA = GDOT(UIP,RTLD,IDDOT)
  207. IF (ABS(SIGMA).LT.SIGTOL) THEN
  208. SIGMA=SIGN(XZPREC,SIGMA)
  209. IWARN=IWARN+1
  210. ENDIF
  211. ALPHA = RHO1 / SIGMA
  212. UI = U.WRK(1)
  213. CALL GAXPY( ALPHA, UI , XTLD)
  214. DO I = 1 , J
  215. UIP = U.WRK(I+1)
  216. RI = R.WRK(I)
  217. CALL GAXPY(-ALPHA,UIP,RI)
  218. ENDDO
  219. RI = R.WRK(J)
  220. RIP = R.WRK(J+1)
  221. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  222. $ ,RHAT,RI)
  223. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,RHAT,ZERO,RIP)
  224. INMV=INMV+1
  225. ENDDO
  226. *
  227. * Perform polynomial part.
  228. *
  229. ZMAX=0.D0
  230. DO I=1,LBCG
  231. DO J=1,I
  232. RI = R.WRK(I+1)
  233. RIP = R.WRK(J+1)
  234. TAU = GDOT (RIP, RI,IDDOT)
  235. Z.IJ(I,J) = TAU
  236. IF (I .NE. J) THEN
  237. Z.IJ(J,I) = TAU
  238. ENDIF
  239. ENDDO
  240. ZMAX=MAX(ZMAX,ABS(Z.IJ(I,I)))
  241. RI = R.WRK(1)
  242. RIP = R.WRK(I+1)
  243. Y.A(I)=GDOT(RI,RIP,IDDOT)
  244. ENDDO
  245. * Ca arrive dans invdiag et invide.dgibi !
  246. IF (ZMAX.LT.(XPETIT**0.5D0)) THEN
  247. ZMAX=XZPREC
  248. ENDIF
  249. * Petite pénalisation des familles pour l'inversibilité
  250. DO I=1,LBCG
  251. Z.IJ(I,I)=Z.IJ(I,I)+PETIT*ZMAX
  252. ENDDO
  253. IIMPR=1
  254. CALL IVMAT(LBCG,Z.IJ,TRAV.LECT,ZM1.IJ,XDETZ,IIMPR,IRET)
  255. * IF (IRET.NE.0) GOTO 25
  256. IF (IRET.NE.0) THEN
  257. * SEGPRT,Z
  258. GOTO 25
  259. ENDIF
  260. *
  261. DO I=1,LBCG
  262. YP.A(I)=ZERO
  263. DO J=1,LBCG
  264. YP.A(I)=YP.A(I)+(ZM1.IJ(I,J)*Y.A(J))
  265. ENDDO
  266. ENDDO
  267. OMEGA=YP.A(LBCG)
  268. DO I=1,LBCG
  269. UI =U.WRK(1)
  270. UIP=U.WRK(I+1)
  271. CALL GAXPY (-YP.A(I),UIP,UI)
  272. RI=R.WRK(I)
  273. CALL GAXPY ( YP.A(I),RI,XTLD)
  274. RI =R.WRK(1)
  275. RIP=R.WRK(I+1)
  276. CALL GAXPY (-YP.A(I),RIP,RI)
  277. ENDDO
  278. *
  279. * Compute residual
  280. *
  281. RI=R.WRK(1)
  282. RNRM2I=GNRM2(RI)
  283. RESID=RNRM2I / BNRM2
  284. IF (RESID.GT.RESMAX) THEN
  285. IGRAN=IGRAN+1
  286. ELSE
  287. IGRAN=0
  288. ENDIF
  289. LRES.PROG(ITER+1)=RESID
  290. LNMV.LECT(ITER+1)=INMV
  291. * WRITE(IOIMP,*) 'ITER=',ITER,' RESID=',RESID
  292. *
  293. * Reliable update
  294. *
  295. ZMXRNX=MAX(RNRM2I,ZMXRNX)
  296. ZMXRNR=MAX(RNRM2I,ZMXRNR)
  297. LUPDAX=(RNRM2I.LT.XDELTA*RNRM2.AND.RNRM2.LT.ZMXRNX)
  298. LCPRES=((RNRM2I.LT.XDELTA*ZMXRNR.AND.RNRM2.LT.ZMXRNR)
  299. $ .OR.LUPDAX)
  300. IF (LCPRES) THEN
  301. RI=R.WRK(1)
  302. CALL GCOPY(BP,RI)
  303. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  304. $ ,UHAT,XTLD)
  305. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,UHAT,ONE,RI)
  306. INMV=INMV+1
  307. ZMXRNR=RNRM2I
  308. *! WRITE(IOIMP,*) 'Compute residual'
  309. IF (LUPDAX) THEN
  310. CALL GAXPY(ONE,UHAT,X)
  311. CALL GSCAL(ZERO,XTLD)
  312. CALL GCOPY(RI,BP)
  313. ZMXRNX=RNRM2I
  314. *! WRITE(IOIMP,*) 'Update approx'
  315. ENDIF
  316. ENDIF
  317. *
  318. * Check for tolerance.
  319. *
  320. IF (RESID.LE.TOL) THEN
  321. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  322. $ ,UHAT,XTLD)
  323. CALL GAXPY(ONE,UHAT,X)
  324. GOTO 30
  325. ENDIF
  326. *
  327. * Iteration fails.
  328. *
  329. * IF (ITER.EQ.MAXIT) THEN
  330. IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  331. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  332. $ ,UHAT,XTLD)
  333. CALL GAXPY(ONE,UHAT,X)
  334. IRET = 1
  335. IF (IMPR.GT.0) THEN
  336. C WRITE(IOIMP,*)
  337. C $ 'bcgg.eso : Convergence to tolerance not achieved'
  338. C WRITE(IOIMP,*) 'INMV=',INMV,
  339. C $ ' TOL=',TOL,
  340. C $ ' RESID=',RESID
  341. IF (IWARN.GT.0) THEN
  342. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  343. $ ,' times.'
  344. ENDIF
  345. ENDIF
  346. RETURN
  347. ENDIF
  348. *
  349. * Do some more
  350. *
  351. GOTO 10
  352. *
  353. * A breakdown in the method has occured.
  354. *
  355. 25 CONTINUE
  356. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  357. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  358. $ ,' times.'
  359. ENDIF
  360. IRET=-1
  361. RETURN
  362. *
  363. * Iteration successful
  364. *
  365. 30 CONTINUE
  366. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  367. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  368. $ ,' times.'
  369. ENDIF
  370. IRET = 0
  371. RETURN
  372. *
  373. * End of BCGG.
  374. *
  375. END
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  

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