Télécharger bicg.eso

Retour à la liste

Numérotation des lignes :

bicg
  1. C BICG SOURCE PV 20/09/26 21:15:22 10724
  2. SUBROUTINE BICG(KMORS,KISA,B,X,
  3. $ KPREC,INVDIA,ILUM,ILUI,
  4. $ LRES,LNMV,
  5. $ R,RTLD,Z,ZTLD,P,PTLD,Q,QTLD,
  6. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. C***********************************************************************
  10. C NOM : BICG
  11. C DESCRIPTION :
  12. C Résolution d'un système linéaire Ax=b
  13. C par la méthode du bi-gradient conjugué préconditionné.
  14. C
  15. C LANGAGE : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  16. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  17. C mél : gounand@semt2.smts.cea.fr
  18. C REFERENCE (bibtex-like) :
  19. C @BOOK{templates,
  20. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  21. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  22. C H. Van der Vorst},
  23. C TITLE={Templates for the Solution of Linear Systems :
  24. C Building Blocks for Iterative Methods},
  25. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  26. C -> URL : http://www.netlib.org/templates/Templates.html
  27. C***********************************************************************
  28. C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
  29. C GNRM2, GDOT (functions)
  30. C APPELES (Calcul) : GMOMV : produit Matrice Morse - Vecteur
  31. C PSOLVE : Preconditionner solve
  32. C***********************************************************************
  33. C ENTREES :
  34. C ENTREES/SORTIES :
  35. C SORTIES :
  36. C CODE RETOUR (IRET) : 0 si ok
  37. C <0 si problème (non utilisé !)
  38. C >0 si l'algorithme n'a pas convergé
  39. C***********************************************************************
  40. C VERSION : v1, 09/02/06
  41. C HISTORIQUE : v1, 09/02/06, création
  42. C HISTORIQUE :
  43. C HISTORIQUE :
  44. C***********************************************************************
  45. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  46. C en cas de modification de ce sous-programme afin de faciliter
  47. C la maintenance !
  48. C***********************************************************************
  49. *
  50. * .. Includes .. a commenter si mise au point ok
  51.  
  52. -INC PPARAM
  53. -INC CCOPTIO
  54. -INC CCREEL
  55. -INC SMLREEL
  56. POINTEUR LRES.MLREEL
  57. -INC SMLENTI
  58. POINTEUR LNMV.MLENTI
  59. POINTEUR KMORS.PMORS
  60. POINTEUR KISA.IZA
  61. POINTEUR INVDIA.IZA
  62. POINTEUR ILUM.PMORS
  63. POINTEUR ILUI.IZA
  64. POINTEUR B.IZA
  65. POINTEUR X.IZA
  66. * .. Work Vectors for CG
  67. POINTEUR R.IZA,P.IZA,Q.IZA,Z.IZA
  68. POINTEUR RTLD.IZA,PTLD.IZA,QTLD.IZA,ZTLD.IZA
  69. * .. Scalar Arguments ..
  70. INTEGER ITER, KPREC, IMPR, IRET
  71. REAL*8 RESID
  72. *
  73. *
  74. * .. Variables locales
  75. * .. Parameters
  76. REAL*8 ZERO , ONE
  77. PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  78. * ..
  79. INTEGER MAXIT
  80. REAL*8 ALPHA,BETA,RHO,RHO1,TOL,GNRM2,OMEGA
  81. REAL*8 GDOT,BNRM2
  82. * ..
  83. * .. External subroutines and functions..
  84. EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
  85. * ..
  86. * .. Executable Statements ..
  87. *
  88.  
  89. IRET = 0
  90. INMV = 0
  91. IWARN = 0
  92. IGRAN = 0
  93. MGRAN = 5
  94. MAXIT = ITER
  95. TOL = RESID
  96. IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bicg.eso'
  97. *
  98. * Set parameter tolerances.
  99. *
  100. RHOTOL = BRTOL
  101. OMETOL = BRTOL
  102. RRMAX = 1.D0/(XZPREC**0.25D0)
  103. *
  104. * Calcul du résidu initial
  105. *
  106. C r(0)=b
  107. CALL GCOPY(B,R)
  108.  
  109. C r(0)=b-Ax
  110. CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  111. INMV=INMV+1
  112. *
  113. ITER = 0
  114. BNRM2 = GNRM2(B)
  115. IF(IMPR.GT.5) THEN
  116. WRITE(IOIMP,*) '||B||=',BNRM2
  117. ENDIF
  118. RNRM2 = GNRM2(R)
  119. IF (ICALRS.EQ.1) BNRM2=RNRM2
  120. IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  121. RESID = RNRM2 / BNRM2
  122. IF(IMPR.GT.5) THEN
  123. WRITE(IOIMP,*) 'Initial residual=',RESID
  124. ENDIF
  125. RESMAX = RESID * RRMAX
  126. LRES.PROG(ITER+1)=RESID
  127. LNMV.LECT(ITER+1)=INMV
  128. *
  129. IF (RESID.LE.TOL) THEN
  130. IRET = 0
  131. GOTO 30
  132. ENDIF
  133. CALL GCOPY(R,RTLD)
  134. *
  135. 10 CONTINUE
  136. *
  137. * Perform (Preconditioned) BiConjugate Gradient Iteration.
  138. *
  139. ITER = ITER +1
  140. IF(IMPR.GT.7) THEN
  141. WRITE(IOIMP,*) 'ITER=',ITER
  142. ENDIF
  143. *
  144. * Preconditioner Solve
  145. C Mz(i-1)=r(i-1)
  146. CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  147. *! CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  148. $ ,Z,R)
  149. C M'z'(i-1)=r'(i-1)
  150. *! CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  151. CALL PSOLVE('T',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  152. $ ,ZTLD,RTLD)
  153. C rho(i-1)=z(i-1)'r'(i-1)
  154. RHO = GDOT(Z,RTLD,IDDOT)
  155. C WRITE(IOIMP,*) ' RHO=',RHO
  156. C XR1 = GDOT(Z,R)
  157. C WRITE(IOIMP,*) ' XR1=',XR1
  158. C XR2 = GDOT(ZTLD,RTLD)
  159. C WRITE(IOIMP,*) ' XR2=',XR2
  160. IF (ABS(RHO).LT.RHOTOL) THEN
  161. RHO=SIGN(XZPREC,RHO)
  162. IWARN=IWARN+1
  163. ENDIF
  164. *
  165. * Compute direction vector P
  166. *
  167. IF (ITER.GT.1) THEN
  168. C beta(i-1)=rho(i-1)/rho(i-2)
  169. BETA = RHO / RHO1
  170. C p(i)=z(i-1)+beta(i-1)p(i-1)
  171. C p'(i)=z'(i-1)+beta(i-1)p'(i-1)
  172. CALL GAXPY(BETA,P,Z)
  173. CALL GCOPY(Z,P)
  174. CALL GAXPY(BETA,PTLD,ZTLD)
  175. CALL GCOPY(ZTLD,PTLD)
  176. ELSE
  177. C p(1)=z(0)
  178. C p'(1)=z'(0)
  179. CALL GCOPY(Z,P)
  180. CALL GCOPY(ZTLD,PTLD)
  181. ENDIF
  182. *
  183. * Compute scalar ALPHA (save A*p to q)
  184. *
  185. C q(i)=Ap(i)
  186. CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,P,ZERO,Q)
  187. INMV=INMV+1
  188. C q'(i)=A'p'(i)
  189. CALL GMOMV(IMVEC,'T',ONE,KMORS,KISA,PTLD,ZERO,QTLD)
  190. INMV=INMV+1
  191. C alpha(i)=rho(i-1)/(p'(i)'q(i))
  192. OMEGA = GDOT(PTLD,Q,IDDOT)
  193. C WRITE(IOIMP,*) ' OMEGA=',OMEGA
  194. C XOM1 = GDOT(P,Q)
  195. C WRITE(IOIMP,*) ' XOM1 =',XOM1
  196. C XOM2 = GDOT(PTLD,QTLD)
  197. C WRITE(IOIMP,*) ' XOM2 =',XOM2
  198. IF (ABS(OMEGA).LT.OMETOL) THEN
  199. OMEGA=SIGN(XZPREC,OMEGA)
  200. IWARN=IWARN+1
  201. ENDIF
  202. ALPHA = RHO / OMEGA
  203. *
  204. * Compute current solution vector x
  205. *
  206. C x(i)=x(i-1)+alpha(i)p(i)
  207. CALL GAXPY(ALPHA,P,X)
  208. *
  209. * Compute residual vector R, find norm,
  210. * then check for tolerance
  211. *
  212. C r(i)=r(i-1)-alpha(i)q(i)
  213. CALL GAXPY(-ALPHA,Q,R)
  214. *
  215. RESID = GNRM2(R) / BNRM2
  216. IF(IMPR.GT.7) THEN
  217. WRITE(IOIMP,*) 'Résidu=',RESID
  218. ENDIF
  219. LRES.PROG(ITER+1)=RESID
  220. LNMV.LECT(ITER+1)=INMV
  221. IF (RESID.GT.RESMAX) THEN
  222. IGRAN=IGRAN+1
  223. ELSE
  224. IGRAN=0
  225. ENDIF
  226. *
  227. * Iteration successful
  228. *
  229. IF (RESID.LE.TOL) THEN
  230. IRET = 0
  231. *
  232. * Iteration fails
  233. *
  234. * ELSE IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  235. ELSE IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  236. IRET = 1
  237. *
  238. * Do some more
  239. *
  240. ELSE
  241. C r'(i)=r'(i-1)-alpha(i)q'(i)
  242. CALL GAXPY(-ALPHA,QTLD,RTLD)
  243. RHO1 = RHO
  244. GOTO 10
  245. ENDIF
  246. *
  247. * Normal termination
  248. *
  249. 30 CONTINUE
  250. IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  251. WRITE(IOIMP,*) 'BiCG nearly failed ',IWARN
  252. $ ,' times.'
  253. ENDIF
  254. RETURN
  255. *
  256. * End of BICG
  257. *
  258. END
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  

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