Télécharger la plaquette Cast3M


fr en

Télécharger kres13.eso

Retour à la liste

Numérotation des lignes : 

kres13

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

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