Télécharger prcgsn.eso

Retour à la liste

Numérotation des lignes :

  1. C PRCGSN SOURCE PV 16/11/17 22:01:03 9180
  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. -INC CCOPTIO
  51. -INC SMLREEL
  52. INTEGER JG
  53. POINTEUR LRES.MLREEL
  54. -INC SMLENTI
  55. POINTEUR LNMV.MLENTI
  56. POINTEUR MAPREC.MATRIK
  57. POINTEUR KMORS.PMORS
  58. POINTEUR KISA.IZA
  59. POINTEUR KS2B.IZA
  60. POINTEUR INCX.IZA
  61. POINTEUR INVDIA.IZA
  62. POINTEUR ILUM.PMORS
  63. POINTEUR ILUI.IZA
  64. * .. Parameters
  65. * This is a breakdown tolerance for BiCGSTAB-type method
  66. REAL*8 BRTOL
  67. * .. Work Vectors for CGS..
  68. POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA
  69. POINTEUR IQ.IZA,IQHAT.IZA,IUHAT.IZA
  70. POINTEUR IXP.IZA,IBP.IZA
  71. *STAT-INC SMSTAT
  72. * .. Scalar Arguments ..
  73. INTEGER ITER, KPREC, IMPR, IRET
  74. REAL*8 RESID
  75. * ..
  76. * Vars reqd for STOPTEST2
  77. * REAL*8 TOL, BNRM2
  78. * ..
  79. * .. External subroutines ..
  80. * EXTERNAL STOPTEST2
  81. INTEGER NBVA,NJA,NTT,NTTPRE
  82. * ..
  83. * .. Executable Statements ..
  84. *
  85. IRET = 0
  86. *
  87. * On récupère les paramètres
  88. *
  89.  
  90. SEGACT MATRIK
  91. SEGACT MAPREC
  92. IF (KSYM.EQ.0) THEN
  93. IF (IMPR.GT.0) THEN
  94. WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ',
  95. $ 'Use a Conjugate Gradient instead !'
  96. C IRET=-2
  97. C GOTO 9999
  98. ENDIF
  99. ENDIF
  100. * Pour le préconditionneur
  101. ILUM =MAPREC.KIDMAT(6)
  102. ILUI =MAPREC.KIDMAT(7)
  103. IDMAT=MAPREC.KIDMAT(1)
  104. SEGACT IDMAT
  105. INVDIA=IDIAG
  106. SEGDES IDMAT
  107. SEGACT KMORS
  108. NTT =KMORS.IA(/1)-1
  109. * NJA =KMORS.JA(/1)
  110. SEGACT KISA
  111. SEGACT KS2B
  112. SEGACT INCX*MOD
  113. IF (KPREC.NE.0) THEN
  114. C Paramètres des préconditionnements diagonaux et D-ILU
  115. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  116. C Est-il compatible avec la matrice ?
  117. SEGACT INVDIA
  118. NTTPRE=INVDIA.A(/1)
  119. IF (NTTPRE.NE.NTT) THEN
  120. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  121. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  122. WRITE(IOIMP,*) 'NTT=',NTT
  123. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  124. IRET=-2
  125. GOTO 9999
  126. ENDIF
  127. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2
  128. C ,ilutp, ilutpg, ilutpg2
  129. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  130. SEGACT ILUM
  131. SEGACT ILUI
  132. * SEGPRT,ILUM
  133. * SEGPRT,ILUI
  134. NTTPRE=ILUM.IA(/1)
  135. IF (NTTPRE.NE.NTT) THEN
  136. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  137. $ 'ne sont pas compatibles...'
  138. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  139. IRET=-2
  140. GOTO 9999
  141. ENDIF
  142. ELSE
  143. INVDIA=0
  144. WRITE(IOIMP,*) 'Erreur de programmation'
  145. GOTO 9999
  146. ENDIF
  147. ENDIF
  148. C
  149. C Initialisation de l'historique de convergence
  150. C
  151. JG=ITER+1
  152. SEGINI LNMV
  153. SEGINI LRES
  154. C
  155. C
  156. C Initialisation des vecteurs de travail pour BiCGSTAB
  157. C
  158. *STAT CALL INMSTA(MSTAT,0)
  159. NBVA=NTT
  160. SEGINI IR,IRTLD,IP,IPHAT
  161. SEGINI IQ,IQHAT,IUHAT
  162. SEGINI IXP,IBP
  163. *STAT CALL PRMSTA(' PRCGSN : objets temporaires',MSTAT,1)
  164. *STAT CALL SUMSTA(MSTAT,0)
  165. C
  166. C Appel de CGSN
  167. C
  168. CALL CGSN(KMORS,KISA,KS2B,INCX,
  169. $ KPREC,INVDIA,ILUM,ILUI,
  170. $ LRES,LNMV,
  171. $ IR,IRTLD,IP,IPHAT,
  172. $ IQ,IQHAT,IUHAT,
  173. $ IXP,IBP,
  174. $ ITER,INMV,BRTOL,RESID,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  175. C
  176. C Désactivation
  177. C
  178. SEGSUP IXP,IBP
  179. SEGSUP IR,IRTLD,IP,IPHAT
  180. SEGSUP IQ,IQHAT,IUHAT
  181. JG=ITER+1
  182. SEGADJ LRES
  183. SEGDES LRES
  184. SEGADJ LNMV
  185. SEGDES LNMV
  186. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  187. SEGDES INVDIA
  188. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  189. SEGDES ILUI
  190. SEGDES ILUM
  191. ENDIF
  192. SEGDES INCX
  193. SEGDES KS2B
  194. SEGDES KISA
  195. SEGDES KMORS
  196. SEGDES MAPREC
  197. SEGDES MATRIK
  198. C
  199. C A breakdown has occured in the CGS method
  200. C
  201. IF (IRET.LT.0) GOTO 9999
  202. *
  203. * Normal termination
  204. *
  205. RETURN
  206. *
  207. * Format handling
  208. *
  209. * 1002 FORMAT(10(1X,1PE11.4))
  210. *
  211. * Error handling
  212. *
  213. 9999 CONTINUE
  214. WRITE(IOIMP,*) 'An error was detected in prcgsn.eso'
  215. RETURN
  216. *
  217. * End of PRCGSN
  218. *
  219. END
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.  

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