Télécharger la plaquette Cast3M


fr en

Télécharger chole.eso

Retour à la liste

Numérotation des lignes : 

chole

  1. C CHOLE     SOURCE    PV090527  24/04/13    21:15:03     11827          
  2. *
  3.       subroutine chole(mmatrx,prec,istab,nbnnma,nligra,xmatri)
  4. *
  5. * octobre 2022 version unifie chole chomod: normal et condensation superelement
  6. *
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8(A-H,O-Z)
  9. C  TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY
  10.       PARAMETER (LPCRAY=10000)
  11.       INTEGER OOOVAL,OOOLEN
  12.       dimension ittime(4)
  13.  
  14.       POINTEUR LILIGN.MILIGN
  15. C
  16. C  ****  MISE SOUS FORME A=Lt D L  DE LA MATRICE MMATRX
  17. C
  18.  
  19. -INC PPARAM
  20. -INC CCOPTIO
  21. -INC CCREEL
  22. -INC SMMATRI
  23. -INC SMRIGID
  24.  
  25. -INC CCASSIS
  26. -INC CCHOLE
  27.       SEGMENT KIVPO(IIMAX)
  28.       SEGMENT KIVLO(IIMAX)
  29.       segment immt(nblig)
  30.       segment ireser(nvstrm)
  31.       external chole3i
  32.       SAVE IPASV
  33.       DATA IPASV/0/
  34. *  ngmpet dit si on tient en memoire (false) ou si on deborde (true)
  35.       logical ngmpet
  36. C     character*8 zen
  37. C     equivalence (zen,izen)
  38.       logical lsgdes,pasfait,ngdyn
  39.       ireser=0
  40.       nbres=0
  41.       matric=xmatri
  42.       maitre=0
  43.       nelrig=1
  44.       nligrp=nligra
  45.       nligrd=nligra
  46.       if (xmatri.ne.0) then
  47.        segini xmatri
  48.       endif
  49.       nbnnmc=nbnnma
  50.       matric=xmatri
  51.  
  52.  
  53.  
  54. 10000 continue
  55.       pasfait=.true.
  56.       lsgdes=.false.
  57. *  faire attention a respecter l'ordre des segdes par la suite
  58.       call ooomru(1)
  59.       if (istab.ne.0.and.istab.ne.1) call erreur(5)
  60.       condmax=0.d0
  61.       condmin=xgrand*xzprec
  62.       ngmpet=.false.
  63.       ngdyn=.true.
  64.       call timespv(ittime,oothrd)
  65.       kcour=(ittime(1)+ittime(2))/10
  66.       kcouri=kcour
  67.       kcour=0
  68.       nbopit=0
  69.  
  70.       iposm=0
  71. C     zen='CPU'//char(0)
  72. C     le=4
  73.       nvaor=0
  74.       nbthro=nbthrs
  75.       ithrd=0
  76.       if (nbthro.gt.1) then
  77.        ithrd=1
  78.        call threadii
  79.        call oooprl(1)
  80.       endif
  81.       nbthr=nbthro
  82.       do ith=1,nbthr
  83.        nbop(ith)=0
  84.       enddo
  85.       stmult=1d-5
  86.  
  87. C  nouvelle methode de gestion de l'espace memoire necessitee par la parallelisation
  88. C  memoire vive totale
  89.       MACTIT=OOOVAL(1,1)
  90. **    write(6,*) ' mactit igrand ',mactit,igrand
  91. C  un bloc de memoire fera au plus macti/2
  92.       call intpdo(inpdo)
  93.       nvstrm=mactit/10
  94.       MMATRI=MMATRX
  95.       SEGACT,MMATRI*MOD
  96.       PRCHLV=PREC
  97.       MILIGN=IILIGN
  98.  
  99.       SEGACT,MILIGN*MOD
  100.       INO=ILIGN(/1)
  101.       MDIAG=IDIAG
  102.       SEGACT,MDIAG*MOD
  103.       NBLIG=INO
  104.       segini immt
  105.       precc=prec
  106.       INC=DIAG(/1)
  107.       if (matric.eq.0) nbnnmc=inc+1
  108.       nvstrm=max(inc*inpdo,nvstrm)
  109. **      write(6,*) ' nvstrm ',nvstrm
  110.       INCC=INC
  111.       MIMIK=IIMIK
  112.       MINCPO=IINCPO
  113.       SEGACT,MINCPO,MIMIK
  114.       IPLUMI=IMIK(/2)*2 +4
  115.       IL2=0
  116.       IIMAX=IJMAX+IPLUMI
  117.       SEGINI KIVPO,KIVLO
  118.       INEG=0
  119.       NBLAG=0
  120.       NENSLX=0
  121.       NVSTOC=0
  122.       NVSTOR=0
  123.       diagmax=XPETIT/XZPREC
  124.       diagmin=xgrand
  125.       do i=1,diag(/1)
  126.        if (ittr(i).eq.0) diagmax=max(diagmax,abs(diag(i)))
  127.        if (ittr(i).eq.0.and.abs(diag(i)).gt.xpetit/xzprec)
  128.      >   diagmin=min(diagmin,abs(diag(i)))
  129.       enddo
  130.       if (diagmax.le.xpetit/xzprec) then
  131.          do i=1,diag(/1)
  132.            diagmax=max(diagmax,abs(diag(i)))
  133.            if (abs(diag(i)).gt.xpetit/xzprec)
  134.      >       diagmin=min(diagmin,abs(diag(i)))
  135.          enddo
  136.       endif
  137.       diagmin=min(diagmin,diagmax)
  138. ***    write (6,*) ' chole diagmin diagmax ',diagmin,diagmax,diag(/1)
  139.  
  140.  
  141. C
  142. C
  143. C
  144. C  **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  145. C  **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  146. C
  147. C  ****  LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  148. C
  149. C SP  indicateurs pour impression message "stabilisation RESO..."
  150.       isr=0
  151.       isrl=0
  152.  
  153.   1   CONTINUE
  154.       IVALMA=IJMAX+IPLUMI
  155.       IL1=IL2+1
  156.       IMIN=IL1
  157. *  reserver de la place ou mettre les lignes superieures dans le cas debordement
  158.       if (ngmpet) then
  159.         if (ireser.eq.0) segini ireser
  160.       endif
  161.       DO 2 I=IL1,INO
  162.       ngdyn=.true.
  163. **    imasq=0
  164.       LIGN=0  
  165.  
  166.       LLIGN=    ILIGN(I)
  167.       SEGACT /ERR=32/LLIGN
  168.       goto 31
  169.   32  continue
  170. **       write(6,*) ' segact llign erreur',i,il1,lsgdes
  171.       if (.not.lsgdes) then
  172. **         write(6,*) ' lsgdes 1 '
  173.         lsgdes=.true.
  174. *****   ngmpet=.true.
  175. **      write(6,*) 'desactivation-1 ',1,il1-1
  176.         do it=il1-1,1,-1
  177.          lign=ilign(it)
  178.          segdes lign
  179.         enddo
  180.       else
  181.         goto 3
  182.       endif
  183.       SEGACT /ERR=3/LLIGN
  184.   31  continue
  185.       NA=  IMMMM(/1)
  186. *     write (6,*) ' chole ligne noeud inconnues ',i,ipno(i),na
  187.       NBPAR=NA+1
  188.       NVALL=  NJMAX
  189.       LMASQ=NVALL/MASDIM+2
  190.       SEGINI /ERR=33/ LIGN
  191.       goto 34
  192.   33  continue
  193. **     write(6,*) ' segini lign erreur',il1
  194.       if (.not.lsgdes) then
  195. **        write(6,*) ' lsgdes 2 '
  196.         lsgdes=.true.
  197. *****   ngmpet=.true.
  198. **      write(6,*) 'desactivation-2 ',1,il1-1
  199.         do it=il1-1,1,-1
  200.          lign=ilign(it)
  201.          segdes lign
  202.         enddo
  203.       else
  204.         goto 3
  205.       endif
  206.       SEGINI /ERR=3/ LIGN
  207. **     write(6,*) 'deuxieme essai lign ok'
  208.   34  continue
  209. C  recuperer la longueur du segment
  210.       lglig=na*(nvall/na)**1.333333333333333333
  211.       NVSTOC=NVSTOC + NVALL
  212.       IVALMA=IVALMA + NVALL
  213.  
  214.  
  215.       nvaor = nvaor + xxva(/1)
  216. C
  217. C  ****  DECOMPACTAGE
  218. C
  219.       IPA=1
  220.       DO 121 JPA=1,NA
  221.             IVPO(JPA)=IPA
  222.         KPA =   IPPO(JPA+1)-  IPPO(JPA)
  223.         IPP =   IPPO(JPA)
  224.             IPPVV(JPA)=IPA-1
  225.         LPA  =   LDEB(JPA)
  226.         LPA1 = LPA-IPA
  227.  
  228.         DO 122 MPA=1,KPA
  229.           LL   =   LINC(MPA+IPP)
  230.           IPLA = LL -LPA1
  231.            xxv=xxva(mpa+ipp)
  232.            if (abs(xxv).gt.xpetit) then
  233.               VAL(IPLA)=  XXV
  234.               IMASQ(IPLA/MASDIM+1)=1
  235.               if (ipla-ipa+1.ge.1) IMASQ((IPLA-ipa+1)/MASDIM+1)=1
  236.            endif
  237. 122     CONTINUE
  238.  
  239.         IPA=IPA+  IMMMM(JPA)-LPA + 1
  240. Cpv         IMMM(JPA)=   IPNO(LPA)
  241.         IMMM(JPA)=LPA
  242.         IF(IMIN  .GT.   IPNO(LPA )) IMIN  =    IPNO(LPA)
  243.  
  244.  
  245.  
  246. 121   CONTINUE
  247.  
  248. * indexation de imasq
  249.         ipln=lmasq/na
  250.         iplp=lmasq/na
  251.         do 123 ipl=lmasq/na,1,-1
  252.           if (imasq(ipl).gt.0) then
  253.             imasq(ipl)=iplp*masdim
  254.             ipln=ipl-1
  255.           else
  256.             imasq(ipl)=-ipln*masdim
  257.             iplp=ipl-1
  258.           endif
  259.  123    continue
  260. **      write (6,*) ' imasq ',lmasq/na
  261. **      write (6,*) (imasq(ipl),ipl=1,lmasq/na)
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270. C*** **** ****
  271.  
  272.  
  273.  
  274.  
  275.  
  276.  
  277.  
  278.  
  279.  
  280.  
  281.  
  282.       if (na.gt.0) then
  283.           IPREL=  IMMMM(1)
  284.           IDERL=  IMMMM(NA)
  285.           lcara(2,i)=iprel
  286.           lcara(3,i)=iderl
  287.       endif
  288.           IPPVV(NA+1)=IPA-1
  289.  
  290.  
  291.  
  292.  
  293.  
  294.       SEGSUP LLIGN
  295.          ILIGN(I)=LIGN
  296.  
  297.  
  298. C*    write (6,*) 'longueur ligne ',nvall
  299. C nb de ligne multiple du nb de threads
  300. C blocage ligne lecture-ecriture pour minimiser le cache
  301. C on note si on est au minimum de lignes
  302.       nbthro=min(nbthrs,lglig/1200+1)
  303.       if (i+1-il1.ge.nbthro.and.(.not.ngmpet)) then
  304.        nbthro=min(nbthrs,i+1-il1)
  305.        ngdyn=.true.
  306.         if(i+1-il1.eq.nbthrs) ngdyn=.false.
  307.          il2=i
  308. **       write(6,*) ' nbthro il1 il2 ',nbthro,il1,il2
  309.  
  310.  
  311.          GOTO 4
  312.       endif
  313.  
  314.  
  315.     2 CONTINUE
  316.       IL2=INO
  317.       GO TO 4
  318.     3 IL2=I-1
  319. **      write(6,*) 'desactivation-4 ',llign
  320.       segdes llign
  321.       if (               lign.ne.0) segsup lign
  322.     4 CONTINUE
  323.       nbthro=min(nbthrs,nbthro)
  324.       nbthr=nbthro
  325.       if(ireser.ne.0) segsup ireser
  326. C
  327.       IF(IL2.GE.IL1) GO TO 40
  328. C
  329. C ****  APPEL AUX ERREURS  MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  330. C
  331. C      ITYP=48
  332.       CALL ERREUR(48)
  333.         call ooodmp(0)
  334.       if (ithrd.eq.1) then
  335.        call threadis
  336.        call oooprl(0)
  337.       endif
  338.       call ooomru(0)
  339.       RETURN
  340.    40 CONTINUE
  341.       IM=INC
  342.  
  343.  
  344.  
  345.       DO 352 IH=IL2 ,IL1,-1
  346.       LIGN=   ILIGN(IH  )
  347.       IL=INC
  348.       DO 354 JH=1,    IMMM(/1)
  349.       IM=MIN(IM,    IMMM(JH))
  350.       IL=MIN(IL,    IMMM(JH))
  351.  354  CONTINUE
  352.       IML=IL
  353.       lcara(1,ih)=iml
  354.       immt(ih)=ipno(im)
  355.   352 CONTINUE
  356. C  353 CONTINUE
  357.       LIGN=ILIGN(IL1)
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.       IL11=    IPREL
  375.  
  376.  
  377. C
  378. C  ****  BOUCLE  *5*  TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  379. C
  380. C    >      il1,il2
  381.       lig1=   ilign(imin)
  382.       ipos=0
  383.       iper=imin
  384.       ider=imin-1
  385.       iderac=imin-1
  386.  
  387.       IL1i=il1
  388.       DO 5 I=IMIN ,IL2
  389.  
  390.       LIG1=   ILIGN(I)
  391.  
  392.       IF(I.LT.IL1)  GO TO 7
  393. C____________
  394. C
  395. C  ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  396. C  ******* IPREL  IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  397. C  ******* LE TERME DIAGONAL
  398. C
  399.       LIGN=LIG1
  400.       DO 156 KHG=1,IMMM(/1)
  401.  
  402.       II=IPREL-1+KHG
  403.       IMMM(KHG)=0
  404.       NN=IPPVV(KHG+1)
  405.  
  406.       NNM1=IPPVV(KHG)
  407.  
  408.       N=NN-NNM1
  409.       DIAG(II)=VAL(NN)
  410.        if (n.eq.1) then
  411.         if(ii-nbnnmc.gt.0) then
  412.          re(ii-nbnnmc,ii-nbnnmc,1)=val(nn)
  413.          goto 41
  414.         else
  415.          goto 8
  416.         endif
  417.        endif
  418.  
  419.  
  420.       NMI=N-II
  421. *  attention subtile difference noeud maitre et non maitre
  422.       if (ii-nbnnmc.le.0) then
  423.        IDEP=MAX(IL11,2-NMI)
  424.       else
  425.        IDEP=MAX(IL11,1-NMI)
  426.       endif
  427.       KIDEP=IDEP+NMI
  428.       KI1=N-1
  429.       KQ=-NMI
  430.       VAL(NN)=VAL(NN)+
  431.      > CHOLE1(ILIGN,IPREL,imasq(1),VAL(1+IPPVV(KHG)),DIAG(1-NMI),
  432.      > IPNO(1-NMI),IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,1+IPPVV(KHG),
  433.      > PREC,nbop(1))
  434.        imasq(nn/masdim+1)=1
  435.        imasq(n/masdim+1)=1
  436.        if(ii-nbnnmc.gt.0) then
  437.          re(ii-nbnnmc,ii-nbnnmc,1)=val(nn)
  438.          goto 41
  439.        endif
  440.    8  CONTINUE
  441.       diagref=max(abs(diag(ii)),diagmin)
  442.       diagcmp=diagref*5d-12
  443.       IF(   ITTR(II).EQ.0.AND.
  444.      & ABS(     VAL(NN)).GT.diagcmp)  GO TO 12
  445.       IF(   ITTR(II).NE.0.AND.
  446.      & ABS(     VAL(NN)).GT.diagcmp)  GO TO 12
  447. **      write(6,*) ' ittr val diagcmp ',ittr(ii),val(nn),diagcmp
  448. C  il faut mettre une valeur plus grande sur les LX car on a un probleme de conditionnement
  449. C  sur le calcul des reactions en cas de 2 relations presque identique
  450. C
  451. C  **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  452. C  **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  453. C  **** AU PREALABLE SUR CETTE INCONNUE.
  454. C
  455. *     write (6,*) ' chole mode d ensemble ittr ligne ',
  456. *    >      ittr(ii),ii,diag(ii),val(nn),diagref
  457. C  on garde le signe car il fau un moins sur les ML
  458.        vmaxi=diagref
  459.        do ipv=1+ippvv(khg),nn
  460.         vmaxi=max(vmaxi,abs(val(ipv)))
  461.        enddo
  462.       if(   ittr(ii).NE.0) then
  463.        VAL(NN)=val(nn)-1.30D0*diagref
  464.        NENSLX=NENSLX+1
  465.       else
  466.        val(nn)=vmaxi
  467.       endif
  468.       NENS=NENS+1
  469.            IMMM(KHG)=NENS
  470.    12 CONTINUE
  471.  
  472.  
  473. *  stabilisation
  474.       IF (ISTAB.NE.0) THEN
  475. *  elimination des Nan
  476.        if (.not.(abs(val(nn)).lt.xgrand*xzprec).and.pasfait) then
  477.          pasfait=.false.
  478.          val(nn)=xgrand*xzprec
  479.          write (6,*) 'Nan dans chole ligne',ii
  480.        endif
  481. *
  482.        diagcmp=abs(diagmin)*1d-5+xpetit/xzprec
  483.        if (val(nn).lt.-diagmax*1d-3.and.ittr(ii).eq.0) then
  484.         val(nn)=abs(val(nn))
  485. ***     val(nn)=diagmax*1d-3
  486.         if (immm(khg).eq.0) NENS=NENS+1
  487.         IMMM(KHG)=NENS
  488.        elseif (val(nn).le.diagcmp.and.ittr(ii).eq.0 ) then
  489.         if (isr.eq.0.or.iimpi.gt.0)
  490.      &     write (6,*) ' stabilisation RESO ',ii,val(nn),diag(ii)
  491.         isr=isr+1
  492.         val(nn)=max(diagcmp,-val(nn))
  493.         if (immm(khg).eq.0) NENS=NENS+1
  494.         IMMM(KHG)=NENS
  495.        else
  496.         if (ittr(ii).eq.0) diagmin=min(diagmin,abs(val(nn)))
  497.        endif
  498.        if (val(nn).ge.abs(diag(ii))*stmult.and.ittr(ii).ne.0) then
  499.         if (isrl.eq.0.or.iimpi.gt.0)
  500.      &     write (6,*) ' stabilisation RESO lagrange ',ii,val(nn)
  501.         isrl=isrl+1
  502.         val(nn)=-abs(diag(ii))*stmult
  503.         if (immm(khg).eq.0) then
  504.            NENS=NENS+1
  505.            NENSLX=NENSLX+1
  506.         endif
  507.         IMMM(KHG)=NENS
  508.        endif
  509.       ENDIF
  510.  
  511.  
  512.  
  513.  
  514.  
  515.       DIAG(II)=     VAL(NN)
  516.       IF(abs(DIAG(II)).gt.xpetit)  GO TO 41
  517.       DIAG(II)=diagmax
  518.       if (ittr(ii).ne.0) diag(ii)=-diagmax
  519.       VAL(NN)=DIAG(II)
  520.       GO TO 41
  521. C
  522. C  ****  ENVOI ERREUR MATRICE SINGUIERE
  523. C
  524. C      ITYP=49
  525.       INTERR(1)=I
  526.       CALL ERREUR(49)
  527.       if (ithrd.eq.1) then
  528.         call threadis
  529.        call oooprl(0)
  530.       endif
  531.       call ooomru(0)
  532.       RETURN
  533. C
  534. C  ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  535. C        ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  536. C
  537.    41 IF(DIAG(II).LT.0.D0)  INEG=INEG+1
  538.       IF(   ITTR(II).NE.0) NBLAG=NBLAG+1
  539.  
  540.       if (ii.le.nbnnmc) condmin=min(condmin,abs(diag(ii)))
  541.       condmax=max(condmax,abs(diag(ii)))
  542.       if (ii.le.nbnnmc) diag(ii)=1.d0/diag(ii)
  543.   156 CONTINUE
  544.       NA=IMMM(/1)
  545. C
  546. C  RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE)
  547. C
  548. *       write(6,*) 'chole ligne lpl',i,na,ippvv(2)-ippvv(1)
  549.       if (na.gt.0)
  550.      >CALL COMPAC(     VAL(1),NBPAR,KIVPO(1),KIVLO(1),
  551.      #    NVALL,     IPPVV(1),IZROSF,NA,PREC,imasq(1),iprel,iderl)
  552.       lmasq=0
  553. C on recree lign car la compacter en place emiette la memoire
  554.        lig1=lign
  555.        segini /err=1431/ lig1
  556.        goto 1432
  557.  1431  continue
  558.        write(6,*) ' segini echoue'
  559.  1432  continue
  560.       nbres=max(nbres,nvall)
  561. *  deplacement fait ici maintenant, avec unrolling
  562.       do 300 nbp=1,nbpar-1
  563.        kdif =kivpo(nbp)-kivlo(nbp)
  564.        do     iv=kivlo(nbp),kivlo(nbp+1)-4,4
  565.           lig1.val(iv)=val(iv+kdif )
  566.           lig1.val(iv+1)=val(iv+1+kdif )
  567.           lig1.val(iv+2)=val(iv+2+kdif )
  568.           lig1.val(iv+3)=val(iv+3+kdif )
  569.        enddo
  570.        do     iv1=iv,kivlo(nbp+1)-1
  571.           lig1.val(iv1)=val(iv1+kdif )
  572.        enddo
  573.  300  continue
  574. *      do it=1,nvall
  575. *       lig1.val(it)=val(it)
  576. *      enddo
  577.        do it=1,na
  578.         lig1.immm(it)=immm(it)
  579.         lig1.ippvv(it)=ippvv(it)
  580.        enddo
  581.         lig1.ippvv(na+1)=ippvv(na+1)
  582.         lig1.iml=iml
  583.         lig1.iprel=iprel
  584.         lig1.iderl=iderl
  585.         lig1.iml=iml
  586.         lcara(2,i)=lig1.iprel
  587.         lcara(3,i)=lig1.iderl
  588.         lcara(1,i)=lig1.iml
  589.        if (lign.ne.lig1) then
  590.         segsup lign
  591.         lign=lig1
  592.        else
  593.         segadj lign
  594.        endif
  595.           ilign(i)=lign
  596.       NVSTOR=NVSTOR+NVALL
  597.       nvstrm=max(nvstrm,nvall*inpdo)
  598.       DO 143 LHG=1,NBPAR
  599.             IVPO(2*LHG-1)=KIVPO(LHG)
  600.             IVPO(2*LHG)  =KIVLO(LHG)
  601.  143  CONTINUE
  602. C   Si la ligne est petite, il n'y a rien a gagner a paralleliser
  603.  
  604.       if (i.gt.1) then
  605.        lig1=ilign(i-1)
  606. **      if(lsgdes) write(6,*) 'desactivation-5 ',i
  607.        if (lsgdes)  segdes lig1
  608.        iderac=min(iderac,i-2)
  609.       endif
  610.  
  611. C
  612. C  ****  ON TRIANGULARISE LES AUTRES LIGNES
  613. C
  614.       il1=il1+1
  615.       if (il1.gt.il2) goto 5
  616.       LIG1=ILIGN(I)
  617.       lign=ilign(il1)
  618.       IL11=IPREL
  619.       goto 7
  620. C  72  continue
  621.   71  continue
  622. *  passage en superlent
  623.       if (ider.lt.il1-1.and..not.ngmpet) then
  624. **      write(6,*) 'ngmpet true ',iper,ider,il1,il2
  625.         ngmpet=.true.
  626. **      do ipv=1,iper-1
  627. **        lig2=ilign(ipv)
  628. **       call oooeta(lig2,ieta,imod)
  629. **         if (ieta.eq.1) write(6,*) 'ligne active ',ipv,iper
  630. **      enddo
  631.       endif
  632.       if (iper.gt.ider) then
  633.         call erreur(48)
  634.         write(6,*) ' 1 chole iper ider i il1 il2 ',
  635.      >         iper,ider,i,il1,il2
  636.         do ipv=1,ino
  637.          lig2=ilign(ipv)
  638.          call oooeta(lig2,ieta,imod)
  639.          if (ipv.lt.il1.or.ipv.gt.il2) then
  640.            if (ieta.eq.1) write(6,*) 'ligne active ',ipv
  641.          endif
  642.         enddo
  643.         call ooodmp(0)
  644.       if (ithrd.eq.1) then
  645.        call threadis
  646.        call oooprl(0)
  647.       endif
  648.       call ooomru(0)
  649.         return
  650.       endif
  651.  
  652. C      soit parce qu'on a fini, soit parce qu'on manque de memoire
  653. C      il faut executer les lignes activees puis les desactiver
  654. C      lancer les chole3 et attendre qu'ils soient finis
  655.       if (ipos.ne.0) then
  656. **      write (6,*) ' lancement thread ',iper,ider,il1,il2
  657.         if (iper.gt.ider) then
  658.           write (6,*) ' erreur interne chole '
  659.           call erreur(5)
  660.         endif
  661.         nbthr=min(nbthr,il2-il1+1)
  662. **      write(6,*) 'nbthr ',nbthr,il2-il1+1
  663.          LILIGN=MILIGN
  664. *  blocage pour rester dans le cache secondaire
  665.         ipers=iper
  666.         iders=ider
  667.         ipas=1500
  668.         if(nbthr.eq.1) ipas=igrand
  669. *       do 400 iper=ipers,iders,ipas
  670.  401    continue
  671.         ider=min(iders,iper+ipas-1)
  672.           do ith=1,nbthr-1
  673.             call threadid(ith,chole3i)
  674.           enddo
  675.           call chole3i(nbthr)
  676.           do ith=nbthr-1,1,-1
  677.             call threadif(ith)
  678.           enddo
  679.         iper=iper+ipas
  680.         ipas=ipas/2
  681.         ipas=max(ipas,750)
  682.         if (iper.le.iders) goto 401
  683.  400    continue
  684.         iper=ipers
  685.         ider=iders
  686.       endif
  687.  
  688. C  test ctrlC
  689.       if (ierr.ne.0) goto 9999
  690.       iposm=max(iposm,ipos)
  691.       ipos=0
  692.       iderf=ider-1
  693.       if (ider.ne.il1-1) iderf=ider
  694.       if (lsgdes) then
  695. **      write(6,*) 'desactivation 7 iper iderf',iper,iderf
  696.       do il=iderf,iper,-1
  697.         lign=ilign(il)
  698.         segdes lign
  699. C       write (6,*) ' desactivation ligne ',il
  700.       enddo
  701.       endif
  702.        iderac=min(iderac,iper-1)
  703.       iper=ider+1
  704. C     write (6,*) ' iper ider il1 ',iper,ider,il1
  705.       if (iper.ne.il1) goto 7
  706.  
  707.       goto 5
  708.     7 CONTINUE
  709.  
  710.       if (lsgdes) SEGACT/err=71/LIG1
  711.       ipos=ipos+1
  712.       ider=i
  713.       if (i.gt.iderac) iderac=i
  714.       if (i.eq.il1-1)   goto 71
  715.     5 CONTINUE
  716. **    write (6,*) ' il1 il2 apres 5 ',il1,il2
  717.       if (lsgdes) then
  718. **      write(6,*) 'desactivation 8 il1 il2',il1,il2
  719.       DO I=min(il1,il2),max(il1,il2)
  720.         LIGN=   ILIGN(I)
  721.         if(lign.ne.0) SEGDES,LIGN
  722.       enddo
  723.       endif
  724.       nbopt=0
  725.       do ith=1,nbthro
  726.        nbopt=nbopt+nbop(ith)
  727.        nbop(ith)=0
  728.       enddo
  729.       nbopin=nbopt
  730.       nbopit=nbopit+nbopin
  731.       call timespv(ittime,oothrd)
  732.       kcour=(ittime(1)+ittime(2))/10
  733.  
  734.  
  735.  
  736.  
  737.  
  738.       iderac=min(iderac,il1-1)
  739.       if (ireser.ne.0) segsup ireser
  740.       IF(IL2.LT.INO)  GO TO 1
  741. C  ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  742. C   ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  743. C       INEG=INEG-NBLAG
  744. C     on ne compte pas 2 fois les multiplicateurs qui vont etre
  745. C     elimines lors de la resolution car mode d'ensemble
  746.       INEG=INEG-(NBLAG-NENSLX)
  747.       if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  748.  4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  749.      1' DE LAGRANGE DETECTES')
  750.  
  751.       IF (IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC
  752.  4821 FORMAT( ' NOMBRE DE VALEURS  DANS LE PROFIL',I12)
  753.       IF (IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR
  754.  4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  755.       IF (IIMPI.EQ.1) WRITE(IOIMP,4825) Nbopit/1000000
  756.  4825 FORMAT( ' NOMBRE DE GIGA OPERATIONS FMA',I40)
  757.       IF (IIMPI.EQ.1) WRITE(IOIMP,4823) NVaor
  758.  4823 FORMAT( ' NOMBRE DE VALEURS initiales',I9)
  759. C     IF (IIMPI.EQ.1) WRITE(IOIMP,4824) nbopit/1d6/(kcour-kcouri)
  760. C 4824 FORMAT( ' Performance en Gigaflops ',F8.1)
  761.       INTERR(1)=NVSTOR
  762.       reaerr(1)=nvstor/inc**(4.D0/3.D0)
  763.       reaerr(2)=2.D0*nbopit/1.D6/max(1,(kcour-kcouri))
  764.       reaerr(3)=condmax/condmin
  765.       IF (IPASV.EQ.0.or.reaerr(3).gt.1.D30) CALL ERREUR(-278)
  766.       IPASV=1
  767.       call ooomru(0)
  768.       if(lsgdes) then
  769.         do ipv=1,ino
  770.           lign=ilign(ipv)
  771.           segdes lign
  772.         enddo
  773.       endif
  774.       SEGDES,MINCPO
  775.       SEGDES,MIMIK
  776.       SEGDES,MMATRI
  777.       SEGDES,MILIGN
  778.       MMATRX=MMATRI
  779.       SEGSUP KIVPO,KIVLO
  780.       segsup immt
  781. C     write (6,*) ' chole ipos max ',iposm
  782.  9999 continue
  783.       if (ithrd.eq.1) then
  784.         call threadis
  785.        call oooprl(0)
  786.       endif
  787.       if(iimpi.ne.0)  write (6,*) 
  788.      >      ' chole condmin condmax ',condmin,condmax,diag(/1)
  789.       SEGDES,MDIAG
  790. **    write(6,*) 'kseq kpar',xkseq,xkpar
  791.  
  792.  
  793.       RETURN
  794.       END
  795.  
  796.  
  797.  
  798.  

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