Télécharger la plaquette Cast3M


fr en

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

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