Télécharger pragmg.eso

Retour à la liste

Numérotation des lignes :

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

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