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    MB234859  26/03/30    21:15:03     12509          
  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,ILIGNS.MILIGN,ILIGNL.MILIGN,ILIGNM.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.       external chole4i
  33.       external shole3i
  34.       SAVE IPASV
  35.       DATA IPASV/0/
  36. *  ngmpet dit si on tient en memoire (false) ou si on deborde (true)
  37.       logical ngmpet
  38.       logical lsgdes,pasfait
  39. C
  40.       SEGDES,MCOORD
  41.       maitre=0                                                                              
  42.       xreser=0
  43.       pasfait=.true.
  44.       lsgdes=.false.
  45. *  faire attention a respecter l'ordre des segdes par la suite
  46.       call ooomru(1)
  47.       if (istab.ne.0.and.istab.ne.1) call erreur(5)
  48.       condmax=0.d0
  49.       condmin=xgrand*xzprec
  50.       ngmpet=.false.
  51.       nvallp=nvall
  52.       call timespv(ittime,oothrd)
  53.       kcour=(ittime(1)+ittime(2))/10
  54.       kcouri=kcour
  55.       kcour=0
  56.       nbopit=0
  57.       nvaor=0
  58.       nbthro=nbthrs
  59.       ithrd=0
  60.       if (nbthro.gt.1) then
  61.        ithrd=1
  62.        call threadii
  63.        call oooprl(1)
  64.       endif
  65.       nbthr=nbthro
  66.       do ith=1,nbthr
  67.        nbop(ith)=0
  68.       enddo
  69.       stmult=1d-5
  70.  
  71.       nvstrm=0             
  72.       ivstrm=oooval(1,4)*2
  73. C     write(6,*) 'nvstrm initial',nvstrm
  74.       MMATRI=MMATRX
  75.       SEGACT,MMATRI*MOD
  76.       PRCHLV=PREC
  77.       MILIGN=IILIGN
  78.       LILIGN=IILIGS
  79.       SEGACT,MILIGN*MOD,LILIGN*MOD
  80. C
  81.       INO=MILIGN.ILIGN(/1)
  82.       NNOE=INO
  83.       INC=0
  84.       SEGINI,ILIGNL,ILIGNM
  85. C
  86.       MDIAG=IDIAG
  87.       SEGACT,MDIAG*MOD
  88.       NBLIG=INO
  89.       segini immt
  90.       precc=prec
  91.       INC=DIAG(/1)
  92. C
  93. C     Cas du super-element
  94.       nbnnmc=inc+1
  95.       if (xmatri.ne.0) then
  96.         nelrig=1
  97.         nligrp=nligra
  98.         nligrd=nligra
  99.         segini xmatri
  100.         nbnnmc=nbnnma
  101.       endif
  102.       matric=xmatri
  103.       INCC=INC
  104.       MIMIK=IIMIK
  105.       MINCPO=IINCPO
  106.       SEGACT,MINCPO,MIMIK
  107.       IPLUMI=IMIK(/2)*2 +4
  108.       IL2=0
  109.       IIMAX=IJMAX+IPLUMI
  110.       SEGINI KIVPO,KIVLO
  111.       INEG=0
  112.       NBLAG=0
  113.       NENSLX=0
  114.       NVSTOC=0
  115.       NVSTOR=0
  116.       diagmax=XPETIT/XZPREC
  117.       diagmin=xgrand
  118.       do i=1,diag(/1)
  119.        if (ittr(i).eq.0) diagmax=max(diagmax,abs(diag(i)))
  120.        if (ittr(i).eq.0.and.abs(diag(i)).gt.xpetit/xzprec)
  121.      >   diagmin=min(diagmin,abs(diag(i)))
  122.       enddo
  123.       if (diagmax.le.xpetit/xzprec) then
  124.          do i=1,diag(/1)
  125.            diagmax=max(diagmax,abs(diag(i)))
  126.            if (abs(diag(i)).gt.xpetit/xzprec)
  127.      >       diagmin=min(diagmin,abs(diag(i)))
  128.          enddo
  129.       endif
  130.       diagmin=min(diagmin,diagmax)
  131. ***    write (6,*) ' ldmts diagmin diagmax ',diagmin,diagmax,diag(/1)
  132. C
  133. C  **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  134. C  **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  135. C
  136. C  ****  LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  137. C
  138. C SP  indicateurs pour impression message "stabilisation RESO..."
  139.       isr=0
  140.       isrl=0
  141.  
  142.   1   CONTINUE
  143.       IL1=IL2+1
  144.       IMIN=IL1
  145. C
  146. C     Reserver de la place ou mettre les lignes superieures dans le cas debordement
  147.       if (ngmpet) then
  148.         if (xreser.eq.0) then
  149. CC      write(6,*)'segini xreser',nvstrm
  150.         segini xreser
  151.         endif
  152.       endif
  153. C
  154.       DO 2 I=IL1,INO
  155. C
  156.       LIGL=0
  157.       LIGM=0
  158.       LGLIG=0
  159. C
  160.       LLIGM=LILIGN.ILIGN(I)
  161.       LLIGL=MILIGN.ILIGN(I)
  162.       SEGACT /ERR=32/ LLIGM
  163.       SEGACT /ERR=32/ LLIGL
  164.       GOTO 31
  165.   32  CONTINUE
  166.       IF (LSGDES) GOTO 3
  167. **    write(6,*) 'desactivation-1 ',1,il1-1
  168.       lsgdes=.true.
  169.       do it=il1-1,1,-1
  170.         lig1=ILIGNM.ILIGN(it)
  171.         segdes lig1
  172.         lig1=ILIGNL.ILIGN(it)
  173.         segdes lig1
  174.         lig1=LILIGN.ilign(it)
  175.         segdes lig1
  176.         lig1=MILIGN.ilign(it)
  177.         segdes lig1
  178.       enddo
  179.       SEGACT /ERR=3/LLIGM
  180.       SEGACT /ERR=3/LLIGL
  181.   31  CONTINUE
  182. ** essai pv
  183.       NVMAX=MAX(LLIGM.NJMAX,LLIGL.NJMAX)
  184.       if (ngmpet) then
  185.         if (nvmax.gt.2*njmaxp.and.i-il1+1.gt.nbthrs/2) then
  186. **        write(6,*) 'fin de bloc il1  i nvmax',il1,i,nvmax
  187.           goto 3
  188.         endif
  189.       endif
  190.       NJMAXP=NVMAX
  191. ** fin essai
  192. C
  193.       NA1  =LLIGM.IMMMM(/1)
  194.       NVAL1=LLIGM.IMMMM(NA1)-LLIGM.LDEB(1)
  195.       NA2  =LLIGL.IMMMM(/1)
  196.       NVAL2=LLIGL.IMMMM(NA2)-LLIGL.LDEB(1)
  197. C
  198.       NVALL=0
  199.       NVALT=MAX(NVAL1,NVAL2)
  200.       if (ngmpet.and.nvalt.lt.nvallp/3) then
  201.         write(6,*) 'passage en lent'
  202.         ngmpet=.false.
  203.       endif
  204.       NVALLP=NVALT
  205.       if (ngmpet) then
  206. C       nvall=(NVMAX*nvstor)/nvstoc
  207.         NVALL=NVMAX
  208.         if (nvall.gt.NVSTRM) then
  209. C       write(*,*)'redimensionnement de xreser',nvstrm,(4*NVALL),il1,il2
  210.           nvstrm=4*nvall
  211.           segsup,xreser
  212.           segini,xreser
  213.         endif
  214.       endif
  215. C
  216.       NA=NA1
  217.       NBPAR=NA+1
  218.       LMASQ=MASQA(NVAL1+MASDIM)
  219.       SEGINI /ERR=33/ LIGM
  220. C
  221.       NA=NA2
  222.       NBPAR=NA+1
  223.       LMASQ=MASQA(NVAL2+MASDIM)
  224.       SEGINI /ERR=33/ LIGL
  225. C
  226.       GOTO 34
  227.   33  CONTINUE
  228.       IF (LSGDES) GOTO 3
  229. **    write(6,*) 'desactivation-2 ',1,il1-1
  230.       lsgdes=.true.
  231.       do it=il1-1,1,-1
  232.         lig1=ILIGNM.ILIGN(it)
  233.         segdes lig1
  234.         lig1=ILIGNL.ILIGN(it)
  235.         segdes lig1
  236.         lig1=LILIGN.ilign(it)
  237.         segdes lig1
  238.         lig1=MILIGN.ilign(it)
  239.         segdes lig1
  240.       enddo
  241.       IF (LIGM.EQ.0) THEN
  242.         NA=NA1
  243.         NBPAR=NA+1
  244.         LMASQ=MASQA(NVAL1+MASDIM)
  245.         SEGINI /ERR=3/ LIGM
  246.       ENDIF
  247.       IF (LIGL.EQ.0) THEN
  248.         NA=NA2
  249.         NBPAR=NA+1
  250.         LMASQ=MASQA(NVAL2+MASDIM)
  251.         SEGINI /ERR=3/ LIGL
  252.       ENDIF
  253.   34  CONTINUE
  254. C
  255. C     Travail sur la partie superieure
  256.       MILIGN=IILIGS
  257.       ILIGNS=ILIGNM
  258.       LLIGN =LLIGM
  259.       LIGN  =LIGM
  260. C
  261.   35  CONTINUE
  262. C
  263. C     Recuperer la longueur du segment
  264.       NA    =LIGN.IMMM(/1)
  265.       LGLIG =LGLIG  + (NA*(NJMAX/NA)**1.333333333333333333)
  266.       NVSTOC=NVSTOC + LLIGN.NJMAX
  267.       NVAOR =NVAOR  + LLIGN.XXVA(/1)
  268. C
  269. C  ****  DECOMPACTAGE
  270. C
  271.       IPA=1
  272.       DO 121 JPA=1,NA
  273.         KPA =LLIGN.IPPO(JPA+1)-LLIGN.IPPO(JPA)
  274.         IPP =LLIGN.IPPO(JPA)
  275.         LPA =LLIGN.LDEB(JPA)
  276.         LPA1=LPA-IPA
  277.         LIGN.IMMM(JPA) =LPA
  278.         LIGN.IVPO(JPA) =IPA
  279.         LIGN.IPPVV(JPA)=IPA-1
  280.         DO 122 MPA=1,KPA
  281.           IPLA=LLIGN.LINC(MPA+IPP)-LPA1
  282.           ICC =IPLA-IPA+1
  283.           ICC1=ICC-MASQD(ICC)+1
  284.           MICC=MASQA(ICC)
  285.           IMSQ=LIGN.IMASQ(MICC)
  286.           IF (IMSQ.LE.0) THEN
  287.             MSQH=icc1
  288.             MSQB=icc1
  289.           ELSE
  290.             MSQH=MAX(MASQH(IMSQ),icc1)
  291.             MSQB=MIN(MASQB(IMSQ),icc1)
  292.           ENDIF
  293.           LIGN.IMASQ(MICC)=MASQV(MSQB,MSQH)
  294. 122     CONTINUE
  295.         IPA=IPA+LLIGN.IMMMM(JPA)-LPA + 1
  296.         IF(IMIN.GT.MILIGN.IPNO(LPA)) IMIN=MILIGN.IPNO(LPA)
  297. 121   CONTINUE
  298.       LIGN.IPPVV(NA+1)=IPA-1
  299.  
  300. C     Indexation de imasq
  301.       LMASQ=LIGN.IMASQ(/1)
  302.       ipln=lmasq
  303.       idec=1
  304.       do 123 ipl=lmasq,1,-1
  305.         if (LIGN.imasq(ipl).gt.0) then
  306.           ipln=ipl-1
  307.           idec=MASQB(LIGN.IMASQ(IPL))
  308.         else
  309.           LIGN.IMASQ(IPL)=-(IPLN*MASDIM+IDEC)
  310.         endif
  311.  123  continue
  312. CCC   write (6,*) ' imasq ',lmasq
  313. CCC   write (6,*) (LIGN.imasq(ipl),ipl=1,lmasq)
  314.  
  315.       IF (NA.GT.0) THEN
  316.         LIGN.IPREL=LLIGN.IMMMM(1)
  317.         LIGN.IDERL=LLIGN.IMMMM(NA)
  318.         MILIGN.LCARA(2,I)=LIGN.IPREL
  319.         MILIGN.LCARA(3,I)=LIGN.IDERL
  320.       ENDIF
  321. C
  322.       ILIGNS.ILIGN(I)=LIGN
  323. C
  324. C     Travail sur la partie inferieure
  325.       IF (ILIGNS.EQ.ILIGNM) THEN
  326.         MILIGN=IILIGN
  327.         ILIGNS=ILIGNL
  328.         LLIGN =LLIGL
  329.         LIGN  =LIGL
  330.         GOTO 35
  331.       ENDIF
  332. C
  333. C*    write (6,*) 'longueur ligne ',nvall
  334. C nb de ligne multiple du nb de threads
  335. C blocage ligne lecture-ecriture pour minimiser le cache
  336. C on note si on est au minimum de lignes
  337.       nbthro=min(nbthrs,lglig/500+1)
  338.       nbthro2=min(nbthrs,lglig/10000+1)
  339.       IF (.not.ngmpet) THEN
  340.         if (i+1-il1.ge.nbthro2) then
  341.           nbthro2=min(nbthrs,i+1-il1)
  342.         endif
  343.         if (i+1-il1.ge.nbthro) then
  344.           nbthro=min(nbthrs,i+1-il1)
  345.           il2=i
  346. **        write(6,*) ' nbthro il1 il2 ',nbthro,il1,il2
  347.           GOTO 4
  348.         endif
  349.       ENDIF
  350.  
  351.     2 CONTINUE
  352.       IL2=INO
  353.       GOTO 4
  354.     3 CONTINUE
  355.       IL2=I-1
  356. **    write(6,*) 'desactivation-4 ',llign
  357.       SEGDES LLIGM,LLIGL
  358.       IF (LIGL.NE.0) SEGSUP,LIGL
  359.       IF (LIGM.NE.0) SEGSUP,LIGM
  360.     4 CONTINUE
  361.       nbthro=min(nbthrs,nbthro)
  362.       nbthro2=min(nbthrs,nbthro2)
  363.       nbthr=nbthro
  364.       if(xreser.ne.0) then
  365. CC      write(6,*)'segsup xreser 1',il1,il2
  366.         segsup xreser
  367.         xreser=0
  368.       endif
  369. C
  370.       IF(IL2.GE.IL1) GOTO 40
  371. C
  372. C ****  MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  373. C
  374.       CALL ERREUR(48)
  375. **    call ooodmp(0)
  376.       if (ithrd.eq.1) then
  377.        call threadis
  378.        call oooprl(0)
  379.       endif
  380.       call ooomru(0)
  381.       RETURN
  382.    40 CONTINUE
  383. C
  384. C     Travail sur la partie superieure
  385.       MILIGN=IILIGS
  386.       ILIGNS=ILIGNM
  387.       LLIGN =LLIGM
  388.       LIGN  =LIGM
  389. C
  390.    39 CONTINUE
  391. C
  392.       IM=INC
  393.       DO 352 IH=IL2,IL1,-1
  394.         LIGN=ILIGNS.ILIGN(IH)
  395.         IL=INC
  396.         DO 354 JH=1,LIGN.IMMM(/1)
  397.           IM=MIN(IM,LIGN.IMMM(JH))
  398.           IL=MIN(IL,LIGN.IMMM(JH))
  399.  354    CONTINUE
  400.         LIGN.IML=IL
  401.         MILIGN.LCARA(1,IH)=LIGN.IML
  402.         IMMT(IH)=MILIGN.IPNO(IM)
  403.   352 CONTINUE
  404. C
  405. C     Travail sur la partie inferieure
  406.       IF (ILIGNS.EQ.ILIGNM) THEN
  407.         MILIGN=IILIGN
  408.         ILIGNS=ILIGNL
  409.         LLIGN =LLIGL
  410.         LIGN  =LIGL
  411.         GOTO 39
  412.       ENDIF
  413. C
  414.       IMAX=IL1-1
  415.       IF (LSGDES) THEN
  416.         DO IT=IMIN,IL1-1
  417.           LIG1=ILIGNM.ILIGN(IT)
  418.           SEGACT /ERR=333/ LIG1
  419.           LIG1=ILIGNL.ILIGN(IT)
  420.           SEGACT /ERR=333/ LIG1
  421.         ENDDO
  422.   333   CONTINUE
  423.         IMAX=IT-1
  424.       ENDIF
  425. C
  426.       nbthrp=nbthr
  427.       nbthr=min(nbthrp,il2-il1+1)
  428.       nbthr=nbthro2 
  429.       IPER=IMIN
  430.       IDER=IMAX
  431.       IL1P=IL1
  432.       DO 13 I=IL1,IL2
  433. C
  434.   336   CONTINUE
  435. C
  436. C       Travail sur la partie inferieure
  437.         LILIGN=IILIGN
  438.         ILIGNS=ILIGNM
  439.         DO ITH=1,NBTHR-1
  440.           CALL THREADID(ITH,SHOLE3I)
  441.         ENDDO
  442.         CALL SHOLE3I(NBTHR)
  443.         DO ITH=NBTHR-1,1,-1
  444.           CALL THREADIF(ITH)
  445.         ENDDO
  446. C
  447. C       Travail sur la partie superieure
  448.         LILIGN=IILIGS
  449.         ILIGNS=ILIGNL
  450.         DO ITH=1,NBTHR-1
  451.           CALL THREADID(ITH,SHOLE3I)
  452.         ENDDO
  453.         CALL SHOLE3I(NBTHR)
  454.         DO ITH=NBTHR-1,1,-1
  455.           CALL THREADIF(ITH)
  456.         ENDDO
  457. C
  458.         IF (I.EQ.IL1P.AND.IDER.NE.IL1P-1) THEN
  459.           DO IT=IDER,IPER,-1
  460.             LIG1=ILIGNM.ILIGN(IT)
  461.             SEGDES LIG1
  462.             LIG1=ILIGNL.ILIGN(IT)
  463.             SEGDES LIG1
  464.           ENDDO
  465.           IPER=IDER+1
  466.           IDER=IL1-1
  467.           DO IT=IPER,IDER
  468.             LIG1=ILIGNM.ILIGN(IT)
  469.             SEGACT /ERR=335/ LIG1
  470.             LIG1=ILIGNL.ILIGN(IT)
  471.             SEGACT /ERR=335/ LIG1
  472.           ENDDO
  473.   335     CONTINUE
  474.           IDER=IT-1
  475.           IF (IDER.LT.IPER) THEN
  476.             CALL ERREUR(48)
  477.             RETURN
  478.           ENDIF
  479.           GOTO 336
  480.         ENDIF
  481. C
  482.         NBTHR=1
  483. C
  484. C       Travail sur la partie superieure
  485.         MILIGN=IILIGS
  486.         ILIGNS=ILIGNM
  487. C
  488.   444   CONTINUE
  489. C
  490.         LLIGN=MILIGN.ILIGN(I)
  491.         LIGN =ILIGNS.ILIGN(I)
  492.         NA=LIGN.IDERL-LIGN.IPREL+1
  493.         CALL SOMPAC(IPPVV(1),IMASQ(1),NA,KIVPO(1),KIVLO(1),NBPAR,IZROSF)
  494. C
  495.         NVALL=KIVLO(NBPAR)-1
  496.         LMASQ=MASQA(LLIGN.NJMAX)+1
  497. C
  498. C       LIGN est deja dimensionne a NJMAX
  499.         IF (ngmpet) THEN
  500.           SEGADJ,LIGN
  501.           LIG5=LIGN
  502.           GOTO 17
  503.         ENDIF
  504. C
  505.         SEGINI /ERR=109/ LIG5
  506.         GOTO 108
  507.  109    CONTINUE
  508.         IF (LSGDES) GOTO 14
  509. C       write(6,*) 'desactivation-3 ',1,il1p-1
  510.         lsgdes=.true.
  511.         do it=il1p-1,1,-1
  512.           lig1=ILIGNS.ILIGN(it)
  513.           segdes lig1
  514.           lig1=MILIGN.ilign(it)
  515.           segdes lig1
  516.         enddo
  517.         SEGINI /ERR=14/ LIG5
  518.  108    CONTINUE
  519.         LIG5.IML  =LIGN.IML
  520.         LIG5.IMM  =NBPAR
  521.         LIG5.IPREL=LIGN.IPREL
  522.         LIG5.IDERL=LIGN.IDERL
  523.         DO 111 LHE=1,NA
  524.           LIG5.IMMM(LHE)=LIGN.IMMM(LHE)
  525.  111    CONTINUE
  526.         DO 112 LHF=1,NA+1
  527.           LIG5.IPPVV(LHF)=LIGN.IPPVV(LHF)
  528.  112    CONTINUE
  529. C
  530.  17     CONTINUE
  531. C
  532.         DO 113 LHG=1,NBPAR
  533.           LIG5.IVPO(2*LHG-1)=KIVPO(LHG)
  534.           LIG5.IVPO(2*LHG)  =KIVLO(LHG)
  535.  113    CONTINUE
  536. C
  537.         DO IT=NBPAR-1,1,-1
  538.           IVI=LIG5.IVPO(2*IT)
  539.           IVF=LIG5.IVPO(2*(IT+1))-1
  540.           ICI=LIG5.IVPO(2*IT-1)
  541.           DO IC=MASQA(IVF-IVI+ICI),MASQA(ICI),-1
  542.             LIG5.IMASQ(IC)=IT
  543.           ENDDO
  544.         ENDDO
  545. C
  546.         IPA=1
  547.         DO 119 JPA=1,NA
  548.           IPP=LLIGN.IPPO(JPA)
  549.           KPA=LLIGN.IPPO(JPA+1)-IPP
  550.           LPA=LLIGN.LDEB(JPA)
  551.           LPA1=LPA-IPA
  552.           DO 120 MPA=1,KPA
  553.             XXV=LLIGN.XXVA(MPA+IPP)
  554.             IF (ABS(XXV).LE.XPETIT) GOTO 120
  555.             IPLA=LLIGN.LINC(MPA+IPP)-LPA1
  556.             IMSQ=LIG5.IMASQ(MASQA(IPLA))
  557.             IF (IMSQ.GT.0) THEN
  558. 114           CONTINUE
  559.               IF (IPLA.GE.LIG5.IVPO(2*(IMSQ+1)-1)) THEN
  560.                 IMSQ=IMSQ+1
  561.                 GOTO 114
  562.               ENDIF
  563.               IDBC=LIG5.IVPO(2*IMSQ-1)
  564.               IDBV=LIG5.IVPO(2*IMSQ  )
  565.               IPOS=IDBV+(IPLA-IDBC)
  566.               LIG5.VAL(IPOS)=XXV
  567.             ENDIF
  568. 120       CONTINUE
  569.           IPA=IPA+LLIGN.IMMMM(JPA)-LPA+1
  570. 119     CONTINUE
  571. C
  572.         ITR=NBPAR
  573.         DO 125 IPL=LMASQ,1,-1
  574.           IF (LIG5.IMASQ(IPL).GT.0) THEN
  575.             ITR=LIG5.IMASQ(IPL)
  576.           ELSE
  577.             LIG5.IMASQ(IPL)=-ITR
  578.           ENDIF
  579.  125    CONTINUE
  580. C
  581.         SEGSUP,LLIGN
  582.         IF (.not.ngmpet) SEGSUP,LIGN
  583.         MILIGN.ILIGN(I)=LIG5
  584. C
  585.         LMASQ=0
  586.         NVALL=0
  587.         SEGINI,LIG3
  588.         LIG3.IML=LIG5.IML
  589.         LIG3.IMM=LIG5.IMM
  590.         LIG3.IPREL=LIG5.IPREL
  591.         LIG3.IDERL=LIG5.IDERL
  592.         DO IPV=1,NA
  593.           LIG3.IMMM(IPV)=LIG5.IMMM(IPV)
  594.         ENDDO
  595.         DO IPV=1,NA+1
  596.           LIG3.IPPVV(IPV)=LIG5.IPPVV(IPV)
  597.         ENDDO
  598.         DO IPV=1,2*NBPAR
  599.           LIG3.IVPO(IPV)=LIG5.IVPO(IPV)
  600.         ENDDO
  601.         ILIGNS.ILIGN(I)=LIG3
  602. C
  603. C       Travail sur la partie inferieure
  604.         IF (ILIGNS.EQ.ILIGNM) THEN
  605.           MILIGN=IILIGN
  606.           ILIGNS=ILIGNL
  607.           GOTO 444
  608.         ENDIF
  609. C
  610.         IDER=I
  611.         IPER=IDER
  612.         IL1=I+1
  613.   13  CONTINUE
  614.       GOTO 15
  615.   14  CONTINUE
  616.       LIGN=ILIGNM.ILIGN(I)
  617.       SEGSUP,LIGN
  618.       LIGN=ILIGNL.ILIGN(I)
  619.       SEGSUP,LIGN
  620.       IL2=I-1
  621.   15  CONTINUE
  622.       IL1=IL1P
  623. C
  624.       if (lsgdes) then
  625. C       write(6,*) 'desactivation-4 ',1,il1-1
  626.         do it=il1-1,imin,-1
  627.           lig1=ILIGNM.ILIGN(it)
  628.           if (lig1.ne.0) segdes lig1
  629.           lig1=ILIGNL.ILIGN(it)
  630.           if (lig1.ne.0) segdes lig1
  631.         enddo
  632.       endif
  633.       nbthr=nbthrp
  634. C
  635. C  ****  BOUCLE  *5*  TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  636. C
  637.       ipos=0
  638.       iper=imin
  639.       ider=imin-1
  640.       DO 5 I=IMIN,IL2
  641.  
  642.       LIGM=LILIGN.ILIGN(I)
  643.       LIGL=MILIGN.ILIGN(I)
  644.       IF(I.LT.IL1) GOTO 7
  645. C
  646. C  ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  647. C  ******* IPREL  IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  648. C  ******* LE TERME DIAGONAL
  649. C
  650.       LIGN=LIGM
  651.       NBCOL=LIGN.IVPO(2*LIGN.IPPVV(2)-1)-1
  652.       DO 156 KHG=1,LIGN.IMMM(/1)
  653. C
  654. C       Travail sur la partie superieure
  655.         LIGN=LIGM
  656.         II=LIGN.IPREL-1+KHG
  657.         LIGN.IMMM(KHG)=0
  658. C
  659.         MM=LIGN.IVPO(2*LIGN.IPPVV(KHG))
  660.         NN=LIGN.IVPO(2*LIGN.IPPVV(KHG+1))
  661.         N=NN-MM
  662.         NN=NN-1
  663. C
  664.         DIAG(II)=LIGN.VAL(NN)
  665.         IF (N.EQ.1) THEN
  666.           IF (II-NBNNMC.GT.0) THEN
  667.             RE(II-NBNNMC,II-NBNNMC,1)=LIGN.VAL(NN)
  668.             GOTO 41
  669.           ENDIF
  670.           GOTO 8
  671.         ENDIF
  672. C
  673.         LIGN.VAL(NN)=LIGN.VAL(NN)+
  674.      &               CHOLE5(LIGN,LIGL,KHG,DIAG(1),NBCOL,0,NBOP(1))
  675.    8    CONTINUE
  676. C
  677. C       Travail sur la partie inferieure + terme diagonal
  678.         LIGN=LIGL
  679.         II=LIGN.IPREL-1+KHG
  680.         LIGN.IMMM(KHG)=0
  681. C
  682.         MM=LIGN.IVPO(2*LIGN.IPPVV(KHG))
  683.         NN=LIGN.IVPO(2*LIGN.IPPVV(KHG+1))
  684.         N=NN-MM
  685.         NN=NN-1
  686. C
  687.         LIGN.VAL(NN)=LIGN.VAL(NN)+
  688.      &               CHOLE5(LIGN,LIGM,KHG,DIAG(1),NBCOL,1,NBOP(1))
  689.         IF (II-NBNNMC.GT.0) THEN
  690.           RE(II-NBNNMC,II-NBNNMC,1)=LIGN.VAL(NN)
  691.           GOTO 41
  692.         ENDIF
  693. C
  694. C       Terme diagonal ok?
  695.         diagref=max(abs(diag(ii)),diagmin)
  696.         diagcmp=diagref*5d-12
  697.         IF (ABS(LIGN.VAL(NN)).GT.diagcmp) GOTO 12
  698. C
  699. C  il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement
  700. C  sur le calcul des reactions en cas de 2 relations presque identique
  701. C
  702. C  **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  703. C  **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  704. C  **** AU PREALABLE SUR CETTE INCONNUE.
  705. C
  706. *     write (6,*) ' ldmts mode d ensemble ittr ligne ',
  707. *    >      ittr(ii),ii,diag(ii),val(nn),diagref
  708. C  on garde le signe car il fau un moins sur les ML
  709.         vmaxi=diagref
  710.         do ipv=MM,NN
  711.           vmaxi=max(vmaxi,abs(lign.val(ipv)))
  712.         enddo
  713.         if(ittr(ii).NE.0) then
  714.           lign.VAL(NN)=lign.val(nn)-1.30D0*diagref
  715.           NENSLX=NENSLX+1
  716.         else
  717.           lign.val(nn)=vmaxi
  718.         endif
  719.         NENS=NENS+1
  720.         LIGL.IMMM(KHG)=NENS
  721.         LIGM.IMMM(KHG)=NENS
  722.    12   CONTINUE
  723. C
  724. C       Stabilisation
  725.         IF (ISTAB.NE.0) THEN
  726. C
  727. C         Elimination des Nan
  728.           if(.not.(abs(LIGN.val(nn)).lt.xgrand*xzprec).and.pasfait) then
  729.             pasfait=.false.
  730.             LIGN.val(nn)=xgrand*xzprec
  731.             write (6,*) 'Nan dans ldmts ligne',ii
  732.           endif
  733. C
  734.           diagcmp=abs(diagmin)*1d-5+xpetit/xzprec
  735.           IF (ITTR(II).EQ.0) THEN
  736.             IF (LIGN.VAL(NN).LT.-DIAGMAX*1D-3) THEN
  737.               LIGN.VAL(NN)=ABS(LIGN.VAL(NN))
  738. ***           val(nn)=diagmax*1d-3
  739.               IF (LIGN.IMMM(KHG).EQ.0) NENS=NENS+1
  740.               LIGN.IMMM(KHG)=NENS
  741.             ELSEIF (LIGN.VAL(NN).LE.DIAGCMP) THEN
  742.               IF (ISR.EQ.0.OR.IIMPI.GT.0)
  743.      &        write (ioimp,*) 'stabilisation RESO',ii,val(nn),diag(ii)
  744.               isr=isr+1
  745.               LIGN.VAL(NN)=MAX(DIAGCMP,-LIGN.VAL(NN))
  746.               IF (LIGN.IMMM(KHG).EQ.0) NENS=NENS+1
  747.               LIGN.IMMM(KHG)=NENS
  748.             ELSE
  749.               DIAGMIN=MIN(DIAGMIN,ABS(LIGN.VAL(NN)))
  750.             ENDIF
  751.           ELSE IF (ITTR(II).NE.0) THEN
  752.             IF (LIGN.VAL(NN).GE.ABS(DIAG(II))*STMULT) THEN
  753.               IF (ISRL.EQ.0.OR.IIMPI.GT.0)
  754.      &        write (ioimp,*) 'stabilisation RESO lagrange',ii,val(nn)
  755.               isrl=isrl+1
  756.               LIGN.VAL(NN)=-ABS(DIAG(II))*STMULT
  757.               IF (LIGN.IMMM(KHG).EQ.0) THEN
  758.                 NENS=NENS+1
  759.                 NENSLX=NENSLX+1
  760.               ENDIF
  761.               LIGN.IMMM(KHG)=NENS
  762.             ENDIF
  763.           ENDIF
  764.         ENDIF
  765. C
  766.         DIAG(II)=LIGN.VAL(NN)
  767.         IF(ABS(DIAG(II)).LE.XPETIT) THEN
  768.           DIAG(II)=diagmax
  769.           IF(ITTR(II).NE.0) DIAG(II)=-DIAGMAX
  770.           LIGL.VAL(NN)=DIAG(II)
  771.         ENDIF
  772.         LIGM.VAL(NN)=DIAG(II)
  773.         GOTO 41
  774. C
  775. C  ****  ERREUR MATRICE SINGUIERE
  776. C
  777.         INTERR(1)=I
  778.         CALL ERREUR(49)
  779.         if (ithrd.eq.1) then
  780.           call threadis
  781.          call oooprl(0)
  782.         endif
  783.         call ooomru(0)
  784.         RETURN
  785. C
  786. C  ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  787. C        ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  788. C
  789.    41   CONTINUE
  790.         IF(DIAG(II).LT.0.D0) INEG=INEG+1
  791.         IF(ITTR(II).NE.0)   NBLAG=NBLAG+1
  792. C
  793.         CONDMAX=MAX(CONDMAX,ABS(DIAG(II)))
  794.         IF (II.LE.NBNNMC) THEN
  795.           CONDMIN=MIN(CONDMIN,ABS(DIAG(II)))
  796.           DIAG(II)=1.D0/DIAG(II)
  797.         ENDIF
  798. C
  799.   156 CONTINUE
  800. C
  801. C     Suppression du masque
  802.       LMASQ=0
  803. C
  804.       LIGN=LIGL
  805.       NVALL=LIGN.VAL(/1)
  806.       NBPAR=LIGN.IVPO(/1)/2
  807.       NA   =LIGN.IMMM(/1)
  808.       SEGADJ,LIGN
  809. C
  810.       LIGN=LIGM
  811.       NVALL=LIGN.VAL(/1)
  812.       NBPAR=LIGN.IVPO(/1)/2
  813.       NA   =LIGN.IMMM(/1)
  814.       SEGADJ,LIGN
  815. C
  816.       IF (I.EQ.IDBGL) THEN
  817.         LIG1=LILIGN.ILIGN(I)
  818.         WRITE(*,*) 'LIGNE FINALE SUP',I,LIG1
  819.         segprt,lig1
  820.         LIG1=MILIGN.ILIGN(I)
  821.         WRITE(*,*) 'LIGNE FINALE INF',I,LIG1
  822.         segprt,lig1
  823.       ENDIF
  824. C
  825.       NVSTOR=NVSTOR+NVALL
  826.       NVSTRM=MAX(NVSTRM,NVALL*4)
  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=ILIGNL.ILIGN(ipv)
  1078.         segsup lig3
  1079.         lig3=ILIGNM.ILIGN(ipv)
  1080.         segsup lig3
  1081.       enddo
  1082.       SEGSUP,ILIGNL,ILIGNM
  1083.       segsup immt
  1084. C     write (6,*) ' ldmts ipos max ',iposm
  1085.  9999 continue
  1086.       if (ithrd.eq.1) then
  1087.         call threadis
  1088.        call oooprl(0)
  1089.       endif
  1090.       if(iimpi.ne.0)  write (6,*) 
  1091.      >      ' ldmts condmin condmax ',condmin,condmax,diag(/1)
  1092.       SEGDES,MDIAG
  1093.       RETURN
  1094.       END
  1095.  
  1096.  

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