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

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