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    PV090527  26/04/30    21:15:49     12529          
  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*(NJMAX/NA)**1.333333333333333333)
  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.         INTERR(1)=I
  760.         CALL ERREUR(49)
  761.         if (ithrd.eq.1) then
  762.           call threadis
  763.          call oooprl(0)
  764.         endif
  765.         call ooomru(0)
  766.         RETURN
  767. C
  768. C  ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  769. C        ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  770. C
  771.    41   CONTINUE
  772.         IF(DIAG(II).LT.0.D0) INEG=INEG+1
  773.         IF(ITTR(II).NE.0)   NBLAG=NBLAG+1
  774. C
  775.         CONDMAX=MAX(CONDMAX,ABS(DIAG(II)))
  776.         IF (II.LE.NBNNMC) THEN
  777.           CONDMIN=MIN(CONDMIN,ABS(DIAG(II)))
  778.           DIAG(II)=1.D0/DIAG(II)
  779.         ENDIF
  780. C
  781.   156 CONTINUE
  782. C
  783. C     Suppression du masque
  784.       LMASQ=0
  785. C
  786.       LIGN=LIGL
  787.       NVALL=LIGN.VAL(/1)
  788.       NBPAR=LIGN.IVPO(/1)/2
  789.       NA   =LIGN.IMMM(/1)
  790.       SEGADJ,LIGN
  791. C
  792.       LIGN=LIGM
  793.       NVALL=LIGN.VAL(/1)
  794.       NBPAR=LIGN.IVPO(/1)/2
  795.       NA   =LIGN.IMMM(/1)
  796.       SEGADJ,LIGN
  797. C
  798.       NVSTOR=NVSTOR+NVALL
  799.       NVSTRM=MAX(NVSTRM,NVALL*4)
  800. C
  801.       NVALL=0
  802.       SEGINI,LIG3
  803.       lig3.iml=iml
  804.       lig3.imm=imm
  805.       lig3.iprel=iprel
  806.       lig3.iderl=iderl
  807.       do ipv=1,immm(/1)
  808.        lig3.immm(ipv)=immm(ipv)
  809.       enddo
  810.       do ipv=1,ippvv(/1)
  811.        lig3.ippvv(ipv)=ippvv(ipv)
  812.       enddo
  813.       do ipv=1,ivpo(/1)
  814.        lig3.ivpo(ipv)=ivpo(ipv)
  815.       enddo
  816.       ILIGNS.ILIGR(I)=LIG3
  817.       IF (lsgdes) SEGDES,LIG3
  818. C
  819.       IF (I.EQ.IDBGL) THEN
  820.         LIG1=LILIGN.ILIGN(I)
  821.         WRITE(*,*) 'LIGNE FINALE SUP',I,LIG1
  822.         segprt,lig1
  823.         LIG1=MILIGN.ILIGN(I)
  824.         WRITE(*,*) 'LIGNE FINALE INF',I,LIG1
  825.         segprt,lig1
  826.       ENDIF
  827. C
  828. C  ****  ON TRIANGULARISE LES AUTRES LIGNES
  829. C
  830.       IL1=IL1+1
  831.       IF (IL1.GT.IL2) GOTO 5
  832.       LIGM=LILIGN.ILIGN(I)
  833.       LIGL=MILIGN.ILIGN(I)
  834.       lign=MILIGN.ilign(il1)
  835.       GOTO 7
  836. C
  837.   71  CONTINUE
  838. C
  839.       IF (IPER.GT.IDER) THEN
  840.         call erreur(48)
  841.         write(6,*) ' 1 ldmts iper ider i il1 il2 ',
  842.      >         iper,ider,i,il1,il2
  843.         do ipv=1,ino
  844.           lig1=MILIGN.ilign(ipv)
  845.           call oooeta(lig1,ieta,imod)
  846.           if (ipv.lt.il1.or.ipv.gt.il2) then
  847. **         if (ieta.eq.1) write(6,*) 'ligne active ',ipv
  848.           endif
  849.         enddo
  850. **      call ooodmp(0)
  851.         if (ithrd.eq.1) then
  852.           call threadis
  853.           call oooprl(0)
  854.         endif
  855.         call ooomru(0)
  856.         RETURN
  857.       ENDIF
  858. C
  859. C     Passage en superlent
  860.       if (ider.lt.il1-1.and..not.ngmpet) then
  861.         write(6,*) 'passage en superlent',ider,i,il1-1
  862.         ngmpet=.true.
  863.       endif
  864. C
  865. C     soit parce qu'on a fini, soit parce qu'on manque de memoire
  866. C     il faut executer les lignes activees puis les desactiver
  867. C     lancer les chole4 et attendre qu'ils soient finis
  868.       if (xreser.ne.0) then
  869. CC      write(6,*)'segsup xreser 2'              
  870.         segsup xreser
  871.         xreser=0
  872.       endif
  873.  
  874.       if (ipos.ne.0) then
  875. **      write (6,*) ' lancement thread ',iper,ider,il1,il2
  876.         nbthr=min(nbthr,il2-il1+1)
  877.         ipers=iper
  878.         iders=ider
  879.         ipert=iper
  880.         ifois=1
  881. *       write(6,*) 'il2 kidepb ider',lcara(1,il2),kidepb,ider
  882. **      write(6,*) 'i il1 il2',i,il1,il2
  883.  401    continue
  884.         ipas=max(1200/ifois,80)
  885.         if (ifois.ge.6 ) ipas= 80
  886.         if (nbthr.eq.1) ipas=igrand
  887. *402    continue
  888.         liper=lcara(2,iper)
  889.         iper2=il1-1
  890. **      call ooosur(MILIGN.ilign(il2))
  891.         do il=il1,il2
  892.           lign=MILIGN.ilign(il)
  893.           kidepb=lcara(1,il)-1
  894.           iperC=liper-kidepb
  895.           iperC=max(1,iperC)
  896.           imsq=LIGN.imasq(masqa(iperC))
  897.           if (imsq.lt.0) then
  898. CC          iperC=max(LIGN.ivpo(2*(-imsq)-1),iperC+1)
  899.             iperC=LIGN.ivpo(2*(-imsq)-1)
  900.           endif
  901.           iper1=ipno(iperC+kidepb)
  902.           iper2=min(iper1,iper2)
  903.         enddo
  904.         iper=iper2
  905.         ipert=iper
  906.         ider=min(iders,iper+ipas-1)
  907.         lider=lcara(3,ider)
  908.         ider2=il1-1
  909.         isaug=igrand
  910.         if (il2.lt.il1) then
  911.           write(6,*) 'il1 > il2',il1,il2
  912.           call erreur(48)
  913.         endif
  914.         do il=il1,il2
  915.           lign=MILIGN.ilign(il)
  916.           kidepb=lcara(1,il)-1
  917.           iderC=lider-kidepb
  918.           iderC=max(1,iderC)
  919.           imsq=LIGN.imasq(masqa(iderC))
  920.           isaut=0
  921.           if (imsq.lt.0) then
  922.             iderC=LIGN.ivpo(2*(-imsq)-1)-1
  923. **          isaut=0
  924. **          do mpv=masqa(iderC),masqa(iper-kidepb),-1
  925. **           if (LIGN.imasq(mpv).lt.0) isaut=isaut+masdim
  926. **          enddo
  927.             isaut=ipas/2
  928. **          isaut=-1
  929.           endif
  930.           ider1=ipno(iderC+kidepb)
  931.           ider2=min(ider1,ider2)
  932.           isaug=min(isaug,isaut)
  933.         enddo
  934.         ider=min(ider2+isaug,iders)
  935.         ifois=ifois+1
  936. C
  937. C       Travail sur la partie superieure
  938.         MILIGN=IILIGS
  939.         LILIGN=IILIGN
  940.         iths=nbthr
  941.         do ith=1,nbthr
  942.           if (ith.ne.iths) call threadid(ith,chole4i)
  943.         enddo
  944.         call chole4i(iths)
  945.         do ith=nbthr,1,-1
  946.           if (ith.ne.iths) call threadif(ith)
  947.         enddo
  948. C
  949. C       Travail sur la partie inferieure
  950.         MILIGN=IILIGN
  951.         LILIGN=IILIGS
  952.         do ith=1,nbthr
  953.           if (ith.ne.iths) call threadid(ith,chole4i)
  954.         enddo
  955.         call chole4i(iths)
  956.         do ith=nbthr,1,-1
  957.           if (ith.ne.iths) call threadif(ith)
  958.         enddo
  959. C
  960.         iper=ider+1
  961.         if (iper.le.iders) goto 401
  962.         iper=ipers
  963.         ider=iders
  964.       ENDIF
  965.  
  966. C  test ctrlC
  967.       if (ierr.ne.0) goto 9999
  968.       ipos=0
  969.       IF (LSGDES) THEN
  970. C       write(6,*) 'desactivation 5 iper ider',iper,ider
  971.         do ipv=ider,iper,-1
  972.           lig1=MILIGN.ilign(ipv)
  973.           segdes lig1
  974.           lig1=LILIGN.ilign(ipv)
  975.           segdes lig1
  976.         enddo
  977.         IF (I.NE.IL1-1) THEN
  978.           LIG1=MILIGN.ILIGN(I)
  979.           SEGACT,LIG1*MOD
  980.           LIG1=LILIGN.ILIGN(I)
  981.           SEGACT,LIG1*MOD
  982.         ENDIF
  983.       ENDIF
  984.       iper=ider+1
  985.       GOTO 5
  986. C
  987.     7 CONTINUE
  988.       IF (LSGDES.AND.(I.LT.IL1)) THEN
  989.         IF (.not.ngmpet) THEN
  990.           iperI=lcara(2,i)
  991.           do il=il1,il2
  992.             lign=MILIGN.ilign(il)
  993.             kidepb=lcara(1,il)-1
  994.             iperC=iperI-kidepb
  995.             iperC=max(1,iperC)
  996.             if (LIGN.imasq(masqa(iperC)).gt.0) goto 432
  997.           enddo
  998.           GOTO 5
  999.         ENDIF
  1000.  432    CONTINUE
  1001.         SEGACT/ERR=71/LIGM
  1002.         SEGACT/ERR=71/LIGL
  1003.       ENDIF
  1004.       ipos=ipos+1
  1005.       ider=i
  1006.       IF (I.EQ.IL1-1) GOTO 71
  1007. C
  1008.     5 CONTINUE
  1009. C
  1010. **    write (6,*) ' il1 il2 apres 5 ',il1,il2
  1011.       if (lsgdes) then
  1012. C       write(6,*) 'desactivation 8 il1 il2',il1,il1p,il2
  1013.         DO I=IL2,IL1P,-1
  1014.           LIGN=MILIGN.ILIGN(I)
  1015.           SEGDES,LIGN
  1016.           LIGN=LILIGN.ILIGN(I)
  1017.           SEGDES,LIGN
  1018.         ENDDO
  1019.       ENDIF
  1020.       nbopt=0
  1021.       do ith=1,nbthro
  1022.        nbopt=nbopt+nbop(ith)
  1023.        nbop(ith)=0
  1024.       enddo
  1025.       nbopin=nbopt
  1026.       nbopit=nbopit+nbopin
  1027.       call timespv(ittime,oothrd)
  1028.       kcour=(ittime(1)+ittime(2))/10
  1029.       if (xreser.ne.0) then 
  1030. CC      write(6,*)'segsup xreser 3'              
  1031.         segsup xreser
  1032.         xreser=0
  1033.       endif
  1034.       IF(IL2.LT.INO) GOTO 1
  1035. C  ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  1036. C   ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  1037. C       INEG=INEG-NBLAG
  1038. C     on ne compte pas 2 fois les multiplicateurs qui vont etre
  1039. C     elimines lors de la resolution car mode d'ensemble
  1040.       INEG=INEG-(NBLAG-NENSLX)
  1041.       if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  1042.  4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  1043.      1' DE LAGRANGE DETECTES')
  1044.  
  1045.       IF (IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC
  1046.  4821 FORMAT( ' NOMBRE DE VALEURS  DANS LE PROFIL',I12)
  1047.       IF (IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR*2
  1048.  4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  1049.       IF (IIMPI.EQ.1) WRITE(IOIMP,4825) Nbopit/1000000
  1050.  4825 FORMAT( ' NOMBRE DE GIGA OPERATIONS FMA',I40)
  1051.       IF (IIMPI.EQ.1) WRITE(IOIMP,4823) NVaor
  1052.  4823 FORMAT( ' NOMBRE DE VALEURS initiales',I9)
  1053. C     IF (IIMPI.EQ.1) WRITE(IOIMP,4824) nbopit/1d6/(kcour-kcouri)
  1054. C 4824 FORMAT( ' Performance en Gigaflops ',F8.1)
  1055.       INTERR(1)=NVSTOR*2
  1056.       REAERR(1)=NVSTOR/INC**(4.D0/3.D0)
  1057.       REAERR(2)=2.D0*NBOPIT/1.D6/MAX(1,(KCOUR-KCOURI))
  1058.       REAERR(3)=CONDMAX/CONDMIN
  1059.       IF (IPASV.EQ.0.or.reaerr(3).gt.1.D30) CALL ERREUR(-278)
  1060.       IPASV=1
  1061.       call ooomru(0)
  1062.       if(lsgdes) then
  1063.         do ipv=1,ino
  1064.           lign=MILIGN.ilign(ipv)
  1065.           segdes lign
  1066.           lign=LILIGN.ilign(ipv)
  1067.           segdes lign
  1068.         enddo
  1069.       endif
  1070.       SEGDES,MINCPO
  1071.       SEGDES,MIMIK
  1072.       SEGDES,MMATRI
  1073.       SEGDES,MILIGN,LILIGN
  1074.       MMATRX=MMATRI
  1075.       SEGSUP KIVPO,KIVLO
  1076.       do ipv=1,nnoe
  1077.         lig3=ILIGNS.ILIGR(ipv)
  1078.         segsup lig3
  1079.       enddo
  1080.       SEGSUP,ILIGNS
  1081.       segsup immt
  1082. C     write (6,*) ' ldmts ipos max ',iposm
  1083.  9999 continue
  1084.       if (ithrd.eq.1) then
  1085.         call threadis
  1086.        call oooprl(0)
  1087.       endif
  1088.       if(iimpi.ne.0)  write (6,*) 
  1089.      >      ' ldmts condmin condmax ',condmin,condmax,diag(/1)
  1090.       SEGDES,MDIAG
  1091.       RETURN
  1092.       END
  1093.  
  1094.  
  1095.  
  1096.  
  1097.  
  1098.  

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