Télécharger pragmg.eso

Retour à la liste

Numérotation des lignes :

  1. C PRAGMG SOURCE PV 16/11/17 22:00:57 9180
  2. SUBROUTINE PRAGMG(KMORS,KTYP,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV,
  3. $ INCX,ITER,INMV,
  4. $ RESID,KPREC,
  5. $ NREST,ICALRS,IDDOT,IMVEC,KTIME,LTIME,IMPR,IRET)
  6. IMPLICIT REAL*8 (A-H,O-Z)
  7. IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM : PRAGMG
  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 Multigrille Algébrique
  13. C
  14. C LANGAGE : ESOPE
  15. C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  16. C mél : gounand@semt2.smts.cea.fr
  17. C REFERENCE (bibtex-like) :
  18. C***********************************************************************
  19. C***********************************************************************
  20. C VERSION : v1, 17/06/08, version initiale
  21. C HISTORIQUE : 17/06/08 : Crétation
  22. C HISTORIQUE :
  23. C***********************************************************************
  24. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  25. C en cas de modification de ce sous-programme afin de faciliter
  26. C la maintenance !
  27. C***********************************************************************
  28. *
  29. * .. Includes et pointeurs associés ..
  30. -INC CCOPTIO
  31. -INC CCREEL
  32. -INC SMLREEL
  33. INTEGER JG
  34. POINTEUR LRES.MLREEL
  35. -INC SMLENTI
  36. POINTEUR LNMV.MLENTI
  37. POINTEUR ATYP.MLENTI
  38. POINTEUR KTYP.MLENTI
  39. POINTEUR MAPREC.MATRIK
  40. POINTEUR KMORS.PMORS
  41. POINTEUR KISA.IZA
  42. POINTEUR AMORS.PMORS
  43. POINTEUR AISA.IZA
  44. POINTEUR KS2B.IZA,RES.IZA
  45. POINTEUR INCX.IZA
  46. POINTEUR INVDIA.IZA
  47. POINTEUR ILUM.PMORS
  48. POINTEUR ILUI.IZA
  49. -INC SMTABLE
  50. POINTEUR KTIME.MTABLE
  51. CHARACTER*8 CHARI
  52. LOGICAL LTIME,LOGII
  53. * .. Parameters
  54. * This is a breakdown tolerance for BiCGSTAB-type method
  55. REAL*8 BRTOL
  56. *STAT-INC SMSTAT
  57. * .. Scalar Arguments ..
  58. INTEGER ITER, KPREC, IMPR, IRET
  59. INTEGER ICNT
  60. REAL*8 RESID,TOL,BNRM2,RNRM2,XFACT
  61. * ..
  62. * Vars reqd for STOPTEST2
  63. * REAL*8 TOL, BNRM2
  64. * ..
  65. * .. External subroutines ..
  66. * EXTERNAL STOPTEST2
  67. INTEGER NBVA,NJA,NTT,NTTPRE
  68. REAL*8 GNRM2
  69. EXTERNAL GNRM2
  70. * ..
  71. * .. Executable Statements ..
  72. *
  73. * WRITE(IOIMP,*) 'Debut de pragmg'
  74. IRET = 0
  75.  
  76. IVALI=0
  77. XVALI=0.D0
  78. LOGII=.FALSE.
  79. IRETI=0
  80. IAT =0
  81. XVALR=0.D0
  82. IRETR=0
  83. IST =0
  84. *
  85. * On récupère les paramètres
  86. *
  87. *
  88. SEGACT MATRIK
  89. SEGINI,AMORS=KMORS
  90. NTT =AMORS.IA(/1)-1
  91. * NJA =KMORS.JA(/1)
  92. SEGINI,ATYP=KTYP
  93. SEGINI,AISA=KISA
  94. SEGACT KS2B
  95. SEGINI,RES=KS2B
  96. SEGACT INCX*MOD
  97. C
  98. C Initialisation de l'historique de convergence
  99. C
  100. JG=ITER+1
  101. SEGINI LNMV
  102. SEGINI LRES
  103. *
  104. * Autres paramètres
  105. *
  106. IJOB=0
  107. IF (IMPR.GT.2) THEN
  108. IPRINT=IOIMP
  109. ELSE
  110. IPRINT=-1
  111. ENDIF
  112. TOL=RESID
  113. ICNT=0
  114. *
  115. * AGMG ne fait que du résidu relatif d'où les petits ajustements
  116. * suivants.
  117. *
  118. * res = b - AX
  119. CALL GMOMV(IMVEC,'N',-1.D0,AMORS,AISA,INCX,1.D0,RES)
  120. BNRM2=GNRM2(KS2B)
  121. RNRM2=GNRM2(RES)
  122. IF (ICALRS.EQ.1) BNRM2=RNRM2
  123. IF (BNRM2.LT.XPETIT) BNRM2=1.D0
  124. RESID=RNRM2 / BNRM2
  125. LRES.PROG(ICNT+1)=RESID
  126. LNMV.LECT(ICNT+1)=1
  127. *
  128. IF (RESID.LE.TOL) THEN
  129. ITER=0
  130. GOTO 30
  131. ENDIF
  132. *
  133. IF (ICALRS.EQ.0) THEN
  134. XFACT=1.D0/RESID
  135. TOL=TOL*XFACT
  136. ENDIF
  137. *
  138. * Changement automatique du signe des lignes de la matrice
  139. * et du second membre si le terme diagonal est négatif.
  140. *
  141. * CALL ECMORS(AMORS,AISA,3)
  142. DO I=1,NTT
  143. IFOUND=0
  144. DO J=AMORS.IA(I),(AMORS.IA(I+1)-1)
  145. * WRITE(IOIMP,*) 'J=',J
  146. * WRITE(IOIMP,*) 'JA(J)=',AMORS.JA(J)
  147. IF (AMORS.JA(J).EQ.I) THEN
  148. IF (AISA.A(J).GT.XZERO) THEN
  149. IFOUND=1
  150. ELSEIF (AISA.A(J).LT.XZERO) THEN
  151. IFOUND=-1
  152. ENDIF
  153. GOTO 10
  154. ENDIF
  155. ENDDO
  156. 10 CONTINUE
  157. * WRITE(IOIMP,*) 'IFOUND=',IFOUND
  158. IF (IFOUND.EQ.0) THEN
  159. WRITE(IOIMP,*) 'The ',I
  160. $ ,'th diagonal term of the matrix is nil'
  161. IRET=-3
  162. GOTO 9999
  163. ELSEIF (IFOUND.EQ.-1) THEN
  164. RES.A(I)=-1.D0*RES.A(I)
  165. DO J=AMORS.IA(I),(AMORS.IA(I+1)-1)
  166. AISA.A(J)=-1.D0*AISA.A(J)
  167. ENDDO
  168. ENDIF
  169. ENDDO
  170. * WRITE(IOIMP,*) '4'
  171. * CALL ECMORS(AMORS,AISA,3)
  172. *
  173. * To use the standard AGMG library, you should make sure that it
  174. * is compiled with same compiler and options than Castem's and that
  175. * the sequential example furnished with AGMG works.
  176. * Then it is sufficient to put agmg_seq.o and agmg.o in the
  177. * current directory and use essai for linking.
  178. *
  179. * No interface to the parallel version of agmg for now
  180. *
  181. WRITE(IOIMP,*) '***********************************'
  182. WRITE(IOIMP,*) ' '
  183. WRITE(IOIMP,*) 'Please uncomment the CALL AGMG(...',
  184. $ ' in the pragmg.eso subroutine if you wish to use'
  185. WRITE(IOIMP,*) 'the AGMG library by Y. Notay'
  186. WRITE(IOIMP,*) ' '
  187. WRITE(IOIMP,*) '***********************************'
  188. * 251 2
  189. *Tentative d'utilisation d'une option non implémentée
  190. CALL ERREUR(251)
  191. IRET=-1
  192. GOTO 9999
  193. *********************************************************************
  194. *
  195. * AGMG v1.1 call (standard sequential version)
  196. *
  197. C ITR2=ITER-1
  198. C CALL AGMG(NTT,AISA.A,AMORS.JA,AMORS.IA,RES.A,
  199. C $ IJOB,IPRINT,NREST,ITR2,TOL)
  200. C* X = X_0 + \delta X
  201. C CALL GAXPY(1.D0,RES,INCX)
  202. C* Résidu final
  203. C CALL GCOPY(KS2B,RES)
  204. C CALL GMOMV(IMVEC,'N',-1.D0,AMORS,AISA,INCX,1.D0,RES)
  205. C RNRM2=GNRM2(RES)
  206. C RESID=RNRM2/BNRM2
  207. C*
  208. C IF (ITR2.LT.0) THEN
  209. C IRET=1
  210. C ITER=-ITR2
  211. C ELSE
  212. C IRET=0
  213. C ITER=ITR2
  214. C ENDIF
  215. C*
  216. C ICNT=ICNT+1
  217. C LRES.PROG(ICNT+1)=RESID
  218. C LNMV.LECT(ICNT+1)=1+ITER
  219. C INMV=1+ITER
  220. *
  221. * End of AGMG call (standard sequential version)
  222. *
  223. *********************************************************************
  224. *********************************************************************
  225. *
  226. * AGMG call (sequential version slightly modified
  227. * for better error handling and back transmission of residual
  228. * norm history)
  229. *
  230. C ITR2=ITER-1
  231. C* CALL ECMORS(AMORS,AISA,3)
  232. C* WRITE(IOIMP,*) 'Appel de agmg'
  233. C CALL AGMG(NTT,AISA.A,AMORS.JA,AMORS.IA,ATYP.LECT,RES.A,
  234. C $ IJOB,IPRINT,NREST,ITR2,TOL,LTIME,IAT,IST,LRES.PROG,IDDOT)
  235. C* X = X_0 + \delta X
  236. C CALL GAXPY(1.D0,RES,INCX)
  237. C*
  238. C IF (ITR2.LT.0) THEN
  239. C IRET=1
  240. C ITER=-ITR2
  241. C ELSE
  242. C IRET=0
  243. C ITER=ITR2
  244. C ENDIF
  245. C IF (ICALRS.EQ.0) THEN
  246. C DO I=1,ITER
  247. C LRES.PROG(I+1)=LRES.PROG(I+1)/XFACT
  248. C ENDDO
  249. C ENDIF
  250. C DO I=1,ITER
  251. C LNMV.LECT(I+1)=I+1
  252. C ENDDO
  253. C INMV=1+ITER
  254. C RESID=LRES.PROG(ITER+1)
  255. C ICNT=ITER
  256. *
  257. * End of AGMG call (sequential version slightly modified
  258. * for better error handling and back transmission of residual
  259. * norm history)
  260. *
  261. *********************************************************************
  262. *
  263. IF (LTIME) THEN
  264. CHARI='MGAGGREG'
  265. CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI,
  266. $ 'ENTIER ',IAT,XVALR,CHARR,LOGIR,IRETR)
  267. CHARI='MGSOLUTI'
  268. CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI,
  269. $ 'ENTIER ',IST,XVALR,CHARR,LOGIR,IRETR)
  270. ENDIF
  271. *
  272. C
  273. C Désactivation
  274. C
  275. 30 CONTINUE
  276. JG=ICNT+1
  277. SEGADJ LRES
  278. SEGDES LRES
  279. SEGADJ LNMV
  280. SEGDES LNMV
  281. SEGDES INCX
  282. SEGDES KS2B
  283. SEGSUP RES
  284. SEGSUP AISA
  285. SEGSUP ATYP
  286. SEGSUP AMORS
  287. SEGDES KISA
  288. SEGDES KTYP
  289. SEGDES KMORS
  290. SEGDES MATRIK
  291. C
  292. C A breakdown has occured in the CGS method
  293. C
  294. IF (IRET.LT.0) GOTO 9999
  295. *
  296. * Normal termination
  297. *
  298. RETURN
  299. *
  300. * Format handling
  301. *
  302. * 1002 FORMAT(10(1X,1PE11.4))
  303. *
  304. * Error handling
  305. *
  306. 9999 CONTINUE
  307. WRITE(IOIMP,*) 'An error was detected in pragmg.eso'
  308. RETURN
  309. *
  310. * End of PRAGMG
  311. *
  312. END
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  

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