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.  
  114. -INC PPARAM
  115. -INC CCOPTIO
  116. -INC SMLREEL
  117. INTEGER JG
  118. POINTEUR LRES.MLREEL
  119. -INC SMLENTI
  120. POINTEUR LNMV.MLENTI
  121. POINTEUR MAPREC.MATRIK
  122. POINTEUR KMORS.PMORS
  123. POINTEUR KISA.IZA
  124. POINTEUR KS2B.IZA
  125. POINTEUR INCX.IZA
  126. POINTEUR INVDIA.IZA
  127. POINTEUR ILUM.PMORS
  128. POINTEUR ILUI.IZA
  129. * .. Work Vectors for CG
  130. POINTEUR IR.IZA,IP.IZA,IQ.IZA,IZ.IZA
  131. POINTEUR IRTLD.IZA,IPTLD.IZA,IQTLD.IZA,IZTLD.IZA
  132. * .. Scalar Arguments ..
  133. INTEGER ITER, KPREC, IMPR, IRET
  134. REAL*8 RESID
  135. INTEGER NBVA,NJA,NTT,NTTPRE
  136. * .. Executable Statements ..
  137. *
  138. IRET = 0
  139. *
  140. * On récupère les paramètres
  141. *
  142. SEGACT MATRIK
  143. SEGACT MAPREC
  144. IF (KSYM.NE.0) THEN
  145. IF (IMPR.GT.0) THEN
  146. WRITE(IOIMP,*) 'MATRIK',MATRIK,'non symétrique : ',
  147. $ 'use Bi-CGSTAB instead !'
  148. * ENDIF
  149. * IRET=-1
  150. * GOTO 9999
  151. ENDIF
  152. ENDIF
  153. * Pour le préconditionneur
  154. ILUM =MAPREC.KIDMAT(6)
  155. ILUI =MAPREC.KIDMAT(7)
  156. IDMAT=MAPREC.KIDMAT(1)
  157. SEGACT IDMAT
  158. INVDIA=IDIAG
  159. SEGDES IDMAT
  160. SEGACT KMORS
  161. NTT =KMORS.IA(/1)-1
  162. * NJA =KMORS.JA(/1)
  163. SEGACT KISA
  164. SEGACT KS2B
  165. SEGACT INCX*MOD
  166. C Paramètres des préconditionnements diagonaux et D-ILU
  167. IF (KPREC.NE.0) THEN
  168. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  169. C Est-il compatible avec la matrice ?
  170. SEGACT INVDIA
  171. NTTPRE=INVDIA.A(/1)
  172. IF (NTTPRE.NE.NTT) THEN
  173. WRITE(IOIMP,*) 'La matrice et son préconditionnement'
  174. WRITE(IOIMP,*) 'ne sont pas compatibles...'
  175. WRITE(IOIMP,*) 'NTT=',NTT
  176. WRITE(IOIMP,*) 'NTTPRE=',NTTPRE
  177. IRET=-2
  178. GOTO 9999
  179. ENDIF
  180. C Paramètres des préconditionnements ILU(0), MILU(0), ILUT et ILUT2
  181. ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN
  182. SEGACT ILUM
  183. SEGACT ILUI
  184. NTTPRE=ILUM.IA(/1)
  185. IF (NTTPRE.NE.NTT) THEN
  186. WRITE(IOIMP,*) 'La matrice et son préconditionnement',
  187. $ 'ne sont pas compatibles...'
  188. WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE
  189. IRET=-2
  190. GOTO 9999
  191. ENDIF
  192. ELSE
  193. WRITE(IOIMP,*) 'Erreur de programmation'
  194. GOTO 9999
  195. ENDIF
  196. ENDIF
  197. C
  198. C Initialisation de l'historique de convergence
  199. C
  200. JG=ITER+1
  201. SEGINI LNMV
  202. SEGINI LRES
  203. C
  204. C
  205. C Initialisation des vecteurs de travail pour le gradient conjugué
  206. C
  207. NBVA=NTT
  208. SEGINI IR,IRTLD,IZ,IZTLD,IP,IPTLD
  209. IQ=IZ
  210. IQTLD=IZTLD
  211. C
  212. C Appel du gradient conjugué
  213. C
  214. CALL BICG(KMORS,KISA,KS2B,INCX,
  215. $ KPREC,INVDIA,ILUM,ILUI,
  216. $ LRES,LNMV,
  217. $ IR,IRTLD,IZ,IZTLD,IP,IPTLD,IQ,IQTLD,
  218. $ ITER,INMV,RESID,BRTOL,ICALRS,IDDOT,IMVEC,IMPR,IRET)
  219. C
  220. C Désactivation
  221. C
  222. SEGSUP IR,IRTLD,IZ,IZTLD,IP,IPTLD
  223. JG=ITER+1
  224. SEGADJ LRES
  225. SEGDES LRES
  226. SEGADJ LNMV
  227. SEGDES LNMV
  228. IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN
  229. SEGDES INVDIA
  230. ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN
  231. SEGDES ILUI
  232. SEGDES ILUM
  233. ENDIF
  234. SEGDES INCX
  235. SEGDES KS2B
  236. SEGDES KISA
  237. SEGDES KMORS
  238. SEGDES MAPREC
  239. SEGDES MATRIK
  240. *
  241. * Normal termination
  242. *
  243. RETURN
  244. *
  245. * Format handling
  246. *
  247. *
  248. * Error handling
  249. *
  250. 9999 CONTINUE
  251. WRITE(IOIMP,*) 'An error was detected in prbcg.eso'
  252. RETURN
  253. *
  254. * End of PRBCG.
  255. *
  256. END
  257.  
  258.  
  259.  
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  

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