Télécharger la plaquette Cast3M


fr en

Télécharger prcgsn.eso

Retour à la liste

Numérotation des lignes : 

prcgsn

  1. C PRCGSN    SOURCE    GOUNAND   22/08/25    21:15:10     11434          
  2.       SUBROUTINE PRCGSN(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,
  3.      $     INCX,ITER,INMV,
  4.      $     RESID,KPREC,
  5.      $     BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.       IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM         : PRCGSN
  10. C DESCRIPTION :
  11. C     Préparation de la résolution d'un système linéaire Ax=b
  12. C     par une méthode de gradient conjugué au carré (CGS) avec la
  13. C     stratégie de stabilisation de Neumaier préconditionnée ou non.
  14. C     On stocke deux vecteurs en plus de ceux de CGS normal
  15. C
  16. C LANGAGE     : ESOPE
  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   pour CGS
  29. C
  30. C@TechReport{vorstsleipjen,
  31. C  author =       {GLG Sleijpen and HA. Van der Vorst},
  32. C  title =        {Hybrid Bi-Conjugate Gradient Methods for CFD Problems},
  33. C  institution =  {Universiteit Utrecht},
  34. C  year =         {1995},
  35. C  type =      {Preprint},
  36. C  number =    {902}}
  37. C   pour la stratégie de Neumaier
  38. C***********************************************************************
  39. C***********************************************************************
  40. C VERSION    : v1, 10/02/06, version initiale
  41. C HISTORIQUE : 10/02/06 : Crétation
  42. C HISTORIQUE :
  43. C***********************************************************************
  44. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  45. C en cas de modification de ce sous-programme afin de faciliter
  46. C la maintenance !
  47. C***********************************************************************
  48. *
  49. *     .. Includes et pointeurs associés ..
  50.  
  51. -INC PPARAM
  52. -INC CCOPTIO
  53. -INC SMLREEL
  54.       INTEGER JG
  55.       POINTEUR LRES.MLREEL
  56. -INC SMLENTI
  57.       POINTEUR LNMV.MLENTI
  59.       POINTEUR MAPREC.MATRIK
  60.       POINTEUR KMORS.PMORS
  61.       POINTEUR KISA.IZA
  62.       POINTEUR KS2B.IZA
  63.       POINTEUR INCX.IZA
  64.       POINTEUR INVDIA.IZA
  65.       POINTEUR ILUM.PMORS
  66.       POINTEUR ILUI.IZA
  67. *     .. Parameters
  68. *     This is a breakdown tolerance for BiCGSTAB-type method
  69.       REAL*8     BRTOL
  70. *     .. Work Vectors for CGS..
  71.       POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
  72.       POINTEUR IQ.IZA,IQHAT.IZA,IUHAT.IZA
  73.       POINTEUR IXP.IZA,IBP.IZA
  74. *STAT-INC SMSTAT
  75. *     .. Scalar Arguments ..
  76.       INTEGER ITER, KPREC, IMPR, IRET
  77.       REAL*8  RESID
  78. *     ..
  79. *     Vars reqd for STOPTEST2
  80. *      REAL*8   TOL, BNRM2
  81. *     ..
  82. *     .. External subroutines ..
  83. *      EXTERNAL STOPTEST2
  84.       INTEGER NBVA,NJA,NTT,NTTPRE
  85. *     ..
  86. *     .. Executable Statements ..
  87. *
  88.       IRET = 0
  89. *
  90. * On récupère les paramètres
  91. *
  92.  
  93.       SEGACT MATRIK
  94.       SEGACT MAPREC
  95.       IF (KSYM.EQ.0) THEN
  96.          IF (IMPR.GT.0) THEN
  97.             WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
  98.      $           'Use a Conjugate Gradient instead !'
  99. C         IRET=-2
  100. C         GOTO 9999
  101.          ENDIF
  102.       ENDIF
  103. *     Pour le préconditionneur
  104.       ILUM =MAPREC.KIDMAT(6)
  105.       ILUI =MAPREC.KIDMAT(7)
  106.       IDMAT=MAPREC.KIDMAT(1)
  107.       SEGACT IDMAT
  108.       INVDIA=IDIAG
  109.       SEGDES IDMAT
  110.       SEGACT KMORS
  111.       NTT =KMORS.IA(/1)-1
  112. *      NJA =KMORS.JA(/1)
  113.       SEGACT KISA
  114.       SEGACT KS2B
  115.       SEGACT INCX*MOD
  116.       IF (KPREC.NE.0) THEN
  117. C Paramètres des préconditionnements diagonaux et D-ILU
  118.          IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  119. C     Est-il compatible avec la matrice ?
  120.             SEGACT INVDIA
  121.             NTTPRE=INVDIA.A(/1)
  122.             IF (NTTPRE.NE.NTT) THEN
  123.                WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  124.                WRITE(IOIMP,*) 'ne sont pas compatibles...'
  125.                WRITE(IOIMP,*) 'NTT=',NTT
  126.                WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  127.                IRET=-2
  128.                GOTO 9999
  129.             ENDIF
  130. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2
  131. C ,ilutp, ilutpg, ilutpg2
  132.          ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  133.             SEGACT ILUM
  134.             SEGACT ILUI
  135. *         SEGPRT,ILUM
  136. *         SEGPRT,ILUI
  137.             NTTPRE=ILUM.IA(/1)
  138.             IF (NTTPRE.NE.NTT) THEN
  139.                WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  140.      $              'ne sont pas compatibles...'
  141.                WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  142.                IRET=-2
  143.                GOTO 9999
  144.             ENDIF
  145.          ELSE
  146.             INVDIA=0
  147.             WRITE(IOIMP,*) 'Erreur de programmation'
  148.             GOTO 9999
  149.          ENDIF
  150.       ENDIF
  151. C
  152. C     Initialisation de l'historique de convergence
  153. C
  154.       JG=ITER+1
  155.       SEGINI LNMV
  156.       SEGINI LRES
  157. C
  158. C
  159. C     Initialisation des vecteurs de travail pour BiCGSTAB
  160. C
  161. *STAT      CALL INMSTA(MSTAT,0)
  162.       NBVA=NTT
  163.       SEGINI IR,IRTLD,IP,IPHAT
  164.       SEGINI IQ,IQHAT,IUHAT
  165.       SEGINI IXP,IBP
  166. *STAT      CALL PRMSTA('    PRCGSN : objets temporaires',MSTAT,1)
  167. *STAT      CALL SUMSTA(MSTAT,0)
  168. C
  169. C   Appel de CGSN
  170. C
  171.       CALL CGSN(KMORS,KISA,KS2B,INCX,
  172.      $     KPREC,INVDIA,ILUM,ILUI,
  173.      $     LRES,LNMV,
  174.      $     IR,IRTLD,IP,IPHAT,
  175.      $     IQ,IQHAT,IUHAT,
  176.      $     IXP,IBP,
  177.      $     ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  178. * Gestion du CTRL-C
  179.       if (ierr.NE.0) return
  180. C
  181. C     Désactivation
  182. C
  183.       SEGSUP IXP,IBP
  184.       SEGSUP IR,IRTLD,IP,IPHAT
  185.       SEGSUP IQ,IQHAT,IUHAT
  186.       JG=ITER+1
  187.       SEGADJ LRES
  188.       SEGDES LRES
  189.       SEGADJ LNMV
  190.       SEGDES LNMV
  191.       IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  192.          SEGDES INVDIA
  193.       ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  194.          SEGDES ILUI
  195.          SEGDES ILUM
  196.       ENDIF
  197.       SEGDES INCX
  198.       SEGDES KS2B
  199.       SEGDES KISA
  200.       SEGDES KMORS
  201.       SEGDES MAPREC
  202.       SEGDES MATRIK
  203. C
  204. C     A breakdown has occured in the CGS method
  205. C
  206.       IF (IRET.LT.0) GOTO 9999
  207. *
  208. *     Normal termination
  209. *
  210.       RETURN
  211. *
  212. * Format handling
  213. *
  214. * 1002 FORMAT(10(1X,1PE11.4))
  215. *
  216. * Error handling
  217. *
  218.  9999 CONTINUE
  219.       WRITE(IOIMP,*) 'An error was detected in prcgsn.eso'
  220.       RETURN
  221. *
  222. *     End of PRCGSN
  223. *
  224.       END
  225.  
  226.  

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