Télécharger bcgg.eso

Retour à la liste

Numérotation des lignes :

bcgg
  1. C BCGG SOURCE PV 20/09/26 21:15:21 10724
  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.  
  55. -INC PPARAM
  56. -INC CCOPTIO
  57. -INC SMLREEL
  58. POINTEUR LRES.MLREEL
  59. -INC SMLENTI
  60. POINTEUR LNMV.MLENTI
  61. POINTEUR TRAV.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 BiCGSTAB(l)
  70. SEGMENT SPACE
  71. REAL*8 IJ(NI,NJ)
  72. ENDSEGMENT
  73. POINTEUR Z.SPACE,ZM1.SPACE
  74. SEGMENT SPACE2
  75. POINTEUR WRK(NIZA).IZA
  76. ENDSEGMENT
  77. POINTEUR R.SPACE2,U.SPACE2
  78. * .. Work Vectors for BiCGStab
  79. POINTEUR RTLD.IZA,XTLD.IZA,BP.IZA
  80. POINTEUR RI.IZA,UI.IZA,UIP.IZA,RIP.IZA
  81. POINTEUR UHAT.IZA,RHAT.IZA,Y.IZA,YP.IZA
  82. *
  83. INTEGER IMPR
  84. * .. Scalar Arguments ..
  85. INTEGER N,NNZ
  86. INTEGER ITER,IRET
  87. REAL*8 BRTOL,RESID
  88. * ..
  89. *
  90. *
  91. * .. Variables locales
  92. * .. Parameters
  93. REAL*8 ZERO,ONE
  94. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  95. * ..
  96. INTEGER MAXIT,III
  97. REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
  98. REAL*8 RHOTOL,OMETOL
  99. REAL*8 GDOT,GNRM2,PETIT,ZMXRNR,ZMXRNX
  100. LOGICAL LCPRES,LUPDAX
  101. * ..
  102. * .. External subroutines and functions..
  103. * ..
  104. * .. Executable Statements ..
  105. *
  106. *
  107. IRET = 0
  108. INMV = 0
  109. IWARN = 0
  110. IGRAN = 0
  111. MGRAN = 3
  112. MAXIT = ITER
  113. TOL = RESID
  114. *
  115. * Same memory space
  116. RHAT=UHAT
  117. *
  118. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bcgg.eso'
  119. *
  120. * Set parameter tolerances.
  121. *
  122. RHOTOL = BRTOL
  123. SIGTOL = BRTOL
  124. RRMAX = 1.D0/(XZPREC**0.25D0)
  125. * XZPREC* 1.D0 plante si tous les termes de la matrice Z sont égaux
  126. * sur IBM
  127. PETIT = XZPREC*10.D0
  128. *
  129. * Set initial residual.
  130. *
  131. C r(1)=b
  132. RI=R.WRK(1)
  133. CALL GCOPY(B,RI)
  134. C r(1)=b-Ax
  135. * IF (GNRM2(X).NE.ZERO) THEN
  136. C r(1)=b-Ax
  137. * CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
  138. * ENDIF
  139. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,RI)
  140. INMV=INMV+1
  141. CALL GCOPY(RI,BP)
  142. *
  143. ITER = 0
  144. BNRM2 = GNRM2(B)
  145. IF(IMPR.GT.6) THEN
  146. WRITE(IOIMP,*) '||B||=',BNRM2
  147. ENDIF
  148. RNRM2 = GNRM2(RI)
  149. IF (ICALRS.EQ.1) BNRM2=RNRM2
  150. *
  151. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  152. RESID = RNRM2 / BNRM2
  153. IF(IMPR.GT.6) THEN
  154. WRITE(IOIMP,*) 'Résidu initial=',RESID
  155. ENDIF
  156. RESMAX = RESID * RRMAX
  157. LRES.PROG(ITER+1)=RESID
  158. LNMV.LECT(ITER+1)=INMV
  159. IF (RESID.LE.TOL) GOTO 30
  160. C
  161. C Variables for reliable updating
  162. C
  163. ZMXRNR=RNRM2
  164. ZMXRNX=RNRM2
  165. LCPRES=.FALSE.
  166. LUPDAX=.FALSE.
  167. XDELTA=1.D-2
  168. C r[tilde]=r(1)
  169. CALL GCOPY(RI,RTLD)
  170. C alpha = rho0 = omega = 1
  171. ALPHA = ONE
  172. RHO0 = ONE
  173. OMEGA = ONE
  174. *
  175. 10 CONTINUE
  176. *
  177. * Perform BiConjugate Gradient part.
  178. *
  179. ITER = ITER + 1
  180. IF(IMPR.GT.7) THEN
  181. WRITE(IOIMP,*) 'ITER=',ITER
  182. ENDIF
  183. *
  184. C rho0 = -omega rho0
  185. RHO0 = -OMEGA * RHO0
  186. DO J = 1 , LBCG
  187. RI = R.WRK(J)
  188. RHO1 = GDOT(RI,RTLD,IDDOT)
  189. IF (ABS(RHO0).LT.RHOTOL) THEN
  190. RHO0=SIGN(XZPREC,RHO0)
  191. IWARN=IWARN+1
  192. ENDIF
  193. BETA = ALPHA * (RHO1 / RHO0)
  194. RHO0 = RHO1
  195. DO I = 1 , J
  196. UI = U.WRK(I)
  197. RI = R.WRK(I)
  198. CALL GSCAL(-BETA,UI)
  199. CALL GAXPY(ONE,RI,UI)
  200. ENDDO
  201. UI = U.WRK(J)
  202. UIP = U.WRK(J+1)
  203. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  204. $ ,UHAT,UI)
  205. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,UHAT,ZERO,UIP)
  206. INMV=INMV+1
  207. C
  208. SIGMA = GDOT(UIP,RTLD,IDDOT)
  209. IF (ABS(SIGMA).LT.SIGTOL) THEN
  210. SIGMA=SIGN(XZPREC,SIGMA)
  211. IWARN=IWARN+1
  212. ENDIF
  213. ALPHA = RHO1 / SIGMA
  214. UI = U.WRK(1)
  215. CALL GAXPY( ALPHA, UI , XTLD)
  216. DO I = 1 , J
  217. UIP = U.WRK(I+1)
  218. RI = R.WRK(I)
  219. CALL GAXPY(-ALPHA,UIP,RI)
  220. ENDDO
  221. RI = R.WRK(J)
  222. RIP = R.WRK(J+1)
  223. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  224. $ ,RHAT,RI)
  225. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,RHAT,ZERO,RIP)
  226. INMV=INMV+1
  227. ENDDO
  228. *
  229. * Perform polynomial part.
  230. *
  231. ZMAX=0.D0
  232. DO I=1,LBCG
  233. DO J=1,I
  234. RI = R.WRK(I+1)
  235. RIP = R.WRK(J+1)
  236. TAU = GDOT (RIP, RI,IDDOT)
  237. Z.IJ(I,J) = TAU
  238. IF (I .NE. J) THEN
  239. Z.IJ(J,I) = TAU
  240. ENDIF
  241. ENDDO
  242. ZMAX=MAX(ZMAX,ABS(Z.IJ(I,I)))
  243. RI = R.WRK(1)
  244. RIP = R.WRK(I+1)
  245. Y.A(I)=GDOT(RI,RIP,IDDOT)
  246. ENDDO
  247. * Ca arrive dans invdiag et invide.dgibi !
  248. IF (ZMAX.LT.(XPETIT**0.5D0)) THEN
  249. ZMAX=XZPREC
  250. ENDIF
  251. * Petite pénalisation des familles pour l'inversibilité
  252. DO I=1,LBCG
  253. Z.IJ(I,I)=Z.IJ(I,I)+PETIT*ZMAX
  254. ENDDO
  255. IIMPR=1
  256. CALL IVMAT(LBCG,Z.IJ,TRAV.LECT,ZM1.IJ,XDETZ,IIMPR,IRET)
  257. * IF (IRET.NE.0) GOTO 25
  258. IF (IRET.NE.0) THEN
  259. * SEGPRT,Z
  260. GOTO 25
  261. ENDIF
  262. *
  263. DO I=1,LBCG
  264. YP.A(I)=ZERO
  265. DO J=1,LBCG
  266. YP.A(I)=YP.A(I)+(ZM1.IJ(I,J)*Y.A(J))
  267. ENDDO
  268. ENDDO
  269. OMEGA=YP.A(LBCG)
  270. DO I=1,LBCG
  271. UI =U.WRK(1)
  272. UIP=U.WRK(I+1)
  273. CALL GAXPY (-YP.A(I),UIP,UI)
  274. RI=R.WRK(I)
  275. CALL GAXPY ( YP.A(I),RI,XTLD)
  276. RI =R.WRK(1)
  277. RIP=R.WRK(I+1)
  278. CALL GAXPY (-YP.A(I),RIP,RI)
  279. ENDDO
  280. *
  281. * Compute residual
  282. *
  283. RI=R.WRK(1)
  284. RNRM2I=GNRM2(RI)
  285. RESID=RNRM2I / BNRM2
  286. IF (RESID.GT.RESMAX) THEN
  287. IGRAN=IGRAN+1
  288. ELSE
  289. IGRAN=0
  290. ENDIF
  291. LRES.PROG(ITER+1)=RESID
  292. LNMV.LECT(ITER+1)=INMV
  293. * WRITE(IOIMP,*) 'ITER=',ITER,' RESID=',RESID
  294. *
  295. * Reliable update
  296. *
  297. ZMXRNX=MAX(RNRM2I,ZMXRNX)
  298. ZMXRNR=MAX(RNRM2I,ZMXRNR)
  299. LUPDAX=(RNRM2I.LT.XDELTA*RNRM2.AND.RNRM2.LT.ZMXRNX)
  300. LCPRES=((RNRM2I.LT.XDELTA*ZMXRNR.AND.RNRM2.LT.ZMXRNR)
  301. $ .OR.LUPDAX)
  302. IF (LCPRES) THEN
  303. RI=R.WRK(1)
  304. CALL GCOPY(BP,RI)
  305. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  306. $ ,UHAT,XTLD)
  307. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,UHAT,ONE,RI)
  308. INMV=INMV+1
  309. ZMXRNR=RNRM2I
  310. *! WRITE(IOIMP,*) 'Compute residual'
  311. IF (LUPDAX) THEN
  312. CALL GAXPY(ONE,UHAT,X)
  313. CALL GSCAL(ZERO,XTLD)
  314. CALL GCOPY(RI,BP)
  315. ZMXRNX=RNRM2I
  316. *! WRITE(IOIMP,*) 'Update approx'
  317. ENDIF
  318. ENDIF
  319. *
  320. * Check for tolerance.
  321. *
  322. IF (RESID.LE.TOL) THEN
  323. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  324. $ ,UHAT,XTLD)
  325. CALL GAXPY(ONE,UHAT,X)
  326. GOTO 30
  327. ENDIF
  328. *
  329. * Iteration fails.
  330. *
  331. * IF (ITER.EQ.MAXIT) THEN
  332. IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  333. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  334. $ ,UHAT,XTLD)
  335. CALL GAXPY(ONE,UHAT,X)
  336. IRET = 1
  337. IF (IMPR.GT.0) THEN
  338. C WRITE(IOIMP,*)
  339. C $ 'bcgg.eso : Convergence to tolerance not achieved'
  340. C WRITE(IOIMP,*) 'INMV=',INMV,
  341. C $ ' TOL=',TOL,
  342. C $ ' RESID=',RESID
  343. IF (IWARN.GT.0) THEN
  344. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  345. $ ,' times.'
  346. ENDIF
  347. ENDIF
  348. RETURN
  349. ENDIF
  350. *
  351. * Do some more
  352. *
  353. GOTO 10
  354. *
  355. * A breakdown in the method has occured.
  356. *
  357. 25 CONTINUE
  358. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  359. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  360. $ ,' times.'
  361. ENDIF
  362. IRET=-1
  363. RETURN
  364. *
  365. * Iteration successful
  366. *
  367. 30 CONTINUE
  368. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  369. WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  370. $ ,' times.'
  371. ENDIF
  372. IRET = 0
  373. RETURN
  374. *
  375. * End of BCGG.
  376. *
  377. END
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  

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