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   25/04/30    21:15:01     12258          
  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**0.5D0
  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.             ZMAX=MAX(ZMAX,ABS(Z.IJ(J,I)))
  246.          ENDDO
  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. * Mettre plus que SQRT(XPETIT) car c'est le critere utilise dans
  267. *     ivmat pour la nullite
  268.          ZMAX2=MAX(ZMAX*PETIT,(XPETIT**0.33D0))
  269. *         write(ioimp,*) 'ZMAX,PETIT,XPETIT,ZMAX2=',ZMAX,PETIT,XPETIT
  270. *     $        ,ZMAX2
  271. *     Petite pénalisation des familles pour l'inversibilité
  272.          DO I=1,LBCG
  273.             Z.IJ(I,I)=Z.IJ(I,I)+ZMAX2
  274.          ENDDO
  275.          CALL IVMAT(LBCG,Z.IJ,TRAV.LECT,ZM1.IJ,XDETZ,IIMPR,IRET)
  276.          IF (IRET.NE.0) THEN
  277.             write(ioimp,*) 'Matrice penalized with ZMAX2=',ZMAX2
  278.             MOTERR(1:8)='bcgg.eso'
  279.             CALL ERREUR(700)
  280.             RETURN
  281.          ENDIF
  282.       ENDIF
  283. *
  284.       DO I=1,LBCG
  285.          YP.A(I)=ZERO
  286.          DO J=1,LBCG
  287.             YP.A(I)=YP.A(I)+(ZM1.IJ(I,J)*Y.A(J))
  288.          ENDDO
  289.       ENDDO
  290.       OMEGA=YP.A(LBCG)
  291.       DO I=1,LBCG
  292.          UI =U.WRK(1)
  293.          UIP=U.WRK(I+1)
  294.          CALL GAXPY (-YP.A(I),UIP,UI)
  295.          RI=R.WRK(I)
  296.          CALL GAXPY ( YP.A(I),RI,XTLD)
  297.          RI =R.WRK(1)
  298.          RIP=R.WRK(I+1)
  299.          CALL GAXPY (-YP.A(I),RIP,RI)
  300.       ENDDO
  301. *
  302. *   Compute residual
  303. *
  304.       RI=R.WRK(1)
  305.       RNRM2I=GNRM2(RI)
  306.       RESID=RNRM2I / BNRM2
  307.       IF (RESID.GT.RESMAX) THEN
  308.          IGRAN=IGRAN+1
  309.       ELSE
  310.          IGRAN=0
  311.       ENDIF
  312.       LRES.PROG(ITER+1)=RESID
  313.       LNMV.LECT(ITER+1)=INMV
  314. *      WRITE(IOIMP,*) 'ITER=',ITER,' RESID=',RESID
  315. *
  316. * Reliable update
  317. *
  318.       ZMXRNX=MAX(RNRM2I,ZMXRNX)
  319.       ZMXRNR=MAX(RNRM2I,ZMXRNR)
  320.       LUPDAX=(RNRM2I.LT.XDELTA*RNRM2.AND.RNRM2.LT.ZMXRNX)
  321.       LCPRES=((RNRM2I.LT.XDELTA*ZMXRNR.AND.RNRM2.LT.ZMXRNR)
  322.      $     .OR.LUPDAX)
  323.       IF (LCPRES) THEN
  324.          RI=R.WRK(1)
  325.          CALL GCOPY(BP,RI)
  326.          CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  327.      $        ,UHAT,XTLD)
  328.          CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,UHAT,ONE,RI)
  329.          INMV=INMV+1
  330.          ZMXRNR=RNRM2I
  331. *!         WRITE(IOIMP,*) 'Compute residual'
  332.          IF (LUPDAX) THEN
  333.             CALL GAXPY(ONE,UHAT,X)
  334.             CALL GSCAL(ZERO,XTLD)
  335.             CALL GCOPY(RI,BP)
  336.             ZMXRNX=RNRM2I
  337. *!            WRITE(IOIMP,*) 'Update approx'
  338.          ENDIF
  339.       ENDIF
  340. *
  341. *   Check for tolerance.
  342. *
  343.       IF (RESID.LE.TOL) THEN
  344.          CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  345.      $        ,UHAT,XTLD)
  346.          CALL GAXPY(ONE,UHAT,X)
  347.          GOTO 30
  348.       ENDIF
  349. *
  350. *     Iteration fails.
  351. *
  352. *      IF (ITER.EQ.MAXIT) THEN
  353.       IF (INMV.GE.MAXIT.OR.IGRAN.GT.MGRAN) THEN
  354.          CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  355.      $        ,UHAT,XTLD)
  356.          CALL GAXPY(ONE,UHAT,X)
  357.          IRET = 1
  358.          IF (IMPR.GT.0) THEN
  359. C            WRITE(IOIMP,*)
  360. C     $           'bcgg.eso : Convergence to tolerance not achieved'
  361. C            WRITE(IOIMP,*) 'INMV=',INMV,
  362. C     $           ' TOL=',TOL,
  363. C     $           ' RESID=',RESID
  364.             IF (IWARN.GT.0) THEN
  365.                WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  366.      $              ,' times.'
  367.             ENDIF
  368.          ENDIF
  369.          RETURN
  370.       ENDIF
  371. *
  372. *     Do some more
  373. *
  374.       GOTO 10
  375. *
  376. *     A breakdown in the method has occured.
  377. *
  378.  25   CONTINUE
  379.       IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  380.          WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  381.      $        ,' times.'
  382.       ENDIF
  383.       IRET=-1
  384.       RETURN
  385. *
  386. *     Iteration successful
  387. *
  388.  30   CONTINUE
  389.       IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  390.          WRITE(IOIMP,*) 'BiCGStab(l)G nearly failed ',IWARN
  391.      $        ,' times.'
  392.       ENDIF
  393.       IRET = 0
  394.       RETURN
  395. *
  396. *     End of BCGG.
  397. *
  398.       END
  399.  
  400.  

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