Télécharger la plaquette Cast3M


fr en

Télécharger kres23.eso

Retour à la liste

Numérotation des lignes : 

kres23

  1. C KRES23    SOURCE    GOUNAND   25/04/30    21:15:13     12258          
  2.       SUBROUTINE KRES23(KMINCT,KMORS,KIZA,NNUTOT,IORINC,KTYPI,
  3.      $     IPERM,IPBLOC)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       IMPLICIT INTEGER (I-N)
  6. C***********************************************************************
  7. C     NOM         : KRES23
  8. C     DESCRIPTION : - Reordonnancer par bloc pour le multigrille
  9. C
  10. C                     Ce source est une adaptation de KRES13
  11. C
  12. C     LANGAGE     : ESOPE
  13. C     AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  14. C     mél : gounand@semt2.smts.cea.fr
  15. C***********************************************************************
  16. C     VERSION    : v1, 22/04/2025, version initiale
  17. C     HISTORIQUE : v1, 22/04/2025, création
  18. C     HISTORIQUE :
  19. C     HISTORIQUE :
  20. C***********************************************************************
  21.  
  22. -INC PPARAM
  23. -INC CCOPTIO
  24. -INC CCREEL
  26.       POINTEUR KMINCT.MINC
  27.       POINTEUR KMORS.PMORS
  28.       POINTEUR KIZA.IZA
  29. -INC SMLENTI
  30.       POINTEUR KTYINC.MLENTI
  31.       POINTEUR KTYDDL.MLENTI
  32.       POINTEUR KNODDL.MLENTI
  33.       POINTEUR KRINC.MLENTI
  34.       POINTEUR IWORK.MLENTI
  35.       POINTEUR LAGRAN.MLENTI
  36.       POINTEUR JORDRE.MLENTI,JORTMP.MLENTI
  37.       POINTEUR IPERM.MLENTI,JPETMP.MLENTI
  38.       POINTEUR JPERM.MLENTI,NNUTOT.MLENTI
  39.       POINTEUR IBLOCK.MLENTI,IPBLOC.MLENTI
  40.       POINTEUR ICPR.MLENTI
  41. -INC SMLOBJE
  42.       POINTEUR IORINC.MLOBJE
  43. -INC SMLMOTS
  44.       POINTEUR JORINC.MLMOTS,JORINU.MLMOTS
  45. -INC SMTABLE
  46.       POINTEUR KTIME.MTABLE
  47.       DIMENSION ITTIME(4)
  48.       CHARACTER*(LOCHPO) CHCOMP
  49.       CHARACTER*16 CHARI
  50.       CHARACTER*1 CCOMP
  51.       LOGICAL LTIME,LOGII,LNOD,LBLOCK,LOK,LNUL,LVERIF
  52.       REAL*8 GNRM2
  53. C
  54. C     Executable statements
  55. C
  56.       LVERIF=.FALSE.
  57.       LNOD=(KTYPI.EQ.11)
  58.       SEGACT KMINCT
  59.       NCOMP=KMINCT.LISINC(/2)
  60.       NNOE=KMINCT.MPOS(/1)
  61.       INC=KMINCT.NPOS(NNOE+1)-1
  62.       JG=INC
  63.       SEGINI KTYDDL
  64.       KNODDL=0
  65.       IF (LNOD) SEGINI KNODDL
  66. *
  67.       WRITE(IOIMP,*) 'NCOMP,IORINC= ',NCOMP,IORINC
  68. *
  69.       IF (IORINC.NE.0) THEN
  70.          SEGACT IORINC
  71.          NTYINC=IORINC.LISOBJ(/1)
  72.          JGN=0
  73.          JGM=0
  74.          DO ITYINC=1,NTYINC
  75.             MLMOTS=IORINC.LISOBJ(ITYINC)
  76.             SEGACT MLMOTS
  77.             JGN=MAX(JGN,MOTS(/1))
  78.             JGM=JGM+MOTS(/2)
  79.          ENDDO
  80.          SEGINI,JORINC
  81.          SEGINI,JORINU
  82.          JG=JGM
  83.          IG=0
  84.          SEGINI KTYINC
  85.          DO ITYINC=1,NTYINC
  86.             MLMOTS=IORINC.LISOBJ(ITYINC)
  87.             DO J=1,MOTS(/2)
  88.                IG=IG+1
  89.                JORINC.MOTS(IG)=MOTS(J)
  90.                KTYINC.LECT(IG)=ITYINC
  91.             ENDDO
  92.          ENDDO
  93. *     write(ioimp,*) 'JGN,JGM=',JGN,JGM
  94.          CALL CUNIQ(JORINC.MOTS,JGN,JGM,JORINU.MOTS,JGMU,IMPR,IRET)
  95.          IF (IRET.NE.0) GOTO 9999
  96.          IF (JGM.NE.JGMU) THEN
  97.             WRITE(IOIMP,*) 'IORINC ne doit pas avoir de doublons'
  98.             GOTO 9999
  99.          ENDIF
  100.          SEGSUP JORINU
  101.          IF (JGM.NE.NCOMP) THEN
  102.             WRITE(IOIMP,*) 'IORINC doit referencer toutes'
  103.      $           ,'les inconnues de la matrice'
  104.             GOTO 9999
  105.          ENDIF
  106.          JG=NCOMP
  107.          SEGINI KRINC
  108.          CALL CREPER(JGN,NCOMP,JGM,KMINCT.LISINC,JORINC.MOTS,KRINC.LECT
  109.      $        ,IMPR,IRET)
  110.          IF (IRET.NE.0) GOTO 9999
  111.       ELSE
  112.          NTYINC=NCOMP
  113.       ENDIF
  114.       IF (IIMPI.NE.0) THEN
  115.          WRITE(IOIMP,*) 'NCOMP= ',NCOMP
  116.          WRITE(IOIMP,*) 'KMINCT.LISINC(1..',NCOMP,')= '
  117.          WRITE(IOIMP,*)(KMINCT.LISINC(II),II=1,NCOMP)
  118.          IF (IORINC.NE.0) THEN
  119.             WRITE(IOIMP,*) 'JORINC.MOTS(1..',NCOMP,')= '
  120.             WRITE(IOIMP,*)(JORINC.MOTS(II),II=1,NCOMP)
  121.             WRITE(IOIMP,*) 'KRINC.LECT(1..',NCOMP,')= '
  122.             WRITE(IOIMP,*)(KRINC.LECT(II),II=1,NCOMP)
  123.          ENDIF
  124.       ENDIF
  125. *      segprt,mincpo
  126.       DO ICOMP=1,NCOMP
  127.          DO INOE=1,NNOE
  128.             IPOS=KMINCT.MPOS(INOE,ICOMP)
  129.             IF (IPOS.NE.0) THEN
  130.                JPOS=KMINCT.NPOS(INOE)+IPOS-1
  131.                IF (IORINC.ne.0) then
  132.                   KTYDDL.LECT(JPOS)=KTYINC.LECT(KRINC.LECT(ICOMP))
  133.                ELSE
  134.                   KTYDDL.LECT(JPOS)=ICOMP
  135.                ENDIF
  136.                IF (LNOD) KNODDL.LECT(JPOS)=INOE
  137.             ENDIF
  138.          ENDDO
  139.       ENDDO
  140. *      segprt,ktyddl
  141. *      if (lnod) segprt,knoddl
  142.       segact kmors
  143.       segact kiza
  144. *
  145.       segact nnutot
  146. *      write(ioimp,*) 'Permutation nnutot'
  147. *      write(ioimp,187)  (nnutot.lect(I),I=1,nnutot.lect(/1))
  148. *     1) Trouver le nombre de blocs
  149.       NBLOCK=NTYINC
  150.       ntt=inc
  151.       if (iimpi.ne.0) then
  152.          write(ioimp,*) 'Nb d''inconnues=',ntt
  153.          write(ioimp,*) 'Nb de termes=',kmors.ja(/1)
  154.          write(ioimp,*) 'Nb de blocs detectes nblock=',nblock
  155.       endif
  156.       JG=NBLOCK+1
  157.       SEGINI IBLOCK
  158.       lblock=nblock.gt.1
  159.  187  FORMAT (5X,10I8)
  160.       if (lblock) then
  161.          jg=ntt
  162. *         segini jperm2
  163.          segini iperm
  164.          segini jperm
  165.          do i=1,ntt
  166.             jperm.lect(nnutot.lect(i))=i
  167.          enddo
  168. *         segprt,jperm
  169. *
  170. *     2) Trouver le nombre d'inconnus par blocs
  171.          DO i=1,ntt
  172.             jblock=ktyddl.lect(i)
  173.             iblock.lect(jblock)=iblock.lect(jblock)+1
  174.          enddo
  175.          if (iimpi.ne.0) then
  176.             write(ioimp,*) 'Nb inconnues par blocs'
  177.             write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  178.          endif
  179. *     3) D'où le segment de repérage
  180.          do i=nblock,1,-1
  181.             iblock.lect(i+1)=iblock.lect(i)
  182.          enddo
  183.          if (iimpi.ne.0) then
  184.             write(ioimp,*) 'Segment de repérage tmp'
  185.             write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  186.          endif
  187.          iblock.lect(1)=1
  188.          do i=1,nblock
  189.             iblock.lect(i+1)=iblock.lect(i+1)+iblock.lect(i)
  190.          enddo
  191.          if (iimpi.ne.0) then
  192.             write(ioimp,*) 'Segment de repérage'
  193.             write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  194.          endif
  195. *     4) Construction des permutations
  196.          ntt=kmors.ia(/1)-1
  197.          if (iimpi.ne.0) then
  198.             write(ioimp,*) 'Nb d''inconnues 2 ntt=',ntt
  199.          endif
  200.          if (ktypi.eq.11) then
  201. *     Petite verif sur la taille des blocs
  202.             nddl1=iblock.lect(2)-iblock.lect(1)
  203.             do i=2,nblock
  204.                nddl=iblock.lect(i+1)-iblock.lect(i)
  205.                if (nddl.ne.nddl1) then
  206.                   write(ioimp,*) 'Le bloc ',i,' de taille =',nddl
  207.                   write(ioimp,*) 'na pas la meme taille que le bloc 1 ='
  208.      $                 ,nddl1
  209.                   goto 9999
  210.                endif
  211.             enddo
  212. *     Le premier bloc donne le tableau de correspondance pour les noeuds
  213.             JG=NNOE
  214.             SEGINI ICPR
  215.             ktt=0
  216.             kblock=1
  217.             do jtt=1,ntt
  218.                itt=jperm.lect(jtt)
  219.                jblock=ktyddl.lect(itt)
  220.                if (jblock.eq.kblock) then
  221.                   ktt=ktt+1
  222.                   iperm.lect(itt)=ktt
  223.                   inod=knoddl.lect(itt)
  224. *     jperm2.lect(ktt)=itt
  225.                   ICPR.LECT(INOD)=ktt
  226.                endif
  227.             enddo
  228. *            segprt,iperm
  229. *            segprt,iCPR
  230.            do kblock=2,nblock
  231.                ideb=iblock.lect(kblock)
  232.                do jtt=1,ntt
  233.                   itt=jperm.lect(jtt)
  234.                   jblock=ktyddl.lect(itt)
  235.                   if (jblock.eq.kblock) then
  236.                      inod=knoddl.lect(itt)
  237.                      ktt1=ICPR.LECT(inod)
  238.                      if (ktt1.eq.0) then
  239.                         write(ioimp,*) 'Le noeud ',inod
  240.      $                       ,' present dans le bloc 1 nexiste pas'
  241.      $                       ,' dans le bloc ',kblock
  242.                         goto 9999
  243.                      else
  244.                         ktti=ideb+ktt1-1
  245.                         iperm.lect(itt)=ktti
  246.                      endif
  247.                   endif
  248.                enddo
  249.             enddo
  250. *            segprt,iperm
  251.          else
  252.             ktt=0
  253.             do kblock=1,nblock
  254.                do jtt=1,ntt
  255.                   itt=jperm.lect(jtt)
  256.                   jblock=ktyddl.lect(itt)
  257.                   if (kblock.eq.jblock) then
  258.                      ktt=ktt+1
  259.                      iperm.lect(itt)=ktt
  260. *                  jperm2.lect(ktt)=itt
  261.                   endif
  262.                enddo
  263.             enddo
  264.          endif
  265. *     Attention, on reutilise jperm pour etre la reciproque de iperm et
  266. *     non plus de nnutot
  267.          do i=1,ntt
  268.             jperm.lect(iperm.lect(i))=i
  269.          enddo
  270. *         write(ioimp,*) 'Permutation i'
  271. *         write(ioimp,187)  (iperm.lect(I),I=1,iperm.lect(/1))
  272. *         write(ioimp,*) 'Permutation j'
  273. *         write(ioimp,187)  (jperm.lect(I),I=1,jperm.lect(/1))
  274.       else
  275.          if (ktypi.eq.10.or.ktypi.eq.11) then
  276.             write(ioimp,*)
  277.      $           'The AGMG Stokes and NS solvers need blocks'
  278.             goto 9999
  279.          endif
  280.          iperm=nnutot
  281.          jperm=0
  282.          iblock.lect(1)=1
  283.          iblock.lect(2)=ntt+1
  284.          if (iimpi.ne.0) then
  285.             write(ioimp,*) 'Segment de repérage'
  286.             write(ioimp,187)  (iblock.lect(I),I=1,iblock.lect(/1))
  287.          endif
  288.       endif
  289. *
  290. * Verif pas de terme diagonal nul dans la matrice (sauf AGMG-Stokes KTYPI=10)
  291. *
  292. *     segprt,kmors
  293. *     segprt,kiza
  294.       DO I=1,NTT
  295.          LNUL=.FALSE.
  296.          DO J=KMORS.IA(I),(KMORS.IA(I+1)-1)
  297. *            WRITE(IOIMP,*) 'I,J=',I,J
  298. *            WRITE(IOIMP,*) 'JA(J)=',KMORS.JA(J)
  299.             IF (KMORS.JA(J).EQ.I) THEN
  300. *            WRITE(IOIMP,*) 'KIZA.A(J)=',KIZA.A(J)
  301.                LNUL=(KIZA.A(J).EQ.XZERO)
  302.                IF (LNUL) GOTO 10
  303.             ENDIF
  304.          ENDDO
  305.  10      CONTINUE
  306. *         WRITE(IOIMP,*) 'IFOUND=',IFOUND
  307.          IF (LNUL) THEN
  308.             LOK=.FALSE.
  309. *     Le pivot a le droit d'etre nul s'il s'agit du dernier bloc pour la
  310. *     methode AGMGStokes
  311.             IF (KTYPI.EQ.10) THEN
  312.                IF (KTYDDL.LECT(I).EQ.NBLOCK) THEN
  313.                   LOK=.TRUE.
  314.                ENDIF
  315.             ENDIF
  316.             IF (.NOT.LOK) THEN
  317.                WRITE(IOIMP,*) 'The ',I
  318.      $              ,'th diagonal term of the matrix is nil'
  319.                IF (LNOD) THEN
  320.                   write(ioimp,*) 'Node =',KNODDL.LECT(I)
  321.                ENDIF
  322.                write(ioimp,*) 'Bloc =',KTYDDL.LECT(I)
  323.                IRET=-3
  324.                GOTO 9999
  325.             ENDIF
  326.          ENDIF
  327.       ENDDO
  328. *
  329. * Verifier que les blocs sont de meme taille pour NS
  330. *
  331.       if (ktypi.eq.10.or.ktypi.eq.11.and.LVERIF) then
  332.          IF (.NOT.lblock) THEN
  333.             write(ioimp,*) 'The AGMG Stokes and NS solvers need blocks'
  334.             goto 9999
  335.          ENDIF
  336.          if (ktypi.eq.11) then
  337.             nddl1=iblock.lect(2)-iblock.lect(1)
  338.             do i=2,nblock
  339.                nddl=iblock.lect(i+1)-iblock.lect(i)
  340.                if (nddl.ne.nddl1) then
  341.                   write(ioimp,*) 'Le bloc ',i,' de taille =',nddl
  342.                   write(ioimp,*) 'na pas la meme taille que le bloc 1 ='
  343.      $                 ,nddl1
  344.                   goto 9999
  345.                endif
  346.             enddo
  347. *
  348. *     Verif que les inconnues sont dans le même ordre à l'intérieur de
  349. *     chaque bloc
  350. *
  351.             do i=2,nblock
  352.                nddl=iblock.lect(i+1)-iblock.lect(i)
  353.                ideb=iblock.lect(i)
  354.                do iddl=1,nddl
  355.                   inod1=knoddl.lect(jperm.lect(iddl))
  356.                   inodi=knoddl.lect(jperm.lect(ideb+iddl-1))
  357.                   if (inodi.ne.inod1) then
  358.                      write(ioimp,*) 'Dans le bloc ',i,' le ddl ',iddl
  359.      $                    ,' porte sur le noeud inodi=',inodi
  360.                      write(ioimp,*) 'Dans le bloc ',1,' le ddl ',iddl
  361.      $                    ,' porte sur le noeud inod1=',inod1
  362.                      write(ioimp,*) 'inod1 != inodi est une erreur !'
  363.                      goto 9999
  364.                   endif
  365.                enddo
  366.             enddo
  367.          endif
  368.       endif
  369. *
  370.       IF (IORINC.NE.0) THEN
  371.          SEGSUP KRINC
  372.          SEGSUP KTYINC
  373.          SEGSUP JORINC
  374.       ENDIF
  375.       IF (LNOD) SEGSUP KNODDL
  376.       SEGSUP KTYDDL
  377.       IF (jperm.ne.0) segsup jperm
  378. C
  379.       IPBLOC=IBLOCK
  380. C
  381. C Normal termination
  382. C
  383.       RETURN
  384. C
  385. C Error Handling
  386. C
  387.  9999 CONTINUE
  388.       MOTERR(1:8)='KRES23  '
  389.       CALL ERREUR(1127)
  390.       RETURN
  391. C
  392. C Format handling
  393. C
  394.  2022 FORMAT(10(1X,1PG12.5))
  395. C
  396. C End of subroutine KRES23
  397. C
  398.       END
  399.  
  400.  

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