Télécharger la plaquette Cast3M


fr en

Télécharger bicg.eso

Retour à la liste

Numérotation des lignes : 

bicg

  1. C BICG      SOURCE    GOUNAND   22/08/25    21:15:02     11434          
  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
  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 CG
  68.       POINTEUR R.IZA,P.IZA,Q.IZA,Z.IZA
  69.       POINTEUR RTLD.IZA,PTLD.IZA,QTLD.IZA,ZTLD.IZA
  70. *     .. Scalar Arguments ..
  71.       INTEGER ITER, KPREC, IMPR, IRET
  72.       REAL*8  RESID
  73. *
  74. *
  75. *     .. Variables locales
  76. *     .. Parameters
  77.       REAL*8     ZERO      , ONE
  78.       PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  79. *     ..
  80.       INTEGER MAXIT
  81.       REAL*8 ALPHA,BETA,RHO,RHO1,TOL,GNRM2,OMEGA
  82.       REAL*8 GDOT,BNRM2
  83. *     ..
  84. *     .. External subroutines and functions..
  85.       EXTERNAL GAXPY,GCOPY,GDOT,GNRM2
  86. *     ..
  87. *     .. Executable Statements ..
  88. *
  89.  
  90.       IRET  = 0
  91.       INMV  = 0
  92.       IWARN = 0
  93.       IGRAN = 0
  94.       MGRAN = 5
  95.       MAXIT = ITER
  96.       TOL   = RESID
  97.       IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans bicg.eso'
  98. *
  99. *     Set parameter tolerances.
  100. *
  101.       RHOTOL = BRTOL
  102.       OMETOL = BRTOL
  103.       RRMAX = 1.D0/(XZPREC**0.25D0)
  104. *
  105. *     Calcul du résidu initial
  106. *
  107. C                    r(0)=b
  108.       CALL GCOPY(B,R)
  109.  
  110. C                    r(0)=b-Ax
  111.       CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  112.       INMV=INMV+1
  113. *
  114.       ITER = 0
  115.       BNRM2 = GNRM2(B)
  116.       IF(IMPR.GT.5) THEN
  117.          WRITE(IOIMP,*) '||B||=',BNRM2
  118.       ENDIF
  119.       RNRM2 = GNRM2(R)
  120.       IF (ICALRS.EQ.1) BNRM2=RNRM2
  121.       IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  122.       RESID = RNRM2 / BNRM2
  123.       IF(IMPR.GT.5) THEN
  124.          WRITE(IOIMP,*) 'Initial residual=',RESID
  125.       ENDIF
  126.       RESMAX = RESID * RRMAX
  127.       LRES.PROG(ITER+1)=RESID
  128.       LNMV.LECT(ITER+1)=INMV
  129. *
  130.       IF (RESID.LE.TOL) THEN
  131.          IRET = 0
  132.          GOTO 30
  133.       ENDIF
  134.       CALL GCOPY(R,RTLD)
  135. *
  136.  10   CONTINUE
  137. * Gestion du CTRL-C
  138.       if (ierr.NE.0) return
  139. *
  140. *     Perform (Preconditioned) BiConjugate Gradient Iteration.
  141. *
  142.       ITER = ITER +1
  143.       IF(IMPR.GT.7) THEN
  144.          WRITE(IOIMP,*) 'ITER=',ITER
  145.       ENDIF
  146. *
  147. *     Preconditioner Solve
  148. C                    Mz(i-1)=r(i-1)
  149.       CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  150. *!      CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  151.      $     ,Z,R)
  152. C                    M'z'(i-1)=r'(i-1)
  153. *!      CALL PSOLVE(KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  154.       CALL PSOLVE('T',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  155.      $     ,ZTLD,RTLD)
  156. C                   rho(i-1)=z(i-1)'r'(i-1)
  157.       RHO = GDOT(Z,RTLD,IDDOT)
  158. C      WRITE(IOIMP,*) '  RHO=',RHO
  159. C      XR1 = GDOT(Z,R)
  160. C      WRITE(IOIMP,*) '  XR1=',XR1
  161. C      XR2 = GDOT(ZTLD,RTLD)
  162. C      WRITE(IOIMP,*) '  XR2=',XR2
  163.       IF (ABS(RHO).LT.RHOTOL) THEN
  164.          RHO=SIGN(XZPREC,RHO)
  165.          IWARN=IWARN+1
  166.       ENDIF
  167. *
  168. *     Compute direction vector P
  169. *
  170.       IF (ITER.GT.1) THEN
  171. C                   beta(i-1)=rho(i-1)/rho(i-2)
  172.          BETA = RHO / RHO1
  173. C                   p(i)=z(i-1)+beta(i-1)p(i-1)
  174. C                   p'(i)=z'(i-1)+beta(i-1)p'(i-1)
  175.          CALL GAXPY(BETA,P,Z)
  176.          CALL GCOPY(Z,P)
  177.          CALL GAXPY(BETA,PTLD,ZTLD)
  178.          CALL GCOPY(ZTLD,PTLD)
  179.       ELSE
  180. C                   p(1)=z(0)
  181. C                   p'(1)=z'(0)
  182.          CALL GCOPY(Z,P)
  183.          CALL GCOPY(ZTLD,PTLD)
  184.       ENDIF
  185. *
  186. *     Compute scalar ALPHA (save A*p to q)
  187. *
  188. C                   q(i)=Ap(i)
  189.       CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,P,ZERO,Q)
  190.       INMV=INMV+1
  191. C                   q'(i)=A'p'(i)
  192.       CALL GMOMV(IMVEC,'T',ONE,KMORS,KISA,PTLD,ZERO,QTLD)
  193.       INMV=INMV+1
  194. C                   alpha(i)=rho(i-1)/(p'(i)'q(i))
  195.       OMEGA = GDOT(PTLD,Q,IDDOT)
  196. C      WRITE(IOIMP,*) '  OMEGA=',OMEGA
  197. C      XOM1 = GDOT(P,Q)
  198. C      WRITE(IOIMP,*) '  XOM1 =',XOM1
  199. C      XOM2 = GDOT(PTLD,QTLD)
  200. C      WRITE(IOIMP,*) '  XOM2 =',XOM2
  201.       IF (ABS(OMEGA).LT.OMETOL) THEN
  202.          OMEGA=SIGN(XZPREC,OMEGA)
  203.          IWARN=IWARN+1
  204.       ENDIF
  205.       ALPHA = RHO / OMEGA
  206. *
  207. *     Compute current solution vector x
  208. *
  209. C                   x(i)=x(i-1)+alpha(i)p(i)
  210.       CALL GAXPY(ALPHA,P,X)
  211. *
  212. *     Compute residual vector R, find norm,
  213. *     then check for tolerance
  214. *
  215. C                   r(i)=r(i-1)-alpha(i)q(i)
  216.       CALL GAXPY(-ALPHA,Q,R)
  217. *
  218.       RESID = GNRM2(R) / BNRM2
  219.       IF(IMPR.GT.7) THEN
  220.          WRITE(IOIMP,*) 'Résidu=',RESID
  221.       ENDIF
  222.       LRES.PROG(ITER+1)=RESID
  223.       LNMV.LECT(ITER+1)=INMV
  224.       IF (RESID.GT.RESMAX) THEN
  225.          IGRAN=IGRAN+1
  226.       ELSE
  227.          IGRAN=0
  228.       ENDIF
  229. *
  230. *     Iteration successful
  231. *
  232.       IF (RESID.LE.TOL) THEN
  233.          IRET = 0
  234. *
  235. *     Iteration fails
  236. *
  237. *      ELSE IF (ITER.EQ.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  238.       ELSE IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  239.          IRET = 1
  240. *
  241. *     Do some more
  242. *
  243.       ELSE
  244. C                   r'(i)=r'(i-1)-alpha(i)q'(i)
  245.          CALL GAXPY(-ALPHA,QTLD,RTLD)
  246.          RHO1 = RHO
  247.          GOTO 10
  248.       ENDIF
  249. *
  250. *     Normal termination
  251. *
  252.  30   CONTINUE
  253.       IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  254.          WRITE(IOIMP,*) 'BiCG nearly failed ',IWARN
  255.      $        ,' times.'
  256.       ENDIF
  257.       RETURN
  258. *
  259. *     End of BICG
  260. *
  261.       END
  262.  
  263.  

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