Télécharger la plaquette Cast3M


fr en

Télécharger shole.eso

Retour à la liste

Numérotation des lignes : 

shole

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

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