Télécharger la plaquette Cast3M


fr en

Télécharger bcgg.eso

Retour à la liste

Numérotation des lignes : 

bcgg

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

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