Télécharger prbcg.eso

Retour à la liste

Numérotation des lignes :

  1. C PRBCG SOURCE PV 16/11/17 22:01:01 9180
  2. SUBROUTINE PRBCG(KMORS,KISA,KS2B,MATRIK,MAPREC,
  3. $ LRES,LNMV,INCX,ITER,INMV,RESID,KPREC,BRTOL,ICALRS,IDDOT,
  4. $ IMVEC,
  5. $ IMPR,IRET)
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8 (A-H,O-Z)
  8. C***********************************************************************
  9. C NOM : PRBCG
  10. C DESCRIPTION :
  11. C Préparation de la résolution d'un système linéaire Ax=b
  12. C par la méthode du gradient conjugué (préconditionnée ou non)
  13. C
  14. C A doit etre symétrique.
  15. C
  16. C Les 4 préconditionnements disponibles pour cette méthode
  17. C sont :
  18. C * Jacobi (diagonal)
  19. C * D-ILU (Diagonal Incomplete LU factorization)
  20. C * ILU(0) (Incomplete LU factorization
  21. C of level zero ie Choleski incomplet)
  22. C * MILU(0) (Modified ILU(0)) relaxé
  23. C 5 : préconditionnement ILUT (dual truncation)
  24. C 6 : préconditionnement ILUT2 (une variante du
  25. C précédent qui remplit mieux la mémoire et
  26. C fonctionne mieux quelquefois)
  27. C
  28. C Ces préconditionneurs sont supposés déjà construits (ie
  29. C factorisés), voir les subroutines MEIDIA, PRDILU, PRILU0,
  30. C PRMILU et PRILUT.
  31. C
  32. C Ce sous-programme est en fait une interface aux sous-programmes :
  33. C cg, cgd, cgdilu, cgilu0
  34. C qui sont en Fortran presque pur (pour raison de rapidité)
  35. C et effectuent la résolution du système linéaire
  36. C proprement dite.
  37. C
  38. C LANGAGE : ESOPE
  39. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  40. C mél : gounand@semt2.smts.cea.fr
  41. C REFERENCE (bibtex-like) :
  42. C @BOOK{templates,
  43. C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato,
  44. C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine,
  45. C H. Van der Vorst},
  46. C TITLE={Templates for the Solution of Linear Systems :
  47. C Building Blocks for Iterative Methods},
  48. C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} }
  49. C -> URL : http://www.netlib.org/templates/Templates.html
  50. C***********************************************************************
  51. C APPELES (Calcul) : CG, CGD, CGDILU, CGILU0
  52. C***********************************************************************
  53. C ENTREES : MATRIK, MAPREC, KPREC, IMPR
  54. C ENTREES/SORTIES : INCX, ITER, RESID
  55. C SORTIES : LRES, IRET
  56. C CODE RETOUR (IRET) : 0 si ok
  57. C <0 si problème
  58. C >0 si l'algorithme appelé n'a pas convergé
  59. C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
  60. C on pioche dedans les informations nécessaires
  61. C (différents pointeurs, nb. de ddl...)
  62. C MAPREC : pointeur sur segment MATRIK de l'include SMMATRIK
  63. C on va l'utiliser comme préconditionneur
  64. C KPREC : type de préconditionnement à effectuer
  65. C 0 : pas de préconditionnement
  66. C 1 : préconditionnement Jacobi (diagonal)
  67. C 2 : préconditionnement D-ILU
  68. C 3 : préconditionnement ILU(0) (Choleski incomplet)
  69. C 4 : préconditionnement MILU(0)
  70. C (Choleski incomplet modifié)
  71. C 5 : préconditionnement ILUT (dual truncation)
  72. C 6 : préconditionnement ILUT2 (une variante du
  73. C précédent qui remplit mieux la mémoire et
  74. C fonctionne mieux quelquefois)
  75. C IMPR : niveau d'impression (0..4)
  76. C INCX : pointeur sur segment IZA de l'include SMMATRIK
  77. C C'est le vecteur des inconnues.
  78. C En entrée, contient l'estimation x(0).
  79. C En sortie, la solution trouvée (si la méthode
  80. C a convergé).
  81. C ITER : type INTEGER.
  82. C En entrée, contient le nombre maximum
  83. C d'itérations à effectuer.
  84. C En sortie, contient le nombre d'itérations
  85. C réellement effectuées.
  86. C RESID : type REAL*8.
  87. C En entrée, la valeur maximale du résidu à
  88. C convergence de l'algorithme mesurée comme suit :
  89. C norme[L2](b - A*x) / norme[L2]( b )
  90. C En sortie, la valeur effective de cette mesure.
  91. C LRES : pointeur sur segment MLREEL (include SMLREEL)
  92. C contient l'historique de la valeur de RESID en
  93. C fonction du nombre d'itérations effectuées.
  94. C
  95. C***********************************************************************
  96. C VERSION : v1, 02/04/98, version initiale
  97. C HISTORIQUE : v1, 02/04/98, création
  98. C HISTORIQUE : 09/02/99 : ajout préconditionneur <> MATRIK
  99. C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau cgilu0
  100. C (factorisations incomplètes stockées au format MSR (Modified Sparse
  101. C Row) (voir la doc de Sparskit version 2+)
  102. C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html
  103. C HISTORIQUE : 22/03/00 : ajout des préconditionneurs ILUT
  104. C HISTORIQUE : 10/02/06 : Nettoyage
  105. C HISTORIQUE :
  106. C***********************************************************************
  107. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  108. C en cas de modification de ce sous-programme afin de faciliter
  109. C la maintenance !
  110. C***********************************************************************
  111. *
  112. * .. Includes et pointeurs associés ..
  113. -INC CCOPTIO
  114. -INC SMLREEL
  115. INTEGER JG
  116. POINTEUR LRES.MLREEL
  117. -INC SMLENTI
  118. POINTEUR LNMV.MLENTI
  119. POINTEUR MAPREC.MATRIK
  120. POINTEUR KMORS.PMORS
  121. POINTEUR KISA.IZA
  122. POINTEUR KS2B.IZA
  123. POINTEUR INCX.IZA
  124. POINTEUR INVDIA.IZA
  125. POINTEUR ILUM.PMORS
  126. POINTEUR ILUI.IZA
  127. * .. Work Vectors for CG
  128. POINTEUR IR.IZA,IP.IZA,IQ.IZA,IZ.IZA
  129. POINTEUR IRTLD.IZA,IPTLD.IZA,IQTLD.IZA,IZTLD.IZA
  130. * .. Scalar Arguments ..
  131. INTEGER ITER, KPREC, IMPR, IRET
  132. REAL*8 RESID
  133. INTEGER NBVA,NJA,NTT,NTTPRE
  134. * .. Executable Statements ..
  135. *
  136. IRET = 0
  137. *
  138. * On récupère les paramètres
  139. *
  140. SEGACT MATRIK
  141. SEGACT MAPREC
  142. IF (KSYM.NE.0) THEN
  143. IF (IMPR.GT.0) THEN
  144. WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
  145. $ 'use Bi-CGSTAB instead !'
  146. * ENDIF
  147. * IRET=-1
  148. * GOTO 9999
  149. ENDIF
  150. ENDIF
  151. * Pour le préconditionneur
  152. ILUM =MAPREC.KIDMAT(6)
  153. ILUI =MAPREC.KIDMAT(7)
  154. IDMAT=MAPREC.KIDMAT(1)
  155. SEGACT IDMAT
  156. INVDIA=IDIAG
  157. SEGDES IDMAT
  158. SEGACT KMORS
  159. NTT =KMORS.IA(/1)-1
  160. * NJA =KMORS.JA(/1)
  161. SEGACT KISA
  162. SEGACT KS2B
  163. SEGACT INCX*MOD
  164. C Paramètres des préconditionnements diagonaux et D-ILU
  165. IF (KPREC.NE.0) THEN
  166. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  167. C Est-il compatible avec la matrice ?
  168. SEGACT INVDIA
  169. NTTPRE=INVDIA.A(/1)
  170. IF (NTTPRE.NE.NTT) THEN
  171. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  172. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  173. WRITE(IOIMP,*) 'NTT=',NTT
  174. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  175. IRET=-2
  176. GOTO 9999
  177. ENDIF
  178. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
  179. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  180. SEGACT ILUM
  181. SEGACT ILUI
  182. NTTPRE=ILUM.IA(/1)
  183. IF (NTTPRE.NE.NTT) THEN
  184. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  185. $ 'ne sont pas compatibles...'
  186. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  187. IRET=-2
  188. GOTO 9999
  189. ENDIF
  190. ELSE
  191. WRITE(IOIMP,*) 'Erreur de programmation'
  192. GOTO 9999
  193. ENDIF
  194. ENDIF
  195. C
  196. C Initialisation de l'historique de convergence
  197. C
  198. JG=ITER+1
  199. SEGINI LNMV
  200. SEGINI LRES
  201. C
  202. C
  203. C Initialisation des vecteurs de travail pour le gradient conjugué
  204. C
  205. NBVA=NTT
  206. SEGINI IR,IRTLD,IZ,IZTLD,IP,IPTLD
  207. IQ=IZ
  208. IQTLD=IZTLD
  209. C
  210. C Appel du gradient conjugué
  211. C
  212. CALL BICG(KMORS,KISA,KS2B,INCX,
  213. $ KPREC,INVDIA,ILUM,ILUI,
  214. $ LRES,LNMV,
  215. $ IR,IRTLD,IZ,IZTLD,IP,IPTLD,IQ,IQTLD,
  216. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  217. C
  218. C Désactivation
  219. C
  220. SEGSUP IR,IRTLD,IZ,IZTLD,IP,IPTLD
  221. JG=ITER+1
  222. SEGADJ LRES
  223. SEGDES LRES
  224. SEGADJ LNMV
  225. SEGDES LNMV
  226. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  227. SEGDES INVDIA
  228. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  229. SEGDES ILUI
  230. SEGDES ILUM
  231. ENDIF
  232. SEGDES INCX
  233. SEGDES KS2B
  234. SEGDES KISA
  235. SEGDES KMORS
  236. SEGDES MAPREC
  237. SEGDES MATRIK
  238. *
  239. * Normal termination
  240. *
  241. RETURN
  242. *
  243. * Format handling
  244. *
  245. *
  246. * Error handling
  247. *
  248. 9999 CONTINUE
  249. WRITE(IOIMP,*) 'An error was detected in prbcg.eso'
  250. RETURN
  251. *
  252. * End of PRBCG.
  253. *
  254. END
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  

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