Télécharger la plaquette Cast3M


fr en

Télécharger ldmt3.eso

Retour à la liste

Numérotation des lignes : 

ldmt3

  1. C LDMT3     SOURCE    PV090527  24/04/13    21:15:03     11827          
  2.       SUBROUTINE LDMT3(MMATRX,PREC)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8(A-H,O-Z)
  5. C  TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY
  6.       PARAMETER (LPCRAY=10000)
  7.       INTEGER OOOVAL,OOOLEN
  8.       dimension ittime(4)
  9.       POINTEUR LIG2.LIGN, LIG3.LIGN
  10.       POINTEUR L.LLIGN,   M.LLIGN
  11.       POINTEUR LL.MILIGN, MM.MILIGN, LILIGN.MILIGN
  12.       POINTEUR LLL.LIGN,   MMM.LIGN
  13.       SEGMENT ITEMP
  14.          REAL*8 P(INC)
  15.       ENDSEGMENT
  16. C      POINTEUR R.ITEMP,W.ITEMP
  17. C
  18. C  ****  MISE SOUS FORME A=L D Mt  DE LA MATRICE MMATRX
  19. C
  20.  
  21. -INC PPARAM
  22. -INC CCOPTIO
  23. -INC CCREEL
  24. -INC SMMATRI
  25.  
  26. -INC CCASSIS
  27. -INC CCHOLE
  28.       SEGMENT KIVPO(IIMAX)
  29.       SEGMENT KIVLO(IIMAX)
  30.       segment immt(nblig)
  31.       segment ireser(nvstrm)
  32.       external chole3i
  33.       SAVE IPASV
  34.       DATA IPASV/0/
  35. *  ngmpet dit si on tient en memoire (false) ou si on deborde (true)
  36.       logical ngmpet
  37. C     character*8 zen
  38. C     equivalence (zen,izen)
  39.       logical lsgdes,pasfait,ngdyn
  40. **    xkpar=0
  41. **    xkseq=0
  42.       ireser=0
  43.       matric=0
  44.       maitre=0
  45.       pasfait=.true.
  46.       lsgdes=.false.
  47. *  faire attention a respecter l'ordre des segdes par la suite
  48.       call ooomru(1)
  49.       condmax=0.d0
  50.       condmin=xgrand
  51.       ngmpet=.false.
  52.       ngdyn=.true.
  53.       call timespv(ittime,oothrd)
  54.       kcour=(ittime(1)+ittime(2))/10
  55.       kcourp=kcour
  56.       kcouri=kcour
  57.       kdiff=0
  58.       kcour=0
  59.       perf=0.d0
  60.       perfp=-1
  61.       nbchan=1
  62.       nbopit=0
  63.       iposm=0
  64. C     zen='CPU'//char(0)
  65. C     le=4
  66.       nvaor=0
  67.       nbthro=nbthrs
  68.       ithrd=0
  69.       if (nbthro.gt.1) then
  70.        ithrd=1
  71.        call threadii
  72.        call oooprl(1)
  73.       endif
  74.       nbthr=nbthro
  75.       do ith=1,nbthr
  76.        nbop(ith)=0
  77.       enddo
  78.       stmult=1d-5
  79.  
  80. C  nouvelle methode de gestion de l'espace memoire necessitee par la parallelisation
  81. C  memoire vive totale
  82.       MACTIT=OOOVAL(1,1)
  83. **    write(6,*) ' mactit igrand ',mactit,igrand
  84. C  un bloc de memoire fera au plus macti/2
  85.       call intpdo(inpdo)
  86.       nvstrm=mactit/10
  87.       MMATRI=MMATRX
  88.       SEGACT,MMATRI*MOD
  89.       PRCHLV=PREC
  90.       LL    =IILIGN
  91.       MM    =IILIGS
  92.       SEGACT,    MM*MOD,LL*MOD
  93.       INO=MM.ILIGN(/1)
  94.       MDIAG=IDIAG
  95.       SEGACT,MDIAG*MOD
  96.       NBLIG=INO
  97.       NBLIGI=INO
  98.       segini immt
  99.       precc=prec
  100.       INC=DIAG(/1)
  101.       nbnnmc=inc+1
  102.       nvstrm=max(inc*inpdo,nvstrm)
  103. **      write(6,*) ' nvstrm ',nvstrm
  104.       INCC=INC
  105.       MIMIK=IIMIK
  106.       MINCPO=IINCPO
  107.       SEGACT,MINCPO,MIMIK
  108.       IPLUMI=IMIK(/2)*2 +4
  109.       IL2=0
  110.       IIMAX=IJMAX+IPLUMI
  111.       SEGINI KIVPO,KIVLO
  112.       INEG=0
  113.       NBLAG=0
  114.       NENSLX=0
  115.       NVSTOC=0
  116.       NVSTOR=0
  117.       NVSTIC=0
  118.       NVSTIR=0
  119.       diagmax=XPETIT/XZPREC
  120.       diagmin=xgrand
  121.       do i=1,diag(/1)
  122.        if (ll.ittr(i).eq.0) diagmax=max(diagmax,abs(diag(i)))
  123.        if (ll.ittr(i).eq.0.and.abs(diag(i)).gt.xpetit/xzprec)
  124.      >   diagmin=min(diagmin,abs(diag(i)))
  125.       enddo
  126.       if (diagmax.le.xpetit/xzprec) then
  127.          do i=1,diag(/1)
  128.            diagmax=max(diagmax,abs(diag(i)))
  129.            if (abs(diag(i)).gt.xpetit/xzprec)
  130.      >       diagmin=min(diagmin,abs(diag(i)))
  131.          enddo
  132.       endif
  133.       diagmin=min(diagmin,diagmax)
  134. ***    write (6,*) ' ldmt3 diagmin diagmax ',diagmin,diagmax,diag(/1)
  135. C
  136.  
  137. C
  138. C
  139. C  **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  140. C  **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  141. C
  142. C  ****  LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  143. C
  144.   1   CONTINUE
  145.       IVALMA=IJMAX+IPLUMI
  146.       IL1=IL2+1
  147.       IVALMI=IJMAX+IPLUMI
  148.       IMINM=IL1
  149.       IMINL=IL1
  150. *  reserver de la place ou mettre les lignes superieures dans le cas debordement
  151.       if (ngmpet) then
  152.        if(ireser.eq.0) segini ireser
  153.       endif
  154.       DO 2 I=IL1,INO
  155.       ngdyn=.true.
  156.       m=0
  157.       l=0
  158.       lll=0
  159.       mmm=0
  160.  
  161.       M    = MM.ILIGN(I)
  162.       L    = LL.ILIGN(I)
  163.       SEGACT /ERR=32/M
  164.       SEGACT /ERR=32/L
  165.       goto 31
  166.   32  continue
  167. **       write(6,*) ' segact llign erreur',i,il1,lsgdes
  168.       if (.not.lsgdes) then
  169. **         write(6,*) ' lsgdes 1 '
  170.         lsgdes=.true.
  171. ****    ngmpet=.true.
  172. **      write(6,*) 'desactivation-1 ',1,il1-1
  173.         do it=il1-1,1,-1
  174.          lign=mm.ilign(it)
  175.          segdes lign
  176.          lign=ll.ilign(it)
  177.          segdes lign
  178.         enddo
  179.       else
  180.         goto 3
  181.       endif
  182.       SEGACT /ERR=3/M
  183.       SEGACT /ERR=3/L 
  184.   31  continue
  185.       NA=  M.IMMMM(/1)
  186. C*    write (6,*) ' chole ligne noeud inconnues ',i,ipno(i),na
  187.       NBPAR=NA+1
  188.       NVALL=M.NJMAX 
  189.       LMASQ=NVALL/MASDIM+2
  190.       lmasqm=lmasq
  191.       NVALLL=NVALL
  192.       mmm=0
  193.       lll=0
  194.       SEGINI /ERR=33/MMM
  195.  
  196.       NA=L.IMMMM(/1)
  197.       NBPAR=NA+1
  198.       NVALL= L.NJMAX
  199.       LMASQ=NVALL/MASDIM+2
  200.       lmasql=lmasq
  201.       NVILL=NVALL
  202.       lll=0
  203.       SEGINI /ERR=33/LLL
  204. C  recuperer la longueur du segment
  205.       lglig=na*(nvall/na)**1.33333333333333333333333333
  206.       goto 34
  207.   33  continue
  208. **       write(6,*) ' segini xx.llign erreur',i,il1
  209.       if (mmm.ne.0) segsup mmm
  210.       if (.not.lsgdes) then
  211. **        write(6,*) ' lsgdes 2 '
  212.         lsgdes=.true.
  213. ****    ngmpet=.true.
  214. **      write(6,*) 'desactivation-2 ',1,il1-1
  215.         do it=il1-1,1,-1
  216.          lign=mm.ilign(it)
  217.          segdes lign
  218.          lign=ll.ilign(it)
  219.          segdes lign
  220.         enddo
  221.       else
  222.         goto 3
  223.       endif
  224.       SEGINI /ERR=33/MMM
  225.       SEGINI /ERR=33/LLL
  226.   34  continue
  227. C  recuperer la longueur du segment
  228.       lglig=na*(nvall/na)**1.33333333333333333333333333
  229.  
  230.  
  231.       NVSTOC=NVSTOC + NVALLL
  232.       IVALMA=IVALMA + NVALLL
  233.       NVSTIC=NVSTIC + NVALL
  234.       IVALMI=IVALMI + NVALL
  235.       NVALL=NVALLL
  236.  
  237.       nvaor = nvaor + M.XXVA(/1)
  238.       nvaori= nvaori+ L.XXVA(/1)
  239. C
  240. C  ****  DECOMPACTAGE
  241. C
  242.       IPA=1
  243.       NA=M.IMMMM(/1)
  244.       DO 121 JPA=1,NA
  245.         MMM.IVPO(JPA)=IPA
  246.         KPA = M.IPPO(JPA+1)-M.IPPO(JPA)
  247.         IPP = M.IPPO(JPA)
  248.         MMM.IPPVV(JPA)=IPA-1
  249.         LPA  = M.LDEB(JPA)
  250.         LPA1 = LPA-IPA
  251.  
  252.         DO 122 MPA=1,KPA
  253.           LLO  = M.LINC(MPA+IPP)
  254.           IPLA = LLO -LPA1
  255.           xxv=m.xxva(mpa+ipp)
  256.           if(abs(xxv).gt.xpetit) then
  257.               MMM.VAL(IPLA)=xxv
  258.               MMM.IMASQ(IPLA/MASDIM+1)=1
  259.            if (ipla-ipa+1.ge.1) MMM.IMASQ((IPLA-ipa+1)/MASDIM+1)=1
  260.           endif
  261. 122     CONTINUE
  262.  
  263.         IPA=IPA+M.IMMMM(JPA)-LPA + 1
  264. Cpv        MMM.IMMM(JPA)=MM.IPNO(LPA)
  265.         MMM.IMMM(JPA)=LPA
  266.         IF(IMINM  .GT.MM.IPNO(LPA )) IMINM  = MM.IPNO(LPA)
  267.  
  268. 121   CONTINUE
  269.  
  270. * indexation de imasq
  271.         ipln=lmasq/na
  272.         iplp=lmasq/na
  273. **      write (6,*) 'ldmt3 271 lmasq imasq ',lmasq,mmm.imasq(/1)
  274.         do 123 ipl=lmasqm/na,1,-1
  275.           if (mmm.imasq(ipl).gt.0) then
  276.             mmm.imasq(ipl)=iplp*masdim
  277.             ipln=ipl-1
  278.           else
  279.             mmm.imasq(ipl)=-ipln*masdim
  280.             iplp=ipl-1
  281.           endif
  282.  123    continue
  283. **      write (6,*) ' imasq ',lmasq/na
  284. **      write (6,*) (imasq(ipl),ipl=1,lmasq/na)
  285.  
  286.       IPI=1
  287.       NA=L.IMMMM(/1)
  288.       DO 1210 JPA=1,NA
  289.         LLL.IVPO(JPA)=IPI
  290.         KPI =L.IPPO(JPA+1)-L.IPPO(JPA)
  291.         IPPI=L.IPPO(JPA)
  292.         LLL.IPPVV(JPA)=IPI-1
  293.         LPAI =L.LDEB(JPA)
  294.         LPA1I=LPAI-IPI
  295.  
  296.         DO 1220 MPA=1,KPI
  297.           LLO=L.LINC(MPA+IPPI)
  298.           IPLA=LLO-LPA1I
  299.           xxv=l.xxva(mpa+ippi)
  300.           if(abs(xxv).gt.xpetit) then
  301.              LLL.VAL(IPLA)=xxv
  302.              LLL.IMASQ(IPLA/MASDIM+1)=1
  303.            if (ipla-ipi+1.ge.1) LLL.IMASQ((IPLA-ipi+1)/MASDIM+1)=1
  304.           endif
  305. 1220    CONTINUE
  306.  
  307.         IPI=IPI+L.IMMMM(JPA)-LPAI+ 1
  308. Cpv        LLL.IMMM(JPA)=LL.IPNO(LPAI)
  309.            LLL.IMMM(JPA)=LPAI
  310.         IF(IMINL.GT.LL.IPNO(LPAI)) IMINL= LL.IPNO(LPAI)
  311. 1210   CONTINUE
  312. C*** **** ****
  313. * indexation de imasq
  314.         ipln=lmasq/na
  315.         iplp=lmasq/na
  316. **      write (6,*) 'ldmt3 314 lmasq imasq ',lmasq,lll.imasq(/1)
  317.         do 1230 ipl=lmasql/na,1,-1
  318.           if (lll.imasq(ipl).gt.0) then
  319.             lll.imasq(ipl)=iplp*masdim
  320.             ipln=ipl-1
  321.           else
  322.             lll.imasq(ipl)=-ipln*masdim
  323.             iplp=ipl-1
  324.           endif
  325.  1230    continue
  326. **      write (6,*) ' imasqi ',lmasq/na
  327. **      write (6,*) (imasqi(ipl),ipl=1,lmasq/na)
  328.  
  329.       NA=M.IMMMM(/1)
  330.       if (NA.gt.0) then
  331.         MMM.IPREL=M.IMMMM(1)
  332.         MMM.IDERL=M.IMMMM(NA)
  333.         mm.lcara(2,i)=mmm.iprel
  334.         mm.lcara(3,i)=mmm.iderl
  335.       endif
  336.       MMM.IPPVV(NA+1)=IPA-1
  337.  
  338.       NA=L.IMMMM(/1)
  339.       if (NA.gt.0) then
  340.         LLL.IPREL=L.IMMMM(1)
  341.         LLL.IDERL=L.IMMMM(NA)
  342.         ll.lcara(2,i)=lll.iprel
  343.         ll.lcara(3,i)=lll.iderl
  344.       endif
  345.       LLL.IPPVV(NA+1)=IPI-1
  346.  
  347.       SEGSUP L,M
  348.       LL.ILIGN(I)=LLL
  349.       MM.ILIGN(I)=MMM
  350. C*    write (6,*) 'longueur ligne ',nvall
  351. C nb de ligne multiple du nb de threads
  352. C blocage ligne lecture-ecriture pour minimiser le cache
  353. C on note si on est au minimum de lignes
  354.       nbthro=min(nbthrs,lglig/1200+1)
  355.       if (i+1-il1.ge.nbthro.and.(.not.ngmpet)) then
  356.        nbthro=min(nbthrs,i+1-il1)
  357.        ngdyn=.true.
  358.         if(i+1-il1.eq.nbthrs) ngdyn=.false.
  359.          il2=i
  360.          GOTO 4
  361.       endif
  362.  
  363.  
  364.     2 CONTINUE
  365.       IL2=INO
  366.       GO TO 4
  367.     3 IL2=I-1
  368.         if(m.ne.0) segdes m
  369.         if(l.ne.0) segdes l
  370.         if (lll.ne.0) segsup lll
  371.         if (mmm.ne.0) segsup mmm        
  372.     4 CONTINUE
  373.       nbthro=min(nbthrs,nbthro)
  374.       nbthr=nbthro
  375.       if(ireser.ne.0) segsup ireser
  376. C      WRITE(IOIMP,*) 'Mactic = ', mactic, nbthr
  377. C
  378.       IF(IL2.GE.IL1) GO TO 40
  379. C
  380. C ****  APPEL AUX ERREURS  MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  381. C
  382. C      ITYP=48
  383.       CALL ERREUR(48)
  384.         call ooodmp(0)
  385.       if (ithrd.eq.1) then
  386.        call threadis
  387.        call oooprl(0)
  388.       endif
  389.       call ooomru(0)
  390.       RETURN
  391.    40 CONTINUE
  392.       IM=INC
  393.       DO 352 IH=IL2,IL1,-1
  394.         MMM=MM.ILIGN(IH)
  395.         IL=INC
  396.         DO 354 JH=1,MMM.IMMM(/1)
  397.           IM=MIN(IM,MMM.IMMM(JH))
  398.           IL=MIN(IL,MMM.IMMM(JH))
  399.   354   CONTINUE
  400.  
  401.         MMM.IML=IL
  402.         mm.lcara(1,iH)=IL
  403.         MMM.IMM=MM.ipno(IM)
  404.         if (immt(ih).ne.0) then
  405.          immt(ih)=min(immt(ih),mm.ipno(IM))
  406.         else
  407.          immt(ih)=mm.ipno(IM)
  408.         endif
  409.   352 CONTINUE
  410. C  353 CONTINUE
  411.       MMM=MM.ILIGN(IL1)
  412.       IL11=MMM.IPREL
  413.  
  414.       IM=INC
  415.       DO 3520 IH=IL2,IL1,-1
  416.         LLL=LL.ILIGN(IH)
  417.         IL=INC
  418.         DO 3540 JH=1,LLL.IMMM(/1)
  419.           IM=MIN(IM,LLL.IMMM(JH))
  420.           IL=MIN(IL,LLL.IMMM(JH))
  421.  3540   CONTINUE
  422.  
  423.         LLL.IML=IL
  424.         ll.lcara(1,iH)=IL
  425.         LLL.IMM=LL.ipno(IM)
  426.         if (immt(ih).ne.0) then
  427.          immt(ih)=min(immt(ih),ll.ipno(IM))
  428.         else
  429.          immt(ih)=ll.ipno(IM)
  430.         endif
  431.  3520 CONTINUE
  432. C 3530 CONTINUE
  433.       LLL=LL.ILIGN(IL1)
  434.       IL22=LLL.IPREL
  435. C
  436. C  ****  BOUCLE  *5*  TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  437. C
  438. C      lig1=MM.ilign(IMINM)
  439. C      lig2=LL.ilign(IMINL)
  440.  
  441.       ipos=0
  442.       iper=IMINM
  443.       ider=IMINM-1
  444.       iderac=IMINM-1
  445.  
  446.  
  447.       IMINA=MIN(IMINM,IMINL)
  448.       IMIN = IMINA
  449.       DO 5 I=IMINA,IL2
  450.         LIG1=MM.ILIGN(I)
  451.         LIG2=LL.ILIGN(I)
  452.         IF(I.LT.IL1)  GO TO 7
  453. C
  454. C       ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  455. C       ******* IPREL  IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  456. C       ******* LE TERME DIAGONAL
  457. C
  458.  
  459. C       on s'occupe d'abord de la partie superieur
  460.         LIGN=LIG1
  461.         NAA= IMMM(/1)
  462.  
  463.         DO 156 KHG=1,NAA
  464.         LIGN=LIG1
  465.         II=IPREL-1+KHG
  466.         IMMM(KHG)=0
  467.         NN1=IPPVV(KHG+1)
  468.         NNM1=IPPVV(KHG)
  469.         NNM1S=NNM1
  470.         N=NN1-NNM1
  471.         DIAG(II)=VAL(NN1)
  472.         diagref=diag(ii)
  473.         IF(N.EQ.1)  GO TO 8
  474.         NMI=N-II
  475.         IDEP=MAX(IL11,2-NMI)
  476.         KIDEP=IDEP+NMI
  477.         KI1=N-1
  478.         KQ=-NMI
  479. C        WRITE(IOIMP,*)'Avant CHOLI1, Imasq(/1)=',Imasq(/1)
  480. *        WRITE(IOIMP,*)'Avant CHOLI1-1 ',val(nn1)
  481.         VAL(NN1)=VAL(NN1)+
  482.      #  CHOLI1(LL.ILIGN,LIG2,LIG1.VAL(1+IPPVV(KHG)),DIAG(1-NMI),
  483.      #   LL.IPNO(1-NMI),LIG1.IPPVV(1),KHG,LIG1.IVPO(1),KIDEP,KI1,
  484.      #   KQ,LIG1.imasq(1),1+IPPVV(KHG),PREC,1,nbop(1))
  485.          lig1.imasq(nn1/masdim+1)=1
  486.          lig1.imasq(n/masdim+1)=1
  487.    8    CONTINUE
  488.  
  489.         LIGN=LIG2
  490.         II=IPREL-1+KHG
  491.         IMMM(KHG)=0
  492.         NN2=IPPVV(KHG+1)
  493.         NNM1=IPPVV(KHG)
  494.         N=NN2-NNM1
  495.         IF(N.EQ.1)  GO TO 88
  496.         NMI=N-II
  497.         IDEP=MAX(IL22,2-NMI)
  498.         KIDEP=IDEP+NMI
  499.         KI1=N-1
  500.         KQ=-NMI
  501. *        WRITE(IOIMP,*)'Avant CHOLI1-2 ',val(nn2)
  502.         VAL(NN2)=VAL(NN2)+
  503.      #  CHOLI1(MM.ILIGN,LIG1,VAL(1+IPPVV(KHG)),DIAG(1-NMI),
  504.      #   MM.IPNO(1-NMI),IPPVV(1),KHG,IVPO(1),KIDEP,KI1,
  505.      #   KQ,LIG2.imasq(1),1+IPPVV(KHG),PREC,2,nbop(1))
  506.          lig2.IMASQ(nn2/masdim+1)=1
  507.          lig2.IMASQ(n/masdim+1)=1
  508.   88    CONTINUE
  509.         LIGN=LIG1
  510.       diagref=max(abs(diag(ii)),diagmin)
  511.       diagcmp=diagref*5d-12
  512.       IF(LL.ITTR(II).EQ.0.AND.
  513.      & ABS(LIG2.VAL(NN2)).GT.diagcmp)  GO TO 12
  514.       IF(LL.ITTR(II).NE.0.AND.
  515.      & ABS(LIG2.VAL(NN2)).GT.diagcmp)  GO TO 12
  516. C  il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement
  517. C  sur le calcul des reactions en cas de 2 relations presque identique
  518. C
  519. C  **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  520. C  **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  521. C  **** AU PREALABLE SUR CETTE INCONNUE.
  522. C
  523. *     write (6,*) ' ldmt3 mode d ensemble ittr ligne ',
  524. *    >      ittr(ii),ii,diag(ii),val(nn2),lig2.val(nn2),diagref,diagcmp
  525. C  on garde le signe car il fau un moins sur les ML
  526.        vmaxi=diagref
  527.        do ipv=1+ippvv(khg),nn2
  528.         vmaxi=max(vmaxi,abs(lig2.val(ipv)))
  529.        enddo
  530.         if(LL.ittr(ii).NE.0) then
  531.          LIG2.VAL(NN2)=LIG2.VAL(NN2)-4.D0*diagref
  532.          NENSLX=NENSLX+1
  533.         else
  534.          LIG2.VAL(NN2)=vmaxi
  535.         endif
  536.         NENS=NENS+1
  537.         LIG2.IMMM(KHG)=NENS
  538.         LIG1.IMMM(KHG)=NENS
  539.    12   CONTINUE
  540.  
  541.  
  542.         DIAG(II)=LIG2.VAL(NN2)
  543.         IF(DIAG(II).NE.0.D0)  GO TO 41
  544.  
  545.         KQ1=1+NNM1S
  546.         KQN=N+NNM1S
  547.         DO 16 LFG=KQ1,KQN
  548.           IF(LIG1.VAL(LFG).NE.0.D0)  GO TO 17
  549.    16   CONTINUE
  550.  
  551.         KQ1=1+NNM1
  552.         KQN=N+NNM1
  553.         DO 160 LFG=KQ1,KQN
  554.           IF(LIG2.VAL(LFG).NE.0.D0)  GO TO 170
  555.   160   CONTINUE
  556.  
  557.         DIAG(II)=1.D0
  558.         if (LL.ittr(ii).ne.0) diag(ii)=-1
  559.         LIG2.VAL(NN2)=DIAG(II)
  560.         GO TO 41
  561.    17   CONTINUE
  562. C       write (6,*) ' ldmt3 apres 17 ',val(lfg)
  563.         diag(ii)=LIG1.VAL(LFG)
  564.         goto 171
  565.   170  CONTINUE
  566.         diag(ii)=LIG2.VAL(LFG)
  567.   171  continue
  568.         if (LL.ittr(ii).ne.0) diag(ii)=-abs(diag(ii))
  569. ***     LIG2.val(nn2)=diag(ii)
  570.         GOTO 41
  571. C
  572. C       ****  ENVOI ERREUR MATRICE SINGUIERE
  573. C
  574. C      ITYP=49
  575.       INTERR(1)=I
  576.       CALL ERREUR(49)
  577.       if (ithrd.eq.1) then
  578.         call threadis
  579.        call oooprl(0)
  580.       endif
  581.       call ooomru(0)
  582.       RETURN
  583. C
  584. C       ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  585. C             ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  586. C
  587.    41   IF(DIAG(II).LT.0.D0)  INEG=INEG+1
  588.         IF(LL.ITTR(II).NE.0) NBLAG=NBLAG+1
  589.         LIG1.VAL(NN1)=DIAG(ii)
  590.        condmin=min(condmin,abs(diag(ii)))
  591.        condmax=max(condmax,abs(diag(ii)))
  592.        diag(ii)=1.d0/diag(ii)
  593.   156   CONTINUE
  594. C
  595. C       RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE)
  596. C       d'abord la triangulaire superieure
  597. C
  598.         NA=LIG1.IMMM(/1)
  599.         NBPAR=NA+1
  600.         if (na.gt.0)
  601.      >  CALL COMPAC(LIG1.VAL(1),NBPAR,KIVPO(1),KIVLO(1),
  602.      #      NVALL,LIG1.IPPVV(1),IZROSF,NA,PREC,lig1.imasq(1),
  603.      #      LIG1.IPREL,LIG1.IDERL)
  604.  
  605. C       on recree lig1 car la compacter en place emiette la memoire
  606.        lmasq=0
  607.        lig3=lig1
  608.        segini /err=1431/ lig3
  609.  1431  continue
  610. *  deplacement fait ici maintenant, avec unrolling
  611.       do 300 nbp=1,nbpar-1
  612.        kdif =kivpo(nbp)-kivlo(nbp)
  613.        do     iv=kivlo(nbp),kivlo(nbp+1)-4,4
  614.           lig3.val(iv)=lig1.val(iv+kdif )
  615.           lig3.val(iv+1)=lig1.val(iv+1+kdif )
  616.           lig3.val(iv+2)=lig1.val(iv+2+kdif )
  617.           lig3.val(iv+3)=lig1.val(iv+3+kdif )
  618.        enddo
  619.        do     iv1=iv,kivlo(nbp+1)-1
  620.           lig3.val(iv1)=lig1.val(iv1+kdif )
  621.        enddo
  622.  300  continue
  623. **     do it=1,nvall
  624. **      lig3.val(it)=lig1.val(it)
  625. **     enddo
  626.        do it=1,na
  627.         lig3.immm(it)=lig1.immm(it)
  628.         lig3.ippvv(it)=lig1.ippvv(it)
  629.        enddo
  630.         lig3.ippvv(na+1)=lig1.ippvv(na+1)
  631.         lig3.iml=lig1.iml
  632.         lig3.iprel=lig1.iprel
  633.         lig3.iderl=lig1.iderl
  634.         mm.lcara(1,i)=lig3.iml
  635.         mm.lcara(2,i)=lig3.iprel
  636.         mm.lcara(3,i)=lig3.iderl
  637.         if (lig3.ne.lig1) then
  638.           segsup lig1
  639.         else
  640.           segadj lig3
  641.         endif
  642.         lig1=lig3
  643.         mm.ilign(i)=lig3
  644.         NVSTOR=NVSTOR+NVALL
  645.         nvstrm=max(nvstrm,nvall)
  646.         DO 143 LHG=1,NBPAR
  647.          LIG1.IVPO(2*LHG-1)=KIVPO(LHG)
  648.          LIG1.IVPO(2*LHG)  =KIVLO(LHG)
  649.  143    CONTINUE
  650. C
  651. C       RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE)
  652. C       puis la triangulaire inférieure
  653. C
  654.         NA=LIG2.IMMM(/1)
  655.         NBPAR=NA+1
  656.         if (na.gt.0)
  657.      >  CALL COMPAC(LIG2.VAL(1),NBPAR,KIVPO(1),KIVLO(1),
  658.      #      NVILL,LIG2.IPPVV(1),IZROSF,NA,PREC,lig2.imasq(1),
  659.      #      LIG2.IPREL,LIG2.IDERL)
  660.  
  661.         NVALL=NVILL
  662. C       on recree lig2 car la compacter en place emiette la memoire
  663. C        WRITE(IOIMP,*) 'Valeur de LIG2', LIG2
  664.        lig3=lig2
  665.        lmasq=0
  666.        segini /err=1432/ lig3
  667.  1432  continue
  668. *  deplacement fait ici maintenant, avec unrolling
  669.       do 301 nbp=1,nbpar-1
  670.        kdif =kivpo(nbp)-kivlo(nbp)
  671.        do     iv=kivlo(nbp),kivlo(nbp+1)-4,4
  672.           lig3.val(iv)=lig2.val(iv+kdif )
  673.           lig3.val(iv+1)=lig2.val(iv+1+kdif )
  674.           lig3.val(iv+2)=lig2.val(iv+2+kdif )
  675.           lig3.val(iv+3)=lig2.val(iv+3+kdif )
  676.        enddo
  677.        do     iv1=iv,kivlo(nbp+1)-1
  678.           lig3.val(iv1)=lig2.val(iv1+kdif )
  679.        enddo
  680.  301  continue
  681. **       do it=1,nvall
  682. **        lig3.val(it)=lig2.val(it)
  683. **       enddo
  684.        do it=1,na
  685.         lig3.immm(it)=lig2.immm(it)
  686.         lig3.ippvv(it)=lig2.ippvv(it)
  687.        enddo
  688.         lig3.ippvv(na+1)=lig2.ippvv(na+1)
  689.         lig3.iml=lig2.iml
  690.         lig3.iprel=lig2.iprel
  691.         lig3.iderl=lig2.iderl
  692.         ll.lcara(1,i)=lig3.iml
  693.         ll.lcara(2,i)=lig3.iprel
  694.         ll.lcara(3,i)=lig3.iderl
  695.         if (lig3.ne.lig2) then
  696.           segsup lig2
  697.         else
  698.           segadj lig3
  699.         endif
  700.         lig2=lig3
  701.         ll.ilign(i)=lig2
  702.         NVSTIR=NVSTIR+NVILL
  703.         nvstrm=max(nvstrm,nvIll*inpdo)
  704.         DO 1430 LHG=1,NBPAR
  705.           LIG2.IVPO(2*LHG-1)=KIVPO(LHG)
  706.           LIG2.IVPO(2*LHG)  =KIVLO(LHG)
  707. 1430    CONTINUE
  708.  
  709.  
  710.  
  711.  
  712.  
  713.       IF (I.GT.1) THEN
  714.        LIG1=MM.ILIGN(I-1)
  715.        LIG2=LL.ILIGN(I-1)
  716.        if (lsgdes) SEGDES LIG1,LIG2
  717.        IDERAC=MIN(IDERAC,I-2)
  718.       ENDIF
  719.  
  720. C
  721. C       ****  ON TRIANGULARISE LES AUTRES LIGNES
  722. C
  723.  
  724.         IL1=IL1+1
  725.         IF (IL1.GT.IL2) GOTO 5
  726.         LIG1=MM.ILIGN(I)
  727.         LIGN=MM.ILIGN(IL1)
  728.         IL11=IPREL
  729.         LIG2=LL.ILIGN(I)
  730.         LIGN=LL.ILIGN(IL1)
  731.         IL22=IPREL
  732.         GOTO 7
  733. C  72    CONTINUE
  734.   71    CONTINUE      
  735. *  passage en superlent
  736.         if (ider.lt.il1-1.and..not.ngmpet) then
  737.           ngmpet=.true.
  738.         endif
  739.       if (iper.gt.ider) then
  740.         call erreur(48)
  741.         call ooodmp(0)
  742.       if (ithrd.eq.1) then
  743.        call threadis
  744.        call oooprl(0)
  745.       endif
  746.       call ooomru(0)
  747.         return
  748.       endif
  749.  
  750.  
  751. ***       IF (I.LT.IL1-10) THEN
  752. C        SOIT PARCE QU'ON A FINI, SOIT PARCE QU'ON MANQUE DE MEMOIRE
  753. C        IL FAUT EXECUTER LES LIGNES ACTIVEES PUIS LES DESACTIVER
  754. C        LANCER LES CHOLE3 ET ATTENDRE QU'ILS SOIENT FINIS
  755.         IF (IPOS.NE.0) THEN
  756. C         WRITE (6,*) ' LANCEMENT THREAD ',IPER,IDER,IL1,IL2
  757.           IF (IPER.GT.IDER) THEN
  758.             WRITE (6,*) ' ERREUR INTERNE LDMT3 '
  759.             CALL ERREUR(5)
  760.           ENDIF
  761. C          WRITE (6,*) ' NBTHR-2 ',NBTHR
  762.  
  763.  
  764.           NBTHR=MIN(NBTHR,IL2-IL1+1)
  765. **         WRITE (6,*) ' NBTHR-3 ',NBTHR
  766.             MILIGN=MM
  767.             LILIGN=LL
  768. C            Write(6,*) 'ldmt3.eso On passe LL dans LILIGN : ', LILIGN
  769. *  blocage pour rester dans le cache secondaire
  770.         ipers=iper
  771.         iders=ider
  772.         ipas=1500
  773.         if(nbthr.eq.1) ipas=igrand
  774.  401    continue
  775.         ider=min(iders,iper+ipas-1)
  776.             DO ITH=1,NBTHR-1
  777.               CALL THREADID(ITH,CHOLE3I)
  778.             ENDDO
  779.             CALL CHOLE3I(NBTHR)
  780.             DO ITH=1,NBTHR-1
  781.               CALL THREADIF(ITH)
  782.             ENDDO
  783.          iper=iper+ipas
  784.          ipas=ipas/2
  785.          ipas=max(ipas,750)
  786.          if(iper.le.iders) goto 401
  787.             MILIGN=LL
  788.             LILIGN=MM
  789. C            Write(6,*) 'ldmt3.eso On passe MM dans LILIGN : ', LILIGN
  790. *  blocage pour rester dans le cache secondaire
  791.         ipas=1500
  792.         if(nbthr.eq.1) ipas=igrand
  793.         iper=ipers
  794.  402    continue
  795.         ider=min(iders,iper+ipas-1)
  796.             DO ITH=1,NBTHR-1
  797.               CALL THREADID(ITH,CHOLE3I)
  798.             ENDDO
  799.             CALL CHOLE3I(NBTHR)
  800.             DO ITH=1,NBTHR-1
  801.               CALL THREADIF(ITH)
  802.             ENDDO
  803.          iper=iper+ipas
  804.          ipas=ipas/2
  805.          ipas=max(ipas,750)
  806.          if(iper.le.iders) goto 402
  807.         iper=ipers
  808.         ider=iders
  809.         ENDIF
  810. *  test ctrlC
  811.       if (ierr.ne.0) goto 9999
  812.  
  813.         IPOSM=MAX(IPOSM,IPOS)
  814.         IPOS=0
  815.         IDERF=IDER-1
  816.         IDERF=IDER
  817.         IF (IDER.NE.IL1-1) IDERF=IDER
  818.         if(lsgdes) then
  819.         DO IL=IDERF,IPER,-1
  820.           LIGN=MM.ILIGN(IL)
  821.           SEGDES LIGN
  822.           LIGN=LL.ILIGN(IL)
  823.           SEGDES LIGN
  824. C         WRITE (6,*) ' DESACTIVATION LIGNE & COLONNE ',IL
  825.         ENDDO
  826.         endif
  827.         IDERAC=MIN(IDERAC,IPER-1)
  828.         IPER=IDER+1
  829. C       WRITE (6,*) ' IPER IDER IL1 ',IPER,IDER,IL1
  830.         IF (IPER.NE.IL1) GOTO 7
  831.         GOTO 5
  832.     7   CONTINUE
  833.  
  834.  
  835.         if(lsgdes) then 
  836.         SEGACT/ERR=71/LIG1
  837.         SEGACT/ERR=71/LIG2
  838.         endif
  839.         IPOS=IPOS+1
  840.         IDER=I
  841.         IF (I.GT.IDERAC) IDERAC=I
  842.         IF (I.EQ.IL1-1)   GOTO 71
  843.     5 CONTINUE
  844.       if(lsgdes) then
  845.       DO  I=min(IL1,IL2),max(il1,il2)
  846.         LLL=LL.ILIGN(I)
  847.         if (lll.ne.0) segdes lll
  848.         MMM=MM.ILIGN(I)
  849.         if (mmm.ne.0) segdes mmm                
  850. C        Write(6,*)'     SEGDES de LLL, MMM : ',LLL,MMM
  851.       enddo
  852.       endif
  853.       nbopt=0
  854.       do ith=1,nbthro
  855.        nbopt=nbopt+nbop(ith)
  856.        nbop(ith)=0
  857.       enddo
  858.       nbopin=nbopt
  859.       nbopit=nbopit+nbopin
  860.       call timespv(ittime,oothrd)
  861.       kcour=(ittime(1)+ittime(2))/10
  862.       kdiff=kcour-kcourp
  863. C*      write (6,*) ' nb operation temps ',nbopin,kdiff
  864.       if (kdiff.ge.1) then
  865.        perf=real(nbopin)/kdiff
  866. C*    if (nbchan.ne.0) perfp=perf
  867.        if (ngdyn) then
  868.         if (perf.lt.perfp*0.90 .and.nbchan.ne.1  ) then
  869.          nbchan=1
  870.          perfp=perf
  871.         elseif (nbchan.eq.0) then
  872.          nbchan=-1
  873.          perfp=max(perf,perfp)
  874.         else
  875.          nbchan=0
  876.         endif
  877.        endif
  878. C*    nbchan=0
  879.       endif
  880.       kcourp=kcour
  881.  
  882.       iderac=min(iderac,il1-1)
  883.       if(ireser.ne.0) segsup ireser
  884.       IF(IL2.LT.INO)  GO TO 1
  885. C     ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  886. C     ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  887. C       INEG=INEG-NBLAG
  888. C     on ne compte pas 2 fois les multiplicateurs qui vont etre
  889. C     elimines lors de la resolution car mode d'ensemble
  890.       INEG=INEG-(NBLAG-NENSLX)
  891.       if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  892.  4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  893.      1' DE LAGRANGE DETECTES')
  894.  
  895.       IF (IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC+NVSTIC
  896.  4821 FORMAT( ' NOMBRE DE VALEURS  DANS LE PROFIL',I12)
  897.       IF (IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR+NVSTIR
  898.  4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  899.       IF (IIMPI.EQ.1) WRITE(IOIMP,4825) Nbopit/1000000
  900.  4825 FORMAT( ' NOMBRE DE GIGA OPERATIONS FMA',I40)
  901.       IF (IIMPI.EQ.1) WRITE(IOIMP,4823) NVaor+NVaori
  902.  4823 FORMAT( ' NOMBRE DE VALEURS initiales',I12)
  903.       INTERR(1)=NVSTOR+NVSTIR
  904.       reaerr(1)=nvstor/inc**(4./3)
  905.       reaerr(2)=2*nbopit/1D6/max(1,(kcour-kcouri))
  906.       reaerr(3)=condmax/condmin
  907.       IF (IPASV.EQ.0.or.reaerr(3).gt.1.D30) CALL ERREUR(-278)
  908.       IPASV=1
  909.       call ooomru(0)
  910.       if(lsgdes) then
  911.         do ipv=1,ino
  912.         lign=mm.ilign(ipv)
  913.         segdes lign
  914.         lign=ll.ilign(ipv)
  915.         segdes lign
  916.         enddo
  917.       endif
  918.       SEGDES,MINCPO
  919.       SEGDES,MIMIK
  920.       SEGDES,MMATRI
  921.       SEGDES,LL,MM
  922.       SEGDES,MDIAG
  923.       MMATRX=MMATRI
  924.       SEGSUP KIVPO,KIVLO
  925.       segsup immt
  926.  9999 continue
  927.       if (ithrd.eq.1) then
  928.         call threadis
  929.        call oooprl(0)
  930.       endif
  931.       RETURN
  932.       END
  933.  
  934.  

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