Télécharger la plaquette Cast3M


fr en

Télécharger shole.eso

Retour à la liste

Numérotation des lignes : 

shole

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

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