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    GOUNAND   25/04/30    21:15:28     12258          
  2.       SUBROUTINE PRAGMG(KMORS,KIZA,KS2B,IPBLOC,KTYPI,LRES,
  3.      $     LNMV,INCX,ITER,INMV,
  4.      $     RESID,KPREC,
  5.      $     NRESTS,ICALRS,IDDOT,IMVEC,KTIME,LTIME,LDUMP,ISMOOT,
  6.      $     IMPR,IRET)
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8.       IMPLICIT INTEGER (I-N)
  9. C***********************************************************************
  10. C NOM         : PRAGMG
  11. C DESCRIPTION :
  12. C     Préparation de la résolution d'un système linéaire Ax=b
  13. C     par une méthode Multigrille Algébrique
  14. C
  15. C LANGAGE     : ESOPE
  16. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  17. C               mél : gounand@semt2.smts.cea.fr
  18. C REFERENCE (bibtex-like) :
  19. C***********************************************************************
  20. C***********************************************************************
  21. C VERSION    : v1, 17/06/08, version initiale
  22. C HISTORIQUE : 17/06/08 : Creation
  23. C HISTORIQUE :
  24. C***********************************************************************
  25. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  26. C en cas de modification de ce sous-programme afin de faciliter
  27. C la maintenance !
  28. C***********************************************************************
  29. *
  30. *     .. Includes et pointeurs associés ..
  31. *
  32. -INC PPARAM
  33. -INC CCOPTIO
  34. -INC CCREEL
  35. -INC SMLREEL
  36.       INTEGER JG
  37.       POINTEUR LRES.MLREEL
  38. -INC SMLENTI
  39.       POINTEUR LNMV.MLENTI
  40.       POINTEUR IBLOCK.MLENTI
  41.       POINTEUR IPBLOC.MLENTI
  42.       POINTEUR IPERM.MLENTI
  43.       POINTEUR JPERM.MLENTI
  45.       POINTEUR KMORS.PMORS
  46.       POINTEUR KIZA.IZA,AIZA.IZA
  47.       POINTEUR KS2B.IZA,RES.IZA
  48.       POINTEUR INCX.IZA,INCDX.IZA
  49. -INC SMTABLE
  50.       POINTEUR KTIME.MTABLE
  51.       CHARACTER*8 CHARI
  52.       LOGICAL LTIME,LOGII,LBLOCK,LOK,LDUMP,LSIGN
  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. * On récupère les paramètres
  77. *
  78. *
  79.       segact kmors
  80.       segact kiza
  81.       NTT =KMORS.IA(/1)-1
  82.       NJA =KMORS.JA(/1)
  83.       NBVA=KIZA.A(/1)
  84.       SEGINI,RES=KS2B
  85.       nbva=INCX.A(/1)
  86.  
  87.       SEGINI INCDX
  88.       IBLOCK=IPBLOC
  89.       NBLOCK=IBLOCK.LECT(/1)-1
  90. C
  91. C  Initialisation de l'historique de convergence
  92. C
  93.       JG=ITER+1
  94.       SEGINI LNMV
  95.       SEGINI LRES
  96. *
  97. * Autres paramètres
  98. *
  99.       IJOB=0
  100.       IF (LDUMP) THEN
  101.          IPRINT=-978
  102.       ELSE
  103.          IF (IMPR.GT.2) THEN
  104.             IPRINT=IOIMP
  105.          ELSE
  106.             IPRINT=-1
  107.          ENDIF
  108.       ENDIF
  109.       TOL=RESID
  110.       ICNT=0
  111. *
  112. * AGMG ne fait que du résidu relatif d'où les petits ajustements
  113. * suivants.
  114. *
  115. *     res = b - AX
  116.       CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,RES)
  117.       BNRM2=GNRM2(KS2B)
  118.       RNRM2=GNRM2(RES)
  119.       WRITE(IOIMP,*) 'RNRM2=',RNRM2
  120.       IF (ICALRS.EQ.1) BNRM2=RNRM2
  121.       IF (BNRM2.LT.XPETIT) BNRM2=1.D0
  122.       RESID=RNRM2 / BNRM2
  123.       ICNT=ICNT+1
  124.       LRES.PROG(ICNT)=RESID
  125.       LNMV.LECT(ICNT)=1
  126. *
  127.       IF (RESID.LE.TOL) THEN
  128.          ITER=0
  129.          GOTO 30
  130.       ENDIF
  131. *
  132.       IF (ICALRS.EQ.0) THEN
  133.          XFACT=1.D0/RESID
  134.          TOL=TOL*XFACT
  135.       ENDIF
  136. *
  137. * Changement automatique du signe des lignes de la matrice
  138. * et du second membre si le terme diagonal est négatif.
  139. *
  140.       SEGACT,KMORS
  141.       segact,KIZA
  142. *
  143.       LSIGN=.FALSE.
  144.       DO I=1,NTT
  145.          DO J=KMORS.IA(I),(KMORS.IA(I+1)-1)
  146.             IF (KMORS.JA(J).EQ.I) THEN
  147.                IF (KIZA.A(J).LT.XZERO) THEN
  148.                   LSIGN=.TRUE.
  149.                   IDEB=I
  150.                   GOTO 20
  151.                ELSE
  152.                   LSIGN=.FALSE.
  153.                ENDIF
  154.             ENDIF
  155.          ENDDO
  156.       ENDDO
  157.  20   CONTINUE
  158. *
  159.       IF (LSIGN) THEN
  160.          NSIGN=0
  161.          SEGINI,AIZA=KIZA
  162.          DO I=IDEB,NTT
  163.             IFOUND=0
  164.             DO J=KMORS.IA(I),(KMORS.IA(I+1)-1)
  165. *     WRITE(IOIMP,*) 'I,J=',I,J
  166. *     WRITE(IOIMP,*) 'JA(J)=',KMORS.JA(J)
  167.                IF (KMORS.JA(J).EQ.I) THEN
  168. *     WRITE(IOIMP,*) 'AISA.A(J)=',AISA.A(J)
  169.                   IF (AIZA.A(J).LT.XZERO) THEN
  170.                      IFOUND=-1
  171.                      GOTO 10
  172.                   ENDIF
  173.                ENDIF
  174.             ENDDO
  175.  10         CONTINUE
  176.             IF (IFOUND.EQ.-1) THEN
  177.                NSIGN=NSIGN+1
  178.                RES.A(I)=-1.D0*RES.A(I)
  179.                DO J=KMORS.IA(I),(KMORS.IA(I+1)-1)
  180.                   AIZA.A(J)=-1.D0*AIZA.A(J)
  181.                ENDDO
  182.             ENDIF
  183.          ENDDO
  184.          WRITE(IOIMP,*) NSIGN,' matrix rows changed sign.'
  185.       ELSE
  186.          AIZA=KIZA
  187.       ENDIF
  188. *
  189. * To use the standard AGMG library, you should make sure that it
  190. * is compiled with same compiler and options than Castem's and that
  191. * the sequential example furnished with AGMG works.
  192. * Then it is sufficient to put agmg_seq.o and agmg.o in the
  193. * current directory and use essai for linking.
  194. *
  195. * No interface to the parallel version of agmg for now
  196. *
  197.       WRITE(IOIMP,*) '***********************************'
  198.       WRITE(IOIMP,*) ' '
  199.       WRITE(IOIMP,*) 'Please uncomment the CALL AGMG(...',
  200.      $     ' in the pragmg.eso subroutine if you wish to use'
  201.       WRITE(IOIMP,*) 'the AGMG library by Y. Notay'
  202.       WRITE(IOIMP,*) ' '
  203.       WRITE(IOIMP,*) '***********************************'
  204. * 251 2
  205. *Tentative d'utilisation d'une option non implémentée
  206.       CALL ERREUR(251)
  207.       IRET=-1
  208.       GOTO 9999
  209. *********************************************************************
  210. *
  211. *   AGMG v3.3.7_block call (standard sequential version)
  212. *
  213. *      WRITE(IOIMP,*) '4'
  214. *      CALL ECMORS(KMORS,AIZA,3)
  215.  
  216. c$$$
  217. c$$$      ITR2=ITER-1
  218. c$$$      instance=0
  219. c$$$      NREST=NRESTS
  220. c$$$      IF (KTYPI.EQ.7.OR.KTYPI.EQ.8) then
  221. c$$$         call DAGMG_INTPARAM('smoothtype',ISMOOT,instance)
  222. c$$$         CALL DAGMG(NTT,AIZA.A,KMORS.JA,KMORS.IA,RES.A,INCDX.A,
  223. c$$$     $        IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT)
  224. c$$$         CALL dagmg_status(istat,instance)
  225. c$$$         CALL dagmg_niter(itr2,instance)
  226. c$$$      ELSE
  227. c$$$         if(ktypi.eq.10) then
  228. c$$$            call DAGMG_STOKES_INTPARAM('smoothtype',-1,instance)
  229. c$$$         endif
  230. c$$$         if(ktypi.eq.11) then
  231. c$$$            call DAGMG_STOKES_INTPARAM('smoothtype',-51,instance)
  232. c$$$            call DAGMG_STOKES_INTPARAM('pointcoarsening',1,instance)
  233. c$$$            call DAGMG_STOKES_INTPARAM('Ltransform',1,instance)
  234. c$$$         endif
  235. c$$$         CALL DAGMG_STOKES(NTT,AIZA.A,KMORS.JA,KMORS.IA,RES.A,INCDX.A,
  236. c$$$     $        IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT)
  237. c$$$         CALL dagmg_stokes_status(istat,instance)
  238. c$$$         CALL dagmg_stokes_niter(itr2,instance)
  239. c$$$      ENDIF
  240. c$$$*  X = X_0 + \delta X
  241. c$$$      CALL GAXPY(1.D0,INCDX,INCX)
  242. c$$$*  ON peut l'eviter car deja calcule par AGMG
  243. c$$$** Résidu final
  244. c$$$*      CALL GCOPY(KS2B,RES)
  245. c$$$*      CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,RES)
  246. c$$$*      RNRM2=GNRM2(RES)
  247. c$$$*      RESID=RNRM2/BNRM2
  248. c$$$*
  249. c$$$      IF (istat.NE.0) THEN
  250. c$$$         IRET=1
  251. c$$$      ELSE
  252. c$$$         IRET=0
  253. c$$$      ENDIF
  254. c$$$      ITER=ITR2
  255. c$$$*
  256. c$$$*      write(ioimp,*) 'ITER,nbva=',ITER,nbva
  257. c$$$      JMAX=MIN(ITER+4,nbva)
  258. c$$$      DO J=5,JMAX
  259. c$$$         ICNT=ICNT+1
  260. c$$$*     write(ioimp,*) 'ICNT=',ICNT
  261. c$$$         RESID=RES.A(J)/BNRM2
  262. c$$$         LRES.PROG(ICNT)=RESID
  263. c$$$         LNMV.LECT(ICNT)=LNMV.LECT(ICNT-1)+1
  264. c$$$      ENDDO
  265. c$$$      INMV=1+ITER
  266. *
  267. *   End of AGMG call (standard sequential version)
  268. *
  269. *********************************************************************
  270. *
  271.       IF (LTIME) THEN
  272. *         CHARI='MGAGGREG'
  273. *         CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  274. *     $        'ENTIER  ',IAT,XVALR,CHARR,LOGIR,IRETR)
  275. *         CHARI='MGSOLUTI'
  276. *         CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  277. *     $        'ENTIER  ',IST,XVALR,CHARR,LOGIR,IRETR)
  278.       ENDIF
  279. *
  280. C
  281. C     Désactivation
  282. C
  283.  30   CONTINUE
  284. *      segprt,lres
  285. *      segprt,lnmv
  286.       JG=ICNT
  287.       SEGADJ LRES
  288.       SEGDES LRES
  289.       SEGADJ LNMV
  290.       SEGDES LNMV
  291. *      segprt,lres
  292. *      segprt,lnmv
  293.       SEGSUP INCDX
  294.       SEGDES INCX
  295.       SEGDES KS2B
  296.       SEGSUP RES
  297.       IF (LSIGN) SEGSUP AIZA
  298.       SEGDES KIZA
  299.       SEGDES KMORS
  300. C
  301. C     A breakdown has occured in the CGS method
  302. C
  303.       IF (IRET.LT.0) GOTO 9999
  304. *
  305. *     Normal termination
  306. *
  307.       RETURN
  308. *
  309. * Format handling
  310. *
  311. * 1002 FORMAT(10(1X,1PE11.4))
  312. *
  313. * Error handling
  314. *
  315.  9999 CONTINUE
  316.       WRITE(IOIMP,*) 'An error was detected in pragmg.eso'
  317.       RETURN
  318. *
  319. *     End of PRAGMG
  320. *
  321.       END
  322.  
  323.  

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