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/03/24    21:15:08     12216          
  2.       SUBROUTINE PRAGMG(KMORS,KTYP,KNOD,KISA,KS2B,MATRIK,MAPREC,LRES,
  3.      $     LNMV,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 KTYP.MLENTI
  40.       POINTEUR KNOD.MLENTI
  41.       POINTEUR IBLOCK.MLENTI
  42.       POINTEUR IPERM.MLENTI
  43.       POINTEUR JPERM.MLENTI
  45.       POINTEUR MAPREC.MATRIK
  46.       POINTEUR KMORS.PMORS
  47.       POINTEUR KISA.IZA
  48.       POINTEUR AMORS.PMORS
  49.       POINTEUR AISA.IZA
  50.       POINTEUR KS2B.IZA,RES.IZA
  51.       POINTEUR INCX.IZA,INCDX.IZA
  52.       POINTEUR INVDIA.IZA
  53.       POINTEUR ILUM.PMORS
  54.       POINTEUR ILUI.IZA
  55. -INC SMTABLE
  56.       POINTEUR KTIME.MTABLE
  57.       CHARACTER*8 CHARI
  58.       LOGICAL LTIME,LOGII,LBLOCK
  59. *     .. Parameters
  60. *     This is a breakdown tolerance for BiCGSTAB-type method
  61.       REAL*8     BRTOL
  62. *STAT-INC SMSTAT
  63. *     .. Scalar Arguments ..
  64.       INTEGER ITER, KPREC, IMPR, IRET
  65.       INTEGER ICNT
  66.       REAL*8  RESID,TOL,BNRM2,RNRM2,XFACT
  67. *     ..
  68. *     Vars reqd for STOPTEST2
  69. *      REAL*8   TOL, BNRM2
  70. *     ..
  71. *     .. External subroutines ..
  72. *      EXTERNAL STOPTEST2
  73.       INTEGER NBVA,NJA,NTT,NTTPRE
  74.       REAL*8 GNRM2
  75.       EXTERNAL GNRM2
  76. *     ..
  77. *     .. Executable Statements ..
  78. *
  79. *      WRITE(IOIMP,*) 'Debut de pragmg'
  80.       IRET = 0
  81.  
  82.       IVALI=0
  83.       XVALI=0.D0
  84.       LOGII=.FALSE.
  85.       IRETI=0
  86.       IAT  =0
  87.       XVALR=0.D0
  88.       IRETR=0
  89.       IST  =0
  90. *
  91. * On récupère les paramètres
  92. *
  93. *
  94.       SEGACT MATRIK
  95.       SEGACT KMORS
  96.       NTT =KMORS.IA(/1)-1
  97.       NJA =KMORS.JA(/1)
  98.       SEGACT KISA
  99.       NBVA=KISA.A(/1)
  100.       SEGACT KS2B
  101.       SEGINI,RES=KS2B
  102.       SEGACT INCX*MOD
  103.       nbva=INCX.A(/1)
  104.       SEGINI INCDX
  105. C
  106. C  Initialisation de l'historique de convergence
  107. C
  108.       JG=ITER+1
  109.       SEGINI LNMV
  110.       SEGINI LRES
  111. *
  112. * Autres paramètres
  113. *
  114.       IJOB=0
  115.       IF (IMPR.GT.2) THEN
  116.          IPRINT=IOIMP
  117.       ELSE
  118.          IPRINT=-1
  119.       ENDIF
  120.       TOL=RESID
  121.       ICNT=0
  122. *
  123. * AGMG ne fait que du résidu relatif d'où les petits ajustements
  124. * suivants.
  125. *
  126. *     res = b - AX
  127.       CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KISA,INCX,1.D0,RES)
  128.       BNRM2=GNRM2(KS2B)
  129.       RNRM2=GNRM2(RES)
  130.       IF (ICALRS.EQ.1) BNRM2=RNRM2
  131.       IF (BNRM2.LT.XPETIT) BNRM2=1.D0
  132.       RESID=RNRM2 / BNRM2
  133.       LRES.PROG(ICNT+1)=RESID
  134.       LNMV.LECT(ICNT+1)=1
  135. *
  136.       IF (RESID.LE.TOL) THEN
  137.          ITER=0
  138.          GOTO 30
  139.       ENDIF
  140. *
  141.       IF (ICALRS.EQ.0) THEN
  142.          XFACT=1.D0/RESID
  143.          TOL=TOL*XFACT
  144.       ENDIF
  145. *
  146. * To use the standard AGMG library, you should make sure that it
  147. * is compiled with same compiler and options than Castem's and that
  148. * the sequential example furnished with AGMG works.
  149. * Then it is sufficient to put agmg_seq.o and agmg.o in the
  150. * current directory and use essai for linking.
  151. *
  152. * No interface to the parallel version of agmg for now
  153. *
  154.       WRITE(IOIMP,*) '***********************************'
  155.       WRITE(IOIMP,*) ' '
  156.       WRITE(IOIMP,*) 'Please uncomment the CALL AGMG(...',
  157.      $     ' in the pragmg.eso subroutine if you wish to use'
  158.       WRITE(IOIMP,*) 'the AGMG library by Y. Notay'
  159.       WRITE(IOIMP,*) ' '
  160.       WRITE(IOIMP,*) '***********************************'
  161. * 251 2
  162. *Tentative d'utilisation d'une option non implémentée
  163.       CALL ERREUR(251)
  164.       IRET=-1
  165.       GOTO 9999
  166. *********************************************************************
  167. *
  168. *   AGMG v3.3.7_block call (standard sequential version)
  169. *
  170. *     1) Trouver le nombre de blocs
  171. c$$$      NBLOCK=0
  172. c$$$      segact ktyp
  173. c$$$      ntt=ktyp.lect(/1)
  174. c$$$      if (nrest.lt.0) then
  175. c$$$         segact knod
  176. c$$$         ntt2=knod.lect(/1)
  177. c$$$         if (ntt2.ne.ntt) THEN
  178. c$$$            call erreur(5)
  179. c$$$            goto 9999
  180. c$$$         endif
  181. c$$$      endif
  182. c$$$      jg=ntt
  183. c$$$      segini iperm
  184. c$$$      segini jperm
  185. c$$$      nja=kisa.a(/1)
  186. c$$$      DO i=1,ntt
  187. c$$$         jblock=ktyp.lect(i)
  188. c$$$         if (jblock.gt.nblock) nblock=jblock
  189. c$$$      enddo
  190. c$$$      if (impr.gt.2) then
  191. c$$$         write(ioimp,*) 'Nb d''inconnues=',ntt
  192. c$$$         write(ioimp,*) 'Nb de termes=',kmors.ja(/1)
  193. c$$$         write(ioimp,*) 'Nb de termes 2=',nja
  194. c$$$         write(ioimp,*) 'Nb de blocs detectes nblock=',nblock
  195. c$$$      endif
  196. c$$$      JG=NBLOCK+1
  197. c$$$      SEGINI IBLOCK
  198. c$$$      lblock=nblock.gt.1
  199. c$$$ 187  FORMAT (5X,10I8)
  200. c$$$      if (lblock) then
  201. c$$$*
  202. c$$$*     2) Trouver le nombre d'inconnus par blocs
  203. c$$$         DO i=1,ktyp.lect(/1)
  204. c$$$            jblock=ktyp.lect(i)
  205. c$$$            iblock.lect(jblock)=iblock.lect(jblock)+1
  206. c$$$         enddo
  207. c$$$         if (impr.gt.2) then
  208. c$$$            write(ioimp,*) 'Nb inconnues par blocs'
  209. c$$$            write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  210. c$$$         endif
  211. c$$$*     3) D'où le segment de repérage
  212. c$$$         do i=nblock,1,-1
  213. c$$$            iblock.lect(i+1)=iblock.lect(i)
  214. c$$$         enddo
  215. c$$$         if (impr.gt.2) then
  216. c$$$            write(ioimp,*) 'Segment de repérage tmp'
  217. c$$$            write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  218. c$$$         endif
  219. c$$$         iblock.lect(1)=1
  220. c$$$         do i=1,nblock
  221. c$$$            iblock.lect(i+1)=iblock.lect(i+1)+iblock.lect(i)
  222. c$$$         enddo
  223. c$$$         if (impr.gt.2) then
  224. c$$$            write(ioimp,*) 'Segment de repérage'
  225. c$$$            write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  226. c$$$         endif
  227. c$$$*     4) Construction des permutations
  228. c$$$         ntt=kmors.ia(/1)-1
  229. c$$$         if (impr.gt.2) then
  230. c$$$            write(ioimp,*) 'Nb d''inconnues 2 ntt=',ntt
  231. c$$$         endif
  232. c$$$         ktt=1
  233. c$$$         do kblock=1,nblock
  234. c$$$            do itt=1,ntt
  235. c$$$               jblock=ktyp.lect(itt)
  236. c$$$               if (kblock.eq.jblock) then
  237. c$$$                  iperm.lect(itt)=ktt
  238. c$$$                  jperm.lect(ktt)=itt
  239. c$$$                  ktt=ktt+1
  240. c$$$               endif
  241. c$$$            enddo
  242. c$$$         enddo
  243. c$$$*         write(ioimp,*) 'Permutation i'
  244. c$$$*         write(ioimp,187)  (iperm.lect(I),I=1,iperm.lect(/1))
  245. c$$$*         write(ioimp,*) 'Permutation j'
  246. c$$$*         write(ioimp,187)  (jperm.lect(I),I=1,jperm.lect(/1))
  247. c$$$*
  248. c$$$*     5) Permutation de la matrice et du second membre
  249. c$$$*
  250. c$$$         SEGINI,AMORS=KMORS
  251. c$$$         SEGINI,AISA=KISA
  252. c$$$*         SEGINI,AMORS
  253. c$$$*         SEGINI,AISA
  254. c$$$*         write(ioimp,*) 'Avant permutation'
  255. c$$$*         CALL ECMORS(AMORS,AISA,3)
  256. c$$$         job=1
  257. c$$$         CALL DPERM(ntt,kisa.a,kmors.ja,kmors.ia,
  258. c$$$     $        aisa.a,amors.ja,amors.ia,iperm.lect,jperm.lect,job)
  259. c$$$         CALL DVPERM(ntt,res.a,iperm.lect)
  260. c$$$         if (nrest.lt.0) CALL IVPERM(ntt,knod.lect,iperm.lect)
  261. c$$$*         write(ioimp,*) 'Apres permutation'
  262. c$$$*         CALL ECMORS(AMORS,AISA,3)
  263. c$$$      else
  264. c$$$         iblock.lect(1)=1
  265. c$$$         iblock.lect(2)=ntt+1
  266. c$$$         if (impr.gt.2) then
  267. c$$$            write(ioimp,*) 'Segment de repérage'
  268. c$$$            write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  269. c$$$         endif
  270. c$$$* Pas de permutation
  271. c$$$         SEGINI,AMORS=KMORS
  272. c$$$         SEGINI,AISA=KISA
  273. c$$$      endif
  274. c$$$*
  275. c$$$* Changement automatique du signe des lignes de la matrice
  276. c$$$* et du second membre si le terme diagonal est négatif.
  277. c$$$*
  278. c$$$*      SEGACT,AMORS*mod
  279. c$$$*      segact,AISA*mod
  280. c$$$      DO I=1,NTT
  281. c$$$         IFOUND=0
  282. c$$$         DO J=AMORS.IA(I),(AMORS.IA(I+1)-1)
  283. c$$$*            WRITE(IOIMP,*) 'I,J=',I,J
  284. c$$$*            WRITE(IOIMP,*) 'JA(J)=',AMORS.JA(J)
  285. c$$$            IF (AMORS.JA(J).EQ.I) THEN
  286. c$$$*               WRITE(IOIMP,*) 'AISA.A(J)=',AISA.A(J)
  287. c$$$               IF (AISA.A(J).GT.XZERO) THEN
  288. c$$$                  IFOUND=1
  289. c$$$               ELSEIF (AISA.A(J).LT.XZERO) THEN
  290. c$$$                  IFOUND=-1
  291. c$$$               ENDIF
  292. c$$$               GOTO 10
  293. c$$$            ENDIF
  294. c$$$         ENDDO
  295. c$$$ 10      CONTINUE
  296. c$$$*         WRITE(IOIMP,*) 'IFOUND=',IFOUND
  297. c$$$         IF (IFOUND.EQ.0) THEN
  298. c$$$            WRITE(IOIMP,*) 'The ',I
  299. c$$$     $           ,'th diagonal term of the matrix is nil'
  300. c$$$            IRET=-3
  301. c$$$            GOTO 9999
  302. c$$$         ELSEIF (IFOUND.EQ.-1) THEN
  303. c$$$            RES.A(I)=-1.D0*RES.A(I)
  304. c$$$            DO J=AMORS.IA(I),(AMORS.IA(I+1)-1)
  305. c$$$               AISA.A(J)=-1.D0*AISA.A(J)
  306. c$$$            ENDDO
  307. c$$$         ENDIF
  308. c$$$      ENDDO
  309. c$$$*      WRITE(IOIMP,*) '4'
  310. c$$$*      CALL ECMORS(AMORS,AISA,3)
  311. c$$$
  312. c$$$
  313. c$$$      ITR2=ITER-1
  314. c$$$      instance=0
  315. c$$$      IF (NREST.GE.0) then
  316. c$$$*      IPRINT=-978
  317. c$$$         CALL DAGMG(NTT,AISA.A,AMORS.JA,AMORS.IA,RES.A,INCDX.A,
  318. c$$$     $        IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT)
  319. c$$$      ELSE
  320. c$$$         IF (.NOT.lblock) THEN
  321. c$$$            write(ioimp,*) 'The AGMG Stokes solver needs blocks'
  322. c$$$            moterr(1:8)='pragmg'
  323. c$$$            call erreur(1127)
  324. c$$$            return
  325. c$$$         ENDIF
  326. c$$$         NRESTA=ABS(NREST)
  327. c$$$         if (nresta.GT.1) then
  328. c$$$            nddl1=iblock.lect(2)-iblock.lect(1)
  329. c$$$            do i=2,nblock
  330. c$$$               nddl=iblock.lect(i+1)-iblock.lect(i)
  331. c$$$               if (nddl.ne.nddl1) then
  332. c$$$                  write(ioimp,*) 'Le bloc ',i,' de taille =',nddl
  333. c$$$                  write(ioimp,*) 'na pas la meme taille que le bloc 1 ='
  334. c$$$     $                 ,nddl1
  335. c$$$                  moterr(1:8)='pragmg'
  336. c$$$                  call erreur(1127)
  337. c$$$                  return
  338. c$$$               endif
  339. c$$$            enddo
  340. c$$$*
  341. c$$$*     Vérif que les inconnues sont dans le même ordre à l'intérieur de
  342. c$$$*     chaque bloc
  343. c$$$*
  344. c$$$            do i=2,nblock
  345. c$$$               nddl=iblock.lect(i+1)-iblock.lect(i)
  346. c$$$               ideb=iblock.lect(i)
  347. c$$$               do iddl=1,nddl
  348. c$$$                  inod1=knod.lect(iddl)
  349. c$$$                  inodi=knod.lect(ideb+iddl-1)
  350. c$$$                  if (inodi.ne.inod1) then
  351. c$$$                     write(ioimp,*) 'Dans le bloc ',i,' le ddl ',iddl
  352. c$$$     $                    ,' porte sur le noeud inodi=',inodi
  353. c$$$                     write(ioimp,*) 'Dans le bloc ',1,' le ddl ',iddl
  354. c$$$     $                    ,' porte sur le noeud inod1=',inod1
  355. c$$$                     write(ioimp,*) 'inod1 != inodi est une erreur !'
  356. c$$$                     moterr(1:8)='pragmg'
  357. c$$$                     call erreur(1127)
  358. c$$$                     return
  359. c$$$                  endif
  360. c$$$               enddo
  361. c$$$            enddo
  362. c$$$
  363. c$$$            call DAGMG_STOKES_INTPARAM('smoothtype',-51,instance)
  364. c$$$            call DAGMG_STOKES_INTPARAM('pointcoarsening',1,instance)
  365. c$$$            call DAGMG_STOKES_INTPARAM('Ltransform',1,instance)
  366. c$$$         endif
  367. c$$$         CALL DAGMG_STOKES(NTT,AISA.A,AMORS.JA,AMORS.IA,RES.A,INCDX.A,
  368. c$$$     $        IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT)
  369. c$$$
  370. c$$$      ENDIF
  371. c$$$      CALL dagmg_status(istat,instance)
  372. c$$$      CALL dagmg_niter(itr2,instance)
  373. c$$$
  374. c$$$      if (lblock) then
  375. c$$$         CALL DVPERM(ntt,incdx.a,jperm.lect)
  376. c$$$      endif
  377. c$$$      SEGSUP IBLOCK
  378. c$$$      segsup iperm
  379. c$$$      segsup jperm
  380. c$$$
  381. c$$$
  382. c$$$
  383. c$$$*  X = X_0 + \delta X
  384. c$$$      CALL GAXPY(1.D0,INCDX,INCX)
  385. c$$$* Résidu final
  386. c$$$      CALL GCOPY(KS2B,RES)
  387. c$$$      CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KISA,INCX,1.D0,RES)
  388. c$$$      RNRM2=GNRM2(RES)
  389. c$$$      RESID=RNRM2/BNRM2
  390. c$$$*
  391. c$$$      IF (istat.NE.0) THEN
  392. c$$$         IRET=1
  393. c$$$      ELSE
  394. c$$$         IRET=0
  395. c$$$      ENDIF
  396. c$$$      ITER=ITR2
  397. c$$$*
  398. c$$$      ICNT=ICNT+1
  399. c$$$      LRES.PROG(ICNT+1)=RESID
  400. c$$$      LNMV.LECT(ICNT+1)=1+ITER
  401. c$$$      INMV=1+ITER
  402. *
  403. *   End of AGMG call (standard sequential version)
  404. *
  405. *********************************************************************
  406. *
  407.       IF (LTIME) THEN
  408.          CHARI='MGAGGREG'
  409.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  410.      $        'ENTIER  ',IAT,XVALR,CHARR,LOGIR,IRETR)
  411.          CHARI='MGSOLUTI'
  412.          CALL ECCTAB(KTIME,'MOT     ',IVALI,XVALI,CHARI,LOGII,IRETI,
  413.      $        'ENTIER  ',IST,XVALR,CHARR,LOGIR,IRETR)
  414.       ENDIF
  415. *
  416. C
  417. C     Désactivation
  418. C
  419.  30   CONTINUE
  420.       JG=ICNT+1
  421.       SEGADJ LRES
  422.       SEGDES LRES
  423.       SEGADJ LNMV
  424.       SEGDES LNMV
  425.       SEGSUP INCDX
  426.       SEGDES INCX
  427.       SEGDES KS2B
  428.       SEGSUP RES
  429.       SEGSUP AISA
  430. *      SEGSUP ATYP
  431.       SEGSUP AMORS
  432.       SEGDES KISA
  433.       SEGDES KTYP
  434.       if (nrest.lt.0) SEGDES KNOD
  435.       SEGDES KMORS
  436.       SEGDES MATRIK
  437. C
  438. C     A breakdown has occured in the CGS method
  439. C
  440.       IF (IRET.LT.0) GOTO 9999
  441. *
  442. *     Normal termination
  443. *
  444.       RETURN
  445. *
  446. * Format handling
  447. *
  448. * 1002 FORMAT(10(1X,1PE11.4))
  449. *
  450. * Error handling
  451. *
  452.  9999 CONTINUE
  453.       WRITE(IOIMP,*) 'An error was detected in pragmg.eso'
  454.       RETURN
  455. *
  456. *     End of PRAGMG
  457. *
  458.       END
  459.  
  460.  

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