Télécharger la plaquette Cast3M


fr en

Télécharger cgsn.eso

Retour à la liste

Numérotation des lignes : 

cgsn

  1. C CGSN      SOURCE    GOUNAND   22/08/25    21:15:03     11434          
  2.       SUBROUTINE CGSN(KMORS,KISA,B,X,
  3.      $     KPREC,INVDIA,ILUM,ILUI,
  4.      $     LRES,LNMV,
  5.      $     R,RTLD,P,PHAT,
  6.      $     Q,QHAT,UHAT,
  7.      $     XP,BP,
  8.      $     ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  9.       IMPLICIT REAL*8 (A-H,O-Z)
  10.       IMPLICIT INTEGER(I-N)
  11. C***********************************************************************
  12. C NOM         : CGSN
  13. C DESCRIPTION :
  14. C     Résolution d'un système linéaire Ax=b
  15. C     par une méthode de gradient conjugué au carré (CGS) avec la
  16. C     stratégie de stabilisation de Neumaier préconditionnée ou non.
  17. C
  18. C LANGAGE     : FORTRAN 77 + chouhia ESOPE (pour les E/S)
  19. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  20. C               mél : gounand@semt2.smts.cea.fr
  21. C REFERENCE (bibtex-like) :
  22. C  @BOOK{templates,
  23. C        AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  24. C                J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  25. C                H. Van der Vorst},
  26. C        TITLE={Templates for the Solution of Linear Systems :
  27. C               Building Blocks for Iterative Methods},
  28. C        PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  29. C  -> URL : http://www.netlib.org/templates/Templates.html
  30. C***********************************************************************
  31. C APPELES (BLAS 1) : GCOPY, GAXPY (subroutines)
  32. C                    GNRM2, GDOT  (functions)
  33. C APPELES (Calcul) : GMOMV  : produit Matrice Morse - Vecteur
  34. C                    PSOLVE : Preconditionner solve
  35. C***********************************************************************
  36. C ENTREES            :
  37. C ENTREES/SORTIES    :
  38. C SORTIES            :
  39. C CODE RETOUR (IRET) : 0 si ok
  40. C                     <0 si problème (non utilisé !)
  41. C                     >0 si l'algorithme n'a pas convergé
  42. C***********************************************************************
  43. C VERSION    : v1, 10/02/06
  44. C HISTORIQUE : v1, 10/02/06, création
  45. C HISTORIQUE :
  46. C HISTORIQUE :
  47. C***********************************************************************
  48. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  49. C en cas de modification de ce sous-programme afin de faciliter
  50. C la maintenance !
  51. C***********************************************************************
  52. *
  53. *     .. Includes .. a commenter si mise au point ok
  54. -INC CCREEL
  55.  
  56. -INC PPARAM
  57. -INC CCOPTIO
  58. -INC SMLREEL
  59.       POINTEUR LRES.MLREEL
  60. -INC SMLENTI
  61.       POINTEUR LNMV.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 CGSN
  71.       POINTEUR R.IZA,RTLD.IZA,P.IZA,PHAT.IZA
  72.       POINTEUR Q.IZA,QHAT.IZA,UHAT.IZA
  73.       POINTEUR U.IZA,VHAT.IZA
  74. *  Note QHAT and U, UHAT,VHAT share the same memory space
  75.       POINTEUR XP.IZA,BP.IZA
  76. *
  77.       INTEGER IMPR
  78. *     .. Scalar Arguments ..
  79.       INTEGER  N,NNZ
  80.       INTEGER  ITER,IRET
  81.       REAL*8   BRTOL,RESID
  82. *     ..
  83. *
  84. *
  85. *     .. Variables locales
  86. *     .. Parameters
  87.       REAL*8 ZERO,ONE
  88.       PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  89. *     ..
  90.       INTEGER MAXIT,III
  91.       REAL*8 TOL,ALPHA,BETA,RHO,RHO1,BNRM2,OMEGA
  92.       REAL*8 RHOTOL
  93.       REAL*8 GDOT,GNRM2,RNRM2
  94.       REAL*8 MU,MUP
  95. *     ..
  96. *     .. External subroutines and functions..
  98. *     ..
  99. *     .. Executable Statements ..
  100. *
  101. *
  102.       IRET  = 0
  103.       INMV  = 0
  104.       IWARN = 0
  105.       MAXIT = ITER
  106.       TOL   = RESID
  107. *
  108. *     Same memory space
  109.       U=QHAT
  110.       VHAT=UHAT
  111. *
  112.       IF(IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans cgsn.eso'
  113. *
  114. *     Set breakdown parameter tolerances.
  115. *
  116.       RHOTOL = BRTOL
  117.       TAUTOL = BRTOL
  118. *
  119. *     Set initial residual.
  120. *
  121. C                    r(0)=b
  122.       CALL GCOPY(B,R)
  123. C                    r(0)=b-Ax
  124. *      IF (GNRM2(X).NE.ZERO) THEN
  125. C                    r(0)=b-Ax
  126. *         CALL GMOMV ('N',-ONE,KMORS,KISA,X,ONE,R)
  127. *      ENDIF
  128.       CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,X,ONE,R)
  129.       INMV=INMV+1
  130. *
  131.       ITER = 0
  132.       BNRM2 = GNRM2(B)
  133.       IF(IMPR.GT.6) THEN
  134.          WRITE(IOIMP,*) '||B||=',BNRM2
  135.       ENDIF
  136.       RNRM2 = GNRM2(R)
  137.       IF (ICALRS.EQ.1) BNRM2=RNRM2
  138. *
  139.       IF (BNRM2.LT.XPETIT) BNRM2 = ONE
  140.       RESID = RNRM2 / BNRM2
  141.       IF(IMPR.GT.6) THEN
  142.          WRITE(IOIMP,*) 'Résidu initial=',RESID
  143.       ENDIF
  144.       LRES.PROG(ITER+1)=RESID
  145.       LNMV.LECT(ITER+1)=INMV
  146.       IF (RESID.LE.TOL) GOTO 30
  147. C                    r[tilde]=r(0)
  148.       CALL GCOPY(R,RTLD)
  149. C                    Neumaier's Strategy
  150. C                    x=x(0)  x'=0
  151. C          Normalement, c'est déjà le cas
  152. C                    b'=r(0)   mu' = ||b'||
  153.       CALL GCOPY(R,BP)
  154.       MUP=RNRM2
  155. *!      WRITE(IOIMP,*) 'MUP0=',MUP
  156. *!      WRITE(IOIMP,*) '||XP||',GNRM2(XP)
  157. *
  158.  10   CONTINUE
  159. * Gestion du CTRL-C
  160.       if (ierr.NE.0) return
  161. *
  162. *     Perform Preconditionned Conjugate Gradient Squared iteration.
  163. *
  164.       ITER = ITER + 1
  165.       IF(IMPR.GT.7) THEN
  166.          WRITE(IOIMP,*) 'ITER=',ITER
  167.       ENDIF
  168. *
  169. C                    rho(i-1)=r[tilde]'r(i-1)
  170.       RHO = GDOT(RTLD,R,IDDOT)
  171. *      IF ((ITER.LE.10).OR.((ITER/100)*100.EQ.ITER)) THEN
  172. *         WRITE(IOIMP,*) 'ITER=',ITER,'  RHO=',RHO
  173. *      ENDIF
  174. C                    if rho(i-1)=0 method fails
  175.       IF (ABS(RHO).LT.RHOTOL) THEN
  176.          RHO=SIGN(XZPREC,RHO)
  177.          IWARN=IWARN+1
  178.       ENDIF
  179. *
  180. *        Compute vector P.
  181. *
  182.       IF (ITER.GT.1) THEN
  183. C                    beta(i-1)=(rho(i-1)/rho(i-2))
  184.          BETA = RHO / RHO1
  185. C                    u(i)=r(i-1)+beta(i-1)*q(i-1)
  186.          CALL GCOPY(R,U)
  187.          CALL GAXPY(BETA,Q,U)
  188. C                    p(i)=u(i)+beta(i-1)*(q(i-1)+beta(i-1)*p(i-1))
  189.          CALL GSCAL(BETA**2,P)
  190.          CALL GAXPY(BETA,Q,P)
  191.          CALL GAXPY(ONE,U,P)
  192.       ELSE
  193. C                    u(1)=r(0)
  194. C                    p(1)=u(1)
  195.          CALL GCOPY(R,U)
  196.          CALL GCOPY(U,P)
  197.       ENDIF
  198. *
  199. *        Compute direction adjusting vector PHAT and scalar ALPHA.
  200. *
  201. C
  202. C        Preconditionner solve
  203. C                    Mp[hat]=p(i)
  204.       CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  205.      $     ,PHAT,P)
  206. C                    vhat=Ap[hat]
  207.       CALL GMOMV(IMVEC,'N',ONE,KMORS,KISA,PHAT,ZERO,VHAT)
  208.       INMV=INMV+1
  209. C                    alpha(i)=rho(i-1)/(r[tilde]'vhat)
  210.       TAU   = GDOT(RTLD,VHAT,IDDOT)
  211.       IF (ABS(TAU).LT.TAUTOL) THEN
  212.          TAU=SIGN(XZPREC,TAU)
  213.          IWARN=IWARN+1
  214.       ENDIF
  215.       ALPHA = RHO / TAU
  216. C                    q(i)=u(i)-alpha(i)vhat
  217.       CALL GCOPY(U,Q)
  218.       CALL GAXPY(-ALPHA,VHAT,Q)
  219. *
  220. *        Compute direction adjusting vector UHAT.
  221. *        PHAT is being used as temporary storage here.
  222. *
  223. C                    M uhat=u(i)+q(i)
  224.       CALL GCOPY(Q,PHAT)
  225.       CALL GAXPY(ONE,U,PHAT)
  226.       CALL PSOLVE('N',KPREC,KMORS,KISA,INVDIA,ILUM,ILUI
  227.      $     ,UHAT,PHAT)
  228. *
  229. *        Compute new solution approximation vector X.
  230. *
  231. *                    x'(i)=x'(i-1)+alpha(i)uhat
  232.       CALL GAXPY(ALPHA,UHAT,XP)
  233. *
  234. *        Compute residual R and check for tolerance
  235. *        using Neumaier's strategy
  236. *
  237. C                    r(i)=b'-Ax'(i)
  238.       CALL GCOPY(BP,R)
  239.       CALL GMOMV(IMVEC,'N',-ONE,KMORS,KISA,XP,ONE,R)
  240.       INMV=INMV+1
  241.       MU=GNRM2(R)
  242.       RESID = MU / BNRM2
  243. CC                    qhat= A uhat
  244. C      CALL GMOMV('N',ONE,KMORS,KISA,UHAT,ZERO,QHAT)
  245. C      INMV=INMV+1
  246. CC                    r(i)=r(i-1)-alpha(i) qhat
  247. C      CALL GAXPY(-ALPHA,QHAT,R)
  248. C      RESID = GNRM2(R) / BNRM2
  249.  
  250.  
  251.       IF(IMPR.GT.6) THEN
  252.          WRITE(IOIMP,*) 'Résidu =',RESID
  253.       ENDIF
  254.       LRES.PROG(ITER+1)=RESID
  255.       LNMV.LECT(ITER+1)=INMV
  256.       IF (RESID.LE.TOL) THEN
  257. *                    x(i)=x(i)+x'(i)
  258.          CALL GAXPY(ONE,XP,X)
  259.          GOTO 30
  260.       ENDIF
  261. *        Neumaier's strategy
  262.       IF (MU.LT.MUP) THEN
  263. *                    x(i)=x(i)+x'(i)
  264.          CALL GAXPY(ONE,XP,X)
  265. *                    x'(i)=0
  266.          CALL GSCAL(ZERO,XP)
  267. *                    b'=r(i)
  268.          CALL GCOPY(R,BP)
  269. *                    mu'=mu
  270.          MUP = MU
  271. *         WRITE(IOIMP,*) 'ITER=',ITER
  272. *         WRITE(IOIMP,*) 'MUP=',MUP
  273.       ENDIF
  274. *
  275. *     Iteration fails.
  276. *
  277. *      IF (ITER.EQ.MAXIT) THEN
  278.       IF (INMV.GE.MAXIT) THEN
  279.          IRET = 1
  280.          IF (IMPR.GT.0) THEN
  281. C            WRITE(IOIMP,*)
  282. C     $           'cgsn.eso : Convergence to tolerance not achieved'
  283. C            WRITE(IOIMP,*) 'INMV=',INMV,
  284. C     $           ' TOL=',TOL,
  285. C     $           ' RESID=',RESID
  286.             IF (IWARN.GT.0) THEN
  287.                WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN
  288.      $              ,' times.'
  289.             ENDIF
  290.          ENDIF
  291.          RETURN
  292.       ENDIF
  293. *
  294. *     Do some more
  295. *
  296.       RHO1 = RHO
  297.       GOTO 10
  298. *
  299. *     Iteration successful
  300. *
  301.  30   CONTINUE
  302.       IF (IWARN.GT.0.AND.IMPR.GT.0) THEN
  303.          WRITE(IOIMP,*) 'CGSN nearly failed ',IWARN,' times.'
  304.       ENDIF
  305.       IRET = 0
  306.       RETURN
  307. *
  308. *     End of CGSN.
  309. *
  310.       END
  311.  
  312.  

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