Télécharger la plaquette Cast3M


fr en

Télécharger ldmts.eso

Retour à la liste

Numérotation des lignes : 

ldmts

  1. C LDMTS     SOURCE    CB215821  26/06/25    21:15:17     12581          
  2.       SUBROUTINE LDMTS(MMATRX,PREC,ISTAB,NBNNMA,NLIGRA,XMATRI)
  3. *
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6. C  TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY
  7.       PARAMETER (LPCRAY=10000)
  8.       INTEGER OOOVAL,OOOLEN
  9.       dimension ittime(4)
  10. C
  11. C  ****  MISE SOUS FORME A=LDM DE LA MATRICE MMATRX
  12. C
  13. -INC PPARAM
  14. -INC CCOPTIO
  15. -INC CCREEL
  16. -INC SMMATRI
  17. -INC SMRIGID
  18. -INC CCASSIS
  19. -INC CCHOLE
  20. -INC SMCOORD
  21. C
  22.       POINTEUR LILIGN.MILIGN
  23.       POINTEUR LLIGL.LLIGN,LLIGM.LLIGN
  24.       POINTEUR LIGL.LIGN,LIGM.LIGN,LIG5.LIGN,LIG3.LIGN
  25.       SEGMENT KIVPO(IIMAX)
  26.       SEGMENT KIVLO(IIMAX)
  27.       segment immt(nblig)
  28.       segment xreser
  29.        integer ireser(ivstrm)
  30.        real yreser(nvstrm)
  31.       endsegment
  32.       SEGMENT ILR
  33.         integer ILIGR(NNOE)
  34.       ENDSEGMENT
  35.       POINTEUR ILIGNS.ILR
  36.       external chole4i
  37.       external shole3i
  38.       SAVE IPASV
  39.       DATA IPASV/0/
  40. *  ngmpet dit si on tient en memoire (false) ou si on deborde (true)
  41.       logical ngmpet
  42.       logical lsgdes,pasfait
  43. C
  44.       SEGDES,MCOORD
  45.       maitre=0
  46.       xreser=0
  47.       pasfait=.true.
  48.       lsgdes=.false.
  49. *  faire attention a respecter l'ordre des segdes par la suite
  50.       call ooomru(1)
  51.       if (istab.ne.0.and.istab.ne.1) call erreur(5)
  52.       condmax=0.d0
  53.       condmin=xgrand*xzprec
  54.       ngmpet=.false.
  55.       nvallp=nvall
  56.       call timespv(ittime,oothrd)
  57.       kcour=(ittime(1)+ittime(2))/10
  58.       kcouri=kcour
  59.       kcour=0
  60.       nbopit=0
  61.       nvaor=0
  62.       nbthro=nbthrs
  63.       ithrd=0
  64.       if (nbthro.gt.1) then
  65.        ithrd=1
  66.        call threadii
  67.        call oooprl(1)
  68.       endif
  69.       nbthr=nbthro
  70.       do ith=1,nbthr
  71.        nbop(ith)=0
  72.       enddo
  73.       stmult=1d-5
  74.  
  75.       nvstrm=0
  76.       ivstrm=oooval(1,4)*2
  77. C     write(6,*) 'nvstrm initial',nvstrm
  78.       MMATRI=MMATRX
  79.       SEGACT,MMATRI*MOD
  80.       PRCHLV=PREC
  81.       MILIGN=IILIGN
  82.       LILIGN=IILIGS
  83.       SEGACT,MILIGN*MOD,LILIGN*MOD
  84. C
  85.       INO=MILIGN.ILIGN(/1)
  86.       NNOE=INO
  87.       SEGINI,ILIGNS
  88. C
  89.       MDIAG=IDIAG
  90.       SEGACT,MDIAG*MOD
  91.       NBLIG=INO
  92.       segini immt
  93.       precc=prec
  94.       INC=DIAG(/1)
  95. C
  96. C     Cas du super-element
  97.       nbnnmc=inc+1
  98.       if (xmatri.ne.0) then
  99.         nelrig=1
  100.         nligrp=nligra
  101.         nligrd=nligra
  102.         rigrel=0
  103.         segini xmatri
  104.         nbnnmc=nbnnma
  105.       endif
  106.       matric=xmatri
  107.       INCC=INC
  108.       MIMIK=IIMIK
  109.       MINCPO=IINCPO
  110.       SEGACT,MINCPO,MIMIK
  111.       IPLUMI=IMIK(/2)*2 +4
  112.       IL2=0
  113.       IIMAX=IJMAX+IPLUMI
  114.       SEGINI KIVPO,KIVLO
  115.       INEG=0
  116.       NBLAG=0
  117.       NENSLX=0
  118.       NVSTOC=0
  119.       NVSTOR=0
  120.       diagmax=XPETIT/XZPREC
  121.       diagmin=xgrand
  122.       do i=1,diag(/1)
  123.        if (ittr(i).eq.0) diagmax=max(diagmax,abs(diag(i)))
  124.        if (ittr(i).eq.0.and.abs(diag(i)).gt.xpetit/xzprec)
  125.      >   diagmin=min(diagmin,abs(diag(i)))
  126.       enddo
  127.       if (diagmax.le.xpetit/xzprec) then
  128.          do i=1,diag(/1)
  129.            diagmax=max(diagmax,abs(diag(i)))
  130.            if (abs(diag(i)).gt.xpetit/xzprec)
  131.      >       diagmin=min(diagmin,abs(diag(i)))
  132.          enddo
  133.       endif
  134.       diagmin=min(diagmin,diagmax)
  135. ***    write (6,*) ' ldmts diagmin diagmax ',diagmin,diagmax,diag(/1)
  136. C
  137. C  **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  138. C  **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  139. C
  140. C  ****  LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  141. C
  142. C SP  indicateurs pour impression message "stabilisation RESO..."
  143.       isr=0
  144.       isrl=0
  145.  
  146.   1   CONTINUE
  147.       IL1=IL2+1
  148.       IMIN=IL1
  149. C
  150. C     Reserver de la place ou mettre les lignes superieures dans le cas debordement
  151.       if (ngmpet) then
  152.         if (xreser.eq.0) then
  153. CC      write(6,*)'segini xreser',nvstrm
  154.         segini xreser
  155.         endif
  156.       endif
  157. C
  158.       DO 2 I=IL1,INO
  159. C
  160.       LIGL=0
  161.       LIGM=0
  162.       LGLIG=0
  163. C
  164.       LLIGM=LILIGN.ILIGN(I)
  165.       LLIGL=MILIGN.ILIGN(I)
  166.       SEGACT /ERR=32/ LLIGM
  167.       SEGACT /ERR=32/ LLIGL
  168.       GOTO 31
  169.   32  CONTINUE
  170.       IF (LSGDES) GOTO 3
  171. **    write(6,*) 'desactivation-1 ',1,il1-1
  172.       lsgdes=.true.
  173.       do it=il1-1,1,-1
  174.         lig1=ILIGNS.ILIGR(it)
  175.         segdes lig1
  176.         lig1=LILIGN.ilign(it)
  177.         segdes lig1
  178.         lig1=MILIGN.ilign(it)
  179.         segdes lig1
  180.       enddo
  181.       SEGACT /ERR=3/LLIGM
  182.       SEGACT /ERR=3/LLIGL
  183.   31  CONTINUE
  184. ** essai pv
  185.       NVMAX=MAX(LLIGM.NJMAX,LLIGL.NJMAX)
  186.       if (ngmpet) then
  187.         if (nvmax.gt.2*njmaxp.and.i-il1+1.gt.nbthrs/2) then
  188. **        write(6,*) 'fin de bloc il1  i nvmax',il1,i,nvmax
  189.           goto 3
  190.         endif
  191.       endif
  192.       NJMAXP=NVMAX
  193. ** fin essai
  194. C
  195.       NA1  =LLIGM.IMMMM(/1)
  196.       NVAL1=LLIGM.IMMMM(NA1)-LLIGM.LDEB(1)
  197.       LDEB1=LLIGM.LDEB(1)
  198.       NA2  =LLIGL.IMMMM(/1)
  199.       NVAL2=LLIGL.IMMMM(NA2)-LLIGL.LDEB(1)
  200.       LDEB2=LLIGL.LDEB(1)
  201. C
  202. C     Il peut arriver que les segments ligne et colonne
  203. C     d'un meme noeud ne demarrent pas de la meme colonne/ligne.
  204. C     Decalage pour que ce ne soit plus le cas
  205.       LDEBM=MIN(LDEB1,LDEB2)
  206.       IF      (LDEB1-LDEBM.NE.0) THEN
  207.         SEGACT,LLIGM*MOD
  208.         LLIGM.NJMAX=LLIGL.NJMAX
  209.         DO IMV=1,NA1
  210.           LLIGM.LDEB(IMV)=LDEBM
  211.         ENDDO
  212.       ELSE IF (LDEB2-LDEBM.NE.0) THEN
  213.         SEGACT,LLIGL*MOD
  214.         LLIGL.NJMAX=LLIGM.NJMAX
  215.         DO IMV=1,NA2
  216.           LLIGL.LDEB(IMV)=LDEBM
  217.         ENDDO
  218.       ENDIF
  219. C
  220.       NVALL=0
  221.       NVALT=MAX(NVAL1,NVAL2)
  222.       if (ngmpet.and.nvalt.lt.nvallp/3) then
  223.         write(6,*) 'passage en lent'
  224.         ngmpet=.false.
  225.       endif
  226.       NVALLP=NVALT
  227.       if (ngmpet) then
  228. C       nvall=(NVMAX*nvstor)/nvstoc
  229.         NVALL=NVMAX
  230.         if (nvall.gt.NVSTRM) then
  231. C       write(*,*)'redimensionnement de xreser',nvstrm,(4*NVALL),il1,il2
  232.           nvstrm=4*nvall
  233.           segsup,xreser
  234.           segini,xreser
  235.         endif
  236.       endif
  237. C
  238. C     Le masque est uniquement present dans LIGM
  239.       NA=NA1
  240.       NBPAR=NA+1
  241.       LMASQ=MASQA(NVALT+MASDIM)
  242.       SEGINI /ERR=33/ LIGM
  243. C
  244.       NA=NA2
  245.       NBPAR=NA+1
  246.       LMASQ=0
  247.       SEGINI /ERR=33/ LIGL
  248. C
  249.       GOTO 34
  250.   33  CONTINUE
  251.       IF (LSGDES) GOTO 3
  252. **    write(6,*) 'desactivation-2 ',1,il1-1
  253.       lsgdes=.true.
  254.       do it=il1-1,1,-1
  255.         lig1=ILIGNS.ILIGR(it)
  256.         segdes lig1
  257.         lig1=LILIGN.ilign(it)
  258.         segdes lig1
  259.         lig1=MILIGN.ilign(it)
  260.         segdes lig1
  261.       enddo
  262.       IF (LIGM.EQ.0) THEN
  263.         NA=NA1
  264.         NBPAR=NA+1
  265.         LMASQ=MASQA(NVAL1+MASDIM)
  266.         SEGINI /ERR=3/ LIGM
  267.       ENDIF
  268.       IF (LIGL.EQ.0) THEN
  269.         NA=NA2
  270.         NBPAR=NA+1
  271.         LMASQ=0
  272.         SEGINI /ERR=3/ LIGL
  273.       ENDIF
  274.   34  CONTINUE
  275. C
  276. C     Travail sur la partie superieure
  277.       MILIGN=IILIGS
  278.       LLIGN =LLIGM
  279.       LIGN  =LIGM
  280.       ILIGNS.ILIGR(I)=LIGM
  281. C
  282.   35  CONTINUE
  283. C
  284. C     Recuperer la longueur du segment
  285.       NA    =LIGN.IMMM(/1)
  286.       LGLIG =LGLIG  + NA*INT(REAL(NJMAX/NA)**(4.D0/3.D0))
  287.       NVSTOC=NVSTOC + LLIGN.NJMAX
  288.       NVAOR =NVAOR  + LLIGN.XXVA(/1)
  289. C
  290. C  ****  DECOMPACTAGE
  291. C
  292.       IPA=1
  293.       DO 121 JPA=1,NA
  294.         KPA =LLIGN.IPPO(JPA+1)-LLIGN.IPPO(JPA)
  295.         IPP =LLIGN.IPPO(JPA)
  296.         LPA =LLIGN.LDEB(JPA)
  297.         LPA1=LPA-IPA
  298.         LIGN.IMMM(JPA) =LPA
  299.         LIGN.IVPO(JPA) =IPA
  300.         LIGN.IPPVV(JPA)=IPA-1
  301.         DO 122 MPA=1,KPA
  302.           IPLA=LLIGN.LINC(MPA+IPP)-LPA1
  303.           ICC =IPLA-IPA+1
  304.           ICC1=ICC-MASQD(ICC)+1
  305.           MICC=MASQA(ICC)
  306.           IMSQ=LIGM.IMASQ(MICC)
  307.           IF (IMSQ.LE.0) THEN
  308.             MSQH=icc1
  309.             MSQB=icc1
  310.           ELSE
  311.             MSQH=MAX(MASQH(IMSQ),icc1)
  312.             MSQB=MIN(MASQB(IMSQ),icc1)
  313.           ENDIF
  314.           LIGM.IMASQ(MICC)=MASQV(MSQB,MSQH)
  315. 122     CONTINUE
  316.         IPA=IPA+LLIGN.IMMMM(JPA)-LPA + 1
  317.         IF(IMIN.GT.MILIGN.IPNO(LPA)) IMIN=MILIGN.IPNO(LPA)
  318. 121   CONTINUE
  319.       LIGN.IPPVV(NA+1)=IPA-1
  320. C
  321.       IF (NA.GT.0) THEN
  322.         LIGN.IPREL=LLIGN.IMMMM(1)
  323.         LIGN.IDERL=LLIGN.IMMMM(NA)
  324.         MILIGN.LCARA(2,I)=LIGN.IPREL
  325.         MILIGN.LCARA(3,I)=LIGN.IDERL
  326.       ENDIF
  327. C
  328. C     Travail sur la partie inferieure
  329.       IF (MILIGN.EQ.IILIGS) THEN
  330.         MILIGN=IILIGN
  331.         LLIGN =LLIGL
  332.         LIGN  =LIGL
  333.         GOTO 35
  334.       ENDIF
  335. C
  336. C     Indexation de imasq
  337.       LMASQ=LIGM.IMASQ(/1)
  338.       ipln=lmasq
  339.       idec=1
  340.       do 123 ipl=lmasq,1,-1
  341.         if (LIGM.imasq(ipl).gt.0) then
  342.           ipln=ipl-1
  343.           idec=MASQB(LIGM.IMASQ(IPL))
  344.         else
  345.           LIGM.IMASQ(IPL)=-(IPLN*MASDIM+IDEC)
  346.         endif
  347.  123  continue
  348. CCC   write (6,*) ' imasq ',lmasq
  349. CCC   write (6,*) (LIGM.imasq(ipl),ipl=1,lmasq)
  350. C
  351. C*    write (6,*) 'longueur ligne ',nvall
  352. C nb de ligne multiple du nb de threads
  353. C blocage ligne lecture-ecriture pour minimiser le cache
  354. C on note si on est au minimum de lignes
  355.       nbthro=min(nbthrs,lglig/800+1)
  356.       nbthro2=min(nbthrs,lglig/10000+1)
  357.       IF (.not.ngmpet) THEN
  358.         if (i+1-il1.ge.nbthro2) then
  359.           nbthro2=min(nbthrs,i+1-il1)
  360.         endif
  361.         if (i+1-il1.ge.nbthro) then
  362.           nbthro=min(nbthrs,i+1-il1)
  363.           il2=i
  364. **        write(6,*) ' nbthro il1 il2 ',nbthro,il1,il2
  365.           GOTO 4
  366.         endif
  367.       ENDIF
  368.  
  369.     2 CONTINUE
  370.       IL2=INO
  371.       GOTO 4
  372.     3 CONTINUE
  373.       IL2=I-1
  374. **    write(6,*) 'desactivation-4 ',llign
  375.       SEGDES LLIGM,LLIGL
  376.       IF (LIGL.NE.0) SEGSUP,LIGL
  377.       IF (LIGM.NE.0) SEGSUP,LIGM
  378.     4 CONTINUE
  379.       nbthro=min(nbthrs,nbthro)
  380.       nbthro2=min(nbthrs,nbthro2)
  381.       nbthr=nbthro
  382.       if(xreser.ne.0) then
  383. CC      write(6,*)'segsup xreser 1',il1,il2
  384.         segsup xreser
  385.         xreser=0
  386.       endif
  387. C
  388.       IF(IL2.GE.IL1) GOTO 40
  389. C
  390. C ****  MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  391. C
  392.       CALL ERREUR(48)
  393. **    call ooodmp(0)
  394.       if (ithrd.eq.1) then
  395.        call threadis
  396.        call oooprl(0)
  397.       endif
  398.       call ooomru(0)
  399.       RETURN
  400.    40 CONTINUE
  401. C
  402. C     Travail sur la partie superieure
  403.       MILIGN=IILIGS
  404. C
  405.    39 CONTINUE
  406. C
  407.       IM=INC
  408.       DO 352 IH=IL2,IL1,-1
  409.         LIGN=ILIGNS.ILIGR(IH)
  410.         IL=INC
  411.         DO 354 JH=1,LIGN.IMMM(/1)
  412.           IM=MIN(IM,LIGN.IMMM(JH))
  413.           IL=MIN(IL,LIGN.IMMM(JH))
  414.  354    CONTINUE
  415.         LIGN.IML=IL
  416.         MILIGN.LCARA(1,IH)=LIGN.IML
  417.         IMMT(IH)=MILIGN.IPNO(IM)
  418.   352 CONTINUE
  419. C
  420. C     Travail sur la partie inferieure
  421.       IF (MILIGN.EQ.IILIGS) THEN
  422.         MILIGN=IILIGN
  423.         GOTO 39
  424.       ENDIF
  425. C
  426.       IMAX=IL1-1
  427.       IF (LSGDES) THEN
  428.         DO IT=IMIN,IL1-1
  429.           LIG1=ILIGNS.ILIGR(IT)
  430.           SEGACT /ERR=333/ LIG1
  431.         ENDDO
  432.   333   CONTINUE
  433.         IMAX=IT-1
  434.       ENDIF
  435. C
  436.       nbthrp=nbthr
  437.       nbthr=min(nbthrp,il2-il1+1)
  438.       nbthr=nbthro2
  439.       IPER=IMIN
  440.       IDER=IMAX
  441.       IL1P=IL1
  442.       DO 13 I=IL1,IL2
  443. C
  444.         MILIGN=IILIGS
  445.   336   CONTINUE
  446. C
  447. C       Travail sur la partie inferieure
  448.         LILIGN=IILIGN
  449.         DO ITH=1,NBTHR-1
  450.           CALL THREADID(ITH,SHOLE3I)
  451.         ENDDO
  452.         CALL SHOLE3I(NBTHR)
  453.         DO ITH=NBTHR-1,1,-1
  454.           CALL THREADIF(ITH)
  455.         ENDDO
  456. C
  457. C       Travail sur la partie superieure
  458.         LILIGN=IILIGS
  459.         DO ITH=1,NBTHR-1
  460.           CALL THREADID(ITH,SHOLE3I)
  461.         ENDDO
  462.         CALL SHOLE3I(NBTHR)
  463.         DO ITH=NBTHR-1,1,-1
  464.           CALL THREADIF(ITH)
  465.         ENDDO
  466. C
  467.         IF (I.EQ.IL1P.AND.IDER.NE.IL1P-1) THEN
  468.           DO IT=IDER,IPER,-1
  469.             LIG1=ILIGNS.ILIGR(IT)
  470.             SEGDES LIG1
  471.           ENDDO
  472.           IPER=IDER+1
  473.           IDER=IL1-1
  474.           DO IT=IPER,IDER
  475.             LIG1=ILIGNS.ILIGR(IT)
  476.             SEGACT /ERR=335/ LIG1
  477.           ENDDO
  478.   335     CONTINUE
  479.           IDER=IT-1
  480.           IF (IDER.LT.IPER) THEN
  481.             CALL ERREUR(48)
  482.             RETURN
  483.           ENDIF
  484.           GOTO 336
  485.         ENDIF
  486. C
  487.         NBTHR=1
  488. C
  489. C       Travail sur la partie superieure
  490.         MILIGN=IILIGS
  491.         LLIGN=MILIGN.ILIGN(I)
  492.         LIGN =ILIGNS.ILIGR(I)
  493.         NA=LIGN.IDERL-LIGN.IPREL+1
  494.         CALL SOMPAC(IPPVV(1),IMASQ(1),NA,KIVPO(1),KIVLO(1),NBPAR,IZROSF)
  495. C
  496.         NVALL=KIVLO(NBPAR)-1
  497.         LMASQ=MASQA(LLIGN.NJMAX)+1
  498. C
  499.   444   CONTINUE
  500. C
  501. C       LIGN est deja dimensionne a NJMAX
  502.         IF (ngmpet) THEN
  503.           SEGADJ,LIGN
  504.           LIG5=LIGN
  505.           GOTO 17
  506.         ENDIF
  507. C
  508.         SEGINI /ERR=109/ LIG5
  509.         GOTO 108
  510.  109    CONTINUE
  511.         IF (LSGDES) GOTO 14
  512. C       write(6,*) 'desactivation-3 ',1,il1p-1
  513.         lsgdes=.true.
  514.         do it=il1p-1,1,-1
  515.           lig1=ILIGNS.ILIGR(it)
  516.           segdes lig1
  517.           lig1=MILIGN.ilign(it)
  518.           segdes lig1
  519.         enddo
  520.         SEGINI /ERR=14/ LIG5
  521.  108    CONTINUE
  522.         LIG5.IML  =LIGN.IML
  523.         LIG5.IMM  =NBPAR
  524.         LIG5.IPREL=LIGN.IPREL
  525.         LIG5.IDERL=LIGN.IDERL
  526.         DO 111 LHE=1,NA
  527.           LIG5.IMMM(LHE)=LIGN.IMMM(LHE)
  528.  111    CONTINUE
  529.         DO 112 LHF=1,NA+1
  530.           LIG5.IPPVV(LHF)=LIGN.IPPVV(LHF)
  531.  112    CONTINUE
  532. C
  533.  17     CONTINUE
  534. C
  535.         DO 113 LHG=1,NBPAR
  536.           LIG5.IVPO(2*LHG-1)=KIVPO(LHG)
  537.           LIG5.IVPO(2*LHG)  =KIVLO(LHG)
  538.  113    CONTINUE
  539. C
  540.         DO IT=NBPAR-1,1,-1
  541.           IVI=LIG5.IVPO(2*IT)
  542.           IVF=LIG5.IVPO(2*(IT+1))-1
  543.           ICI=LIG5.IVPO(2*IT-1)
  544.           DO IC=MASQA(IVF-IVI+ICI),MASQA(ICI),-1
  545.             LIG5.IMASQ(IC)=IT
  546.           ENDDO
  547.         ENDDO
  548. C
  549.         IPA=1
  550.         DO 119 JPA=1,NA
  551.           IPP=LLIGN.IPPO(JPA)
  552.           KPA=LLIGN.IPPO(JPA+1)-IPP
  553.           LPA=LLIGN.LDEB(JPA)
  554.           LPA1=LPA-IPA
  555.           DO 120 MPA=1,KPA
  556.             XXV=LLIGN.XXVA(MPA+IPP)
  557.             IF (ABS(XXV).LE.XPETIT) GOTO 120
  558.             IPLA=LLIGN.LINC(MPA+IPP)-LPA1
  559.             IMSQ=LIG5.IMASQ(MASQA(IPLA))
  560.             IF (IMSQ.GT.0) THEN
  561. 114           CONTINUE
  562.               IF (IPLA.GE.LIG5.IVPO(2*(IMSQ+1)-1)) THEN
  563.                 IMSQ=IMSQ+1
  564.                 GOTO 114
  565.               ENDIF
  566.               IDBC=LIG5.IVPO(2*IMSQ-1)
  567.               IDBV=LIG5.IVPO(2*IMSQ  )
  568.               IPOS=IDBV+(IPLA-IDBC)
  569.               LIG5.VAL(IPOS)=XXV
  570.             ENDIF
  571. 120       CONTINUE
  572.           IPA=IPA+LLIGN.IMMMM(JPA)-LPA+1
  573. 119     CONTINUE
  574. C
  575.         ITR=NBPAR
  576.         DO 125 IPL=LMASQ,1,-1
  577.           IF (LIG5.IMASQ(IPL).GT.0) THEN
  578.             ITR=LIG5.IMASQ(IPL)
  579.           ELSE
  580.             LIG5.IMASQ(IPL)=-ITR
  581.           ENDIF
  582.  125    CONTINUE
  583. C
  584.         SEGSUP,LLIGN
  585.         MILIGN.ILIGN(I)=LIG5
  586.         ILIGNS.ILIGR(I)=0
  587. C
  588. C       Travail sur la partie inferieure
  589.         IF (MILIGN.EQ.IILIGS) THEN
  590.           MILIGN=IILIGN
  591.           LLIGN=MILIGN.ILIGN(I)
  592.           GOTO 444
  593.         ENDIF
  594.         IF (.not.ngmpet) SEGSUP,LIGN
  595. C
  596.         IDER=I
  597.         IPER=IDER
  598.         IL1=I+1
  599.   13  CONTINUE
  600.       GOTO 15
  601.   14  CONTINUE
  602.       LIGN=ILIGNS.ILIGR(I)
  603.       SEGSUP,LIGN
  604.       IL2=I-1
  605.   15  CONTINUE
  606.       IL1=IL1P
  607. C
  608.       if (lsgdes) then
  609. C       write(6,*) 'desactivation-4 ',1,il1-1
  610.         do it=il1-1,imin,-1
  611.           lig1=ILIGNS.ILIGR(it)
  612.           if (lig1.ne.0) segdes lig1
  613.         enddo
  614.       endif
  615.       nbthr=nbthrp
  616. C
  617. C  ****  BOUCLE  *5*  TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  618. C
  619.       ipos=0
  620.       iper=imin
  621.       ider=imin-1
  622.       DO 5 I=IMIN,IL2
  623.  
  624.       LIGM=LILIGN.ILIGN(I)
  625.       LIGL=MILIGN.ILIGN(I)
  626.       IF(I.LT.IL1) GOTO 7
  627. C
  628. C  ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  629. C  ******* IPREL  IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  630. C  ******* LE TERME DIAGONAL
  631. C
  632.       LIGN=LIGM
  633.       NBCOL=LIGN.IVPO(2*LIGN.IPPVV(2)-1)-1
  634.       DO 156 KHG=1,LIGN.IMMM(/1)
  635. C
  636. C       Travail sur la partie superieure
  637.         LIGN=LIGM
  638.         II=LIGN.IPREL-1+KHG
  639.         LIGN.IMMM(KHG)=0
  640. C
  641.         MM=LIGN.IVPO(2*LIGN.IPPVV(KHG))
  642.         NN=LIGN.IVPO(2*LIGN.IPPVV(KHG+1))
  643.         N=NN-MM
  644.         NN=NN-1
  645. C
  646.         DIAG(II)=LIGN.VAL(NN)
  647.         IF (N.EQ.1) THEN
  648.           IF (II-NBNNMC.GT.0) THEN
  649.             RE(II-NBNNMC,II-NBNNMC,1)=LIGN.VAL(NN)
  650.             GOTO 41
  651.           ENDIF
  652.           GOTO 8
  653.         ENDIF
  654. C
  655.         LIGN.VAL(NN)=LIGN.VAL(NN)+
  656.      &               CHOLE5(LIGN,LIGL,KHG,DIAG(1),NBCOL,0,NBOP(1))
  657.    8    CONTINUE
  658. C
  659. C       Travail sur la partie inferieure + terme diagonal
  660.         LIGN=LIGL
  661.         II=LIGN.IPREL-1+KHG
  662.         LIGN.IMMM(KHG)=0
  663. C
  664.         MM=LIGN.IVPO(2*LIGN.IPPVV(KHG))
  665.         NN=LIGN.IVPO(2*LIGN.IPPVV(KHG+1))
  666.         N=NN-MM
  667.         NN=NN-1
  668. C
  669.         LIGN.VAL(NN)=LIGN.VAL(NN)+
  670.      &               CHOLE5(LIGN,LIGM,KHG,DIAG(1),NBCOL,1,NBOP(1))
  671.         IF (II-NBNNMC.GT.0) THEN
  672.           RE(II-NBNNMC,II-NBNNMC,1)=LIGN.VAL(NN)
  673.           GOTO 41
  674.         ENDIF
  675. C
  676. C       Terme diagonal ok?
  677.         diagref=max(abs(diag(ii)),diagmin)
  678.         diagcmp=diagref*5d-12
  679.         IF (ABS(LIGN.VAL(NN)).GT.diagcmp) GOTO 12
  680. C
  681. C  il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement
  682. C  sur le calcul des reactions en cas de 2 relations presque identique
  683. C
  684. C  **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  685. C  **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  686. C  **** AU PREALABLE SUR CETTE INCONNUE.
  687. C
  688. *     write (6,*) ' ldmts mode d ensemble ittr ligne ',
  689. *    >      ittr(ii),ii,diag(ii),val(nn),diagref
  690. C  on garde le signe car il fau un moins sur les ML
  691.         vmaxi=diagref
  692.         do ipv=MM,NN
  693.           vmaxi=max(vmaxi,abs(lign.val(ipv)))
  694.         enddo
  695.         if(ittr(ii).NE.0) then
  696.           lign.VAL(NN)=lign.val(nn)-1.30D0*diagref
  697.           NENSLX=NENSLX+1
  698.         else
  699.           lign.val(nn)=vmaxi
  700.         endif
  701.         NENS=NENS+1
  702.         LIGL.IMMM(KHG)=NENS
  703.         LIGM.IMMM(KHG)=NENS
  704.    12   CONTINUE
  705. C
  706. C       Stabilisation
  707.         IF (ISTAB.NE.0) THEN
  708. C
  709. C         Elimination des Nan
  710.           if(.not.(abs(LIGN.val(nn)).lt.xgrand*xzprec).and.pasfait) then
  711.             pasfait=.false.
  712.             LIGN.val(nn)=xgrand*xzprec
  713.             write (6,*) 'Nan dans ldmts ligne',ii
  714.           endif
  715. C
  716.           diagcmp=abs(diagmin)*1d-5+xpetit/xzprec
  717.           IF (ITTR(II).EQ.0) THEN
  718.             IF (LIGN.VAL(NN).LT.-DIAGMAX*1D-3) THEN
  719.               LIGN.VAL(NN)=ABS(LIGN.VAL(NN))
  720. ***           val(nn)=diagmax*1d-3
  721.               IF (LIGN.IMMM(KHG).EQ.0) NENS=NENS+1
  722.               LIGN.IMMM(KHG)=NENS
  723.             ELSEIF (LIGN.VAL(NN).LE.DIAGCMP) THEN
  724.               IF (ISR.EQ.0.OR.IIMPI.GT.0)
  725.      &        write (ioimp,*) 'stabilisation RESO',ii,val(nn),diag(ii)
  726.               isr=isr+1
  727.               LIGN.VAL(NN)=MAX(DIAGCMP,-LIGN.VAL(NN))
  728.               IF (LIGN.IMMM(KHG).EQ.0) NENS=NENS+1
  729.               LIGN.IMMM(KHG)=NENS
  730.             ELSE
  731.               DIAGMIN=MIN(DIAGMIN,ABS(LIGN.VAL(NN)))
  732.             ENDIF
  733.           ELSE IF (ITTR(II).NE.0) THEN
  734.             IF (LIGN.VAL(NN).GE.ABS(DIAG(II))*STMULT) THEN
  735.               IF (ISRL.EQ.0.OR.IIMPI.GT.0)
  736.      &        write (ioimp,*) 'stabilisation RESO lagrange',ii,val(nn)
  737.               isrl=isrl+1
  738.               LIGN.VAL(NN)=-ABS(DIAG(II))*STMULT
  739.               IF (LIGN.IMMM(KHG).EQ.0) THEN
  740.                 NENS=NENS+1
  741.                 NENSLX=NENSLX+1
  742.               ENDIF
  743.               LIGN.IMMM(KHG)=NENS
  744.             ENDIF
  745.           ENDIF
  746.         ENDIF
  747. C
  748.         DIAG(II)=LIGN.VAL(NN)
  749.         IF(ABS(DIAG(II)).LE.XPETIT) THEN
  750.           DIAG(II)=diagmax
  751.           IF(ITTR(II).NE.0) DIAG(II)=-DIAGMAX
  752.           LIGL.VAL(NN)=DIAG(II)
  753.         ENDIF
  754.         LIGM.VAL(NN)=DIAG(II)
  755.         GOTO 41
  756. C
  757. C  ****  ERREUR MATRICE SINGUIERE
  758. C
  759. C CB215821 : No PATH to do this STATEMENT
  760.         INTERR(1)=I
  761.         CALL ERREUR(49)
  762.         if (ithrd.eq.1) then
  763.           call threadis
  764.          call oooprl(0)
  765.         endif
  766.         call ooomru(0)
  767.         RETURN
  768. C
  769. C  ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  770. C        ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  771. C
  772.    41   CONTINUE
  773.         IF(DIAG(II).LT.0.D0) INEG=INEG+1
  774.         IF(ITTR(II).NE.0)   NBLAG=NBLAG+1
  775. C
  776.         CONDMAX=MAX(CONDMAX,ABS(DIAG(II)))
  777.         IF (II.LE.NBNNMC) THEN
  778.           CONDMIN=MIN(CONDMIN,ABS(DIAG(II)))
  779.           DIAG(II)=1.D0/DIAG(II)
  780.         ENDIF
  781. C
  782.   156 CONTINUE
  783. C
  784. C     Suppression du masque
  785.       LMASQ=0
  786. C
  787.       LIGN=LIGL
  788.       NVALL=LIGN.VAL(/1)
  789.       NBPAR=LIGN.IVPO(/1)/2
  790.       NA   =LIGN.IMMM(/1)
  791.       SEGADJ,LIGN
  792. C
  793.       LIGN=LIGM
  794.       NVALL=LIGN.VAL(/1)
  795.       NBPAR=LIGN.IVPO(/1)/2
  796.       NA   =LIGN.IMMM(/1)
  797.       SEGADJ,LIGN
  798. C
  799.       NVSTOR=NVSTOR+NVALL
  800.       NVSTRM=MAX(NVSTRM,NVALL*4)
  801. C
  802.       NVALL=0
  803.       SEGINI,LIG3
  804.       lig3.iml=iml
  805.       lig3.imm=imm
  806.       lig3.iprel=iprel
  807.       lig3.iderl=iderl
  808.       do ipv=1,immm(/1)
  809.        lig3.immm(ipv)=immm(ipv)
  810.       enddo
  811.       do ipv=1,ippvv(/1)
  812.        lig3.ippvv(ipv)=ippvv(ipv)
  813.       enddo
  814.       do ipv=1,ivpo(/1)
  815.        lig3.ivpo(ipv)=ivpo(ipv)
  816.       enddo
  817.       ILIGNS.ILIGR(I)=LIG3
  818.       IF (lsgdes) SEGDES,LIG3
  819. C
  820. C      IF (I.EQ.IDBGL) THEN
  821. C        LIG1=LILIGN.ILIGN(I)
  822. C        WRITE(*,*) 'LIGNE FINALE SUP',I,LIG1
  823. C        segprt,lig1
  824. C        LIG1=MILIGN.ILIGN(I)
  825. C        WRITE(*,*) 'LIGNE FINALE INF',I,LIG1
  826. C        segprt,lig1
  827. C      ENDIF
  828. C
  829. C  ****  ON TRIANGULARISE LES AUTRES LIGNES
  830. C
  831.       IL1=IL1+1
  832.       IF (IL1.GT.IL2) GOTO 5
  833.       LIGM=LILIGN.ILIGN(I)
  834.       LIGL=MILIGN.ILIGN(I)
  835.       lign=MILIGN.ilign(il1)
  836.       GOTO 7
  837. C
  838.   71  CONTINUE
  839. C
  840.       IF (IPER.GT.IDER) THEN
  841.         call erreur(48)
  842.         write(6,*) ' 1 ldmts iper ider i il1 il2 ',
  843.      >         iper,ider,i,il1,il2
  844.         do ipv=1,ino
  845.           lig1=MILIGN.ilign(ipv)
  846.           call oooeta(lig1,ieta,imod)
  847.           if (ipv.lt.il1.or.ipv.gt.il2) then
  848. **         if (ieta.eq.1) write(6,*) 'ligne active ',ipv
  849.           endif
  850.         enddo
  851. **      call ooodmp(0)
  852.         if (ithrd.eq.1) then
  853.           call threadis
  854.           call oooprl(0)
  855.         endif
  856.         call ooomru(0)
  857.         RETURN
  858.       ENDIF
  859. C
  860. C     Passage en superlent
  861.       if (ider.lt.il1-1.and..not.ngmpet) then
  862.         write(6,*) 'passage en superlent',ider,i,il1-1
  863.         ngmpet=.true.
  864.       endif
  865. C
  866. C     soit parce qu'on a fini, soit parce qu'on manque de memoire
  867. C     il faut executer les lignes activees puis les desactiver
  868. C     lancer les chole4 et attendre qu'ils soient finis
  869.       if (xreser.ne.0) then
  870. CC      write(6,*)'segsup xreser 2'
  871.         segsup xreser
  872.         xreser=0
  873.       endif
  874.  
  875.       if (ipos.ne.0) then
  876. **      write (6,*) ' lancement thread ',iper,ider,il1,il2
  877.         nbthr=min(nbthr,il2-il1+1)
  878.         ipers=iper
  879.         iders=ider
  880.         ipert=iper
  881.         ifois=1
  882. *       write(6,*) 'il2 kidepb ider',lcara(1,il2),kidepb,ider
  883. **      write(6,*) 'i il1 il2',i,il1,il2
  884.  401    continue
  885.         ipas=max(1200/ifois,80)
  886.         if (ifois.ge.6 ) ipas= 80
  887.         if (nbthr.eq.1) ipas=igrand
  888. *402    continue
  889.         liper=lcara(2,iper)
  890.         iper2=il1-1
  891. **      call ooosur(MILIGN.ilign(il2))
  892.         do il=il1,il2
  893.           lign=MILIGN.ilign(il)
  894.           kidepb=lcara(1,il)-1
  895.           iperC=liper-kidepb
  896.           iperC=max(1,iperC)
  897.           imsq=LIGN.imasq(masqa(iperC))
  898.           if (imsq.lt.0) then
  899. CC          iperC=max(LIGN.ivpo(2*(-imsq)-1),iperC+1)
  900.             iperC=LIGN.ivpo(2*(-imsq)-1)
  901.           endif
  902.           iper1=ipno(iperC+kidepb)
  903.           iper2=min(iper1,iper2)
  904.         enddo
  905.         iper=iper2
  906.         ipert=iper
  907.         ider=min(iders,iper+ipas-1)
  908.         lider=lcara(3,ider)
  909.         ider2=il1-1
  910.         isaug=igrand
  911.         if (il2.lt.il1) then
  912.           write(6,*) 'il1 > il2',il1,il2
  913.           call erreur(48)
  914.         endif
  915.         do il=il1,il2
  916.           lign=MILIGN.ilign(il)
  917.           kidepb=lcara(1,il)-1
  918.           iderC=lider-kidepb
  919.           iderC=max(1,iderC)
  920.           imsq=LIGN.imasq(masqa(iderC))
  921.           isaut=0
  922.           if (imsq.lt.0) then
  923.             iderC=LIGN.ivpo(2*(-imsq)-1)-1
  924. **          isaut=0
  925. **          do mpv=masqa(iderC),masqa(iper-kidepb),-1
  926. **           if (LIGN.imasq(mpv).lt.0) isaut=isaut+masdim
  927. **          enddo
  928.             isaut=ipas/2
  929. **          isaut=-1
  930.           endif
  931.           ider1=ipno(iderC+kidepb)
  932.           ider2=min(ider1,ider2)
  933.           isaug=min(isaug,isaut)
  934.         enddo
  935.         ider=min(ider2+isaug,iders)
  936.         ifois=ifois+1
  937. C
  938. C       Travail sur la partie superieure
  939.         MILIGN=IILIGS
  940.         LILIGN=IILIGN
  941.         iths=nbthr
  942.         do ith=1,nbthr
  943.           if (ith.ne.iths) call threadid(ith,chole4i)
  944.         enddo
  945.         call chole4i(iths)
  946.         do ith=nbthr,1,-1
  947.           if (ith.ne.iths) call threadif(ith)
  948.         enddo
  949. C
  950. C       Travail sur la partie inferieure
  951.         MILIGN=IILIGN
  952.         LILIGN=IILIGS
  953.         do ith=1,nbthr
  954.           if (ith.ne.iths) call threadid(ith,chole4i)
  955.         enddo
  956.         call chole4i(iths)
  957.         do ith=nbthr,1,-1
  958.           if (ith.ne.iths) call threadif(ith)
  959.         enddo
  960. C
  961.         iper=ider+1
  962.         if (iper.le.iders) goto 401
  963.         iper=ipers
  964.         ider=iders
  965.       ENDIF
  966.  
  967. C  test ctrlC
  968.       if (ierr.ne.0) goto 9999
  969.       ipos=0
  970.       IF (LSGDES) THEN
  971. C       write(6,*) 'desactivation 5 iper ider',iper,ider
  972.         do ipv=ider,iper,-1
  973.           lig1=MILIGN.ilign(ipv)
  974.           segdes lig1
  975.           lig1=LILIGN.ilign(ipv)
  976.           segdes lig1
  977.         enddo
  978.         IF (I.NE.IL1-1) THEN
  979.           LIG1=MILIGN.ILIGN(I)
  980.           SEGACT,LIG1*MOD
  981.           LIG1=LILIGN.ILIGN(I)
  982.           SEGACT,LIG1*MOD
  983.         ENDIF
  984.       ENDIF
  985.       iper=ider+1
  986.       GOTO 5
  987. C
  988.     7 CONTINUE
  989.       IF (LSGDES.AND.(I.LT.IL1)) THEN
  990.         IF (.not.ngmpet) THEN
  991.           iperI=lcara(2,i)
  992.           do il=il1,il2
  993.             lign=MILIGN.ilign(il)
  994.             kidepb=lcara(1,il)-1
  995.             iperC=iperI-kidepb
  996.             iperC=max(1,iperC)
  997.             if (LIGN.imasq(masqa(iperC)).gt.0) goto 432
  998.           enddo
  999.           GOTO 5
  1000.         ENDIF
  1001.  432    CONTINUE
  1002.         SEGACT/ERR=71/LIGM
  1003.         SEGACT/ERR=71/LIGL
  1004.       ENDIF
  1005.       ipos=ipos+1
  1006.       ider=i
  1007.       IF (I.EQ.IL1-1) GOTO 71
  1008. C
  1009.     5 CONTINUE
  1010. C
  1011. **    write (6,*) ' il1 il2 apres 5 ',il1,il2
  1012.       if (lsgdes) then
  1013. C       write(6,*) 'desactivation 8 il1 il2',il1,il1p,il2
  1014.         DO I=IL2,IL1P,-1
  1015.           LIGN=MILIGN.ILIGN(I)
  1016.           SEGDES,LIGN
  1017.           LIGN=LILIGN.ILIGN(I)
  1018.           SEGDES,LIGN
  1019.         ENDDO
  1020.       ENDIF
  1021.       nbopt=0
  1022.       do ith=1,nbthro
  1023.        nbopt=nbopt+nbop(ith)
  1024.        nbop(ith)=0
  1025.       enddo
  1026.       nbopin=nbopt
  1027.       nbopit=nbopit+nbopin
  1028.       call timespv(ittime,oothrd)
  1029.       kcour=(ittime(1)+ittime(2))/10
  1030.       if (xreser.ne.0) then
  1031. CC      write(6,*)'segsup xreser 3'
  1032.         segsup xreser
  1033.         xreser=0
  1034.       endif
  1035.       IF(IL2.LT.INO) GOTO 1
  1036. C  ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  1037. C   ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  1038. C       INEG=INEG-NBLAG
  1039. C     on ne compte pas 2 fois les multiplicateurs qui vont etre
  1040. C     elimines lors de la resolution car mode d'ensemble
  1041.       INEG=INEG-(NBLAG-NENSLX)
  1042.       if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  1043.  4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  1044.      1' DE LAGRANGE DETECTES')
  1045.  
  1046.       IF (IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC
  1047.  4821 FORMAT( ' NOMBRE DE VALEURS  DANS LE PROFIL',I12)
  1048.       IF (IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR*2
  1049.  4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  1050.       IF (IIMPI.EQ.1) WRITE(IOIMP,4825) Nbopit/1000000
  1051.  4825 FORMAT( ' NOMBRE DE GIGA OPERATIONS FMA',I40)
  1052.       IF (IIMPI.EQ.1) WRITE(IOIMP,4823) NVaor
  1053.  4823 FORMAT( ' NOMBRE DE VALEURS initiales',I9)
  1054. C     IF (IIMPI.EQ.1) WRITE(IOIMP,4824) nbopit/1d6/(kcour-kcouri)
  1055. C 4824 FORMAT( ' Performance en Gigaflops ',F8.1)
  1056.       INTERR(1)=NVSTOR*2
  1057.       REAERR(1)=NVSTOR/INC**(4.D0/3.D0)
  1058.       REAERR(2)=2.D0*NBOPIT/1.D6/MAX(1,(KCOUR-KCOURI))
  1059.       REAERR(3)=CONDMAX/CONDMIN
  1060.       IF (IPASV.EQ.0.or.reaerr(3).gt.1.D30) CALL ERREUR(-278)
  1061.       IPASV=1
  1062.       call ooomru(0)
  1063.       if(lsgdes) then
  1064.         do ipv=1,ino
  1065.           lign=MILIGN.ilign(ipv)
  1066.           segdes lign
  1067.           lign=LILIGN.ilign(ipv)
  1068.           segdes lign
  1069.         enddo
  1070.       endif
  1071.       SEGDES,MINCPO
  1072.       SEGDES,MIMIK
  1073.       SEGDES,MMATRI
  1074.       SEGDES,MILIGN,LILIGN
  1075.       MMATRX=MMATRI
  1076.       SEGSUP KIVPO,KIVLO
  1077.       do ipv=1,nnoe
  1078.         lig3=ILIGNS.ILIGR(ipv)
  1079.         segsup lig3
  1080.       enddo
  1081.       SEGSUP,ILIGNS
  1082.       segsup immt
  1083. C     write (6,*) ' ldmts ipos max ',iposm
  1084.  9999 continue
  1085.       if (ithrd.eq.1) then
  1086.         call threadis
  1087.        call oooprl(0)
  1088.       endif
  1089.       if(iimpi.ne.0)  write (6,*)
  1090.      >      ' ldmts condmin condmax ',condmin,condmax,diag(/1)
  1091.       SEGDES,MDIAG
  1092.       RETURN
  1093.       END
  1094.  
  1095.  
  1096.  
  1097.  
  1098.  
  1099.  
  1100.  
  1101.  

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