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

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