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

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