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

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