Télécharger la plaquette Cast3M


fr en

Télécharger graco2.eso

Retour à la liste

Numérotation des lignes : 

graco2

  1. C GRACO2    SOURCE    PV090527  24/06/12    21:15:05     11940          
  2.       SUBROUTINE GRACO2(MMATRX)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8(A-H,O-Z)
  5. C  TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY
  6. C
  7. C COPIE DE CHOLE POUR FAIRE UN CHOLEVSKI INCOMPLET DE PLUS ON GARDE LA
  8. C MATRICE ASSEMBLEE ( SEGMENT LLIGN)
  9. C
  10.       PARAMETER (LPCRAY=10000)
  11.       INTEGER OOOVAL,OOOLEN
  12.       dimension ittime(4)
  13. C
  14. C  ****  MISE SOUS FORME A=Lt D L  DE LA MATRICE MMATRX
  15. C        SEULES LES VALEURS INITIALEMENT NON NULLES SONT CALCULEES
  16. C
  17.  
  18. -INC PPARAM
  19. -INC CCOPTIO
  20. -INC SMMATRI
  21. -INC CCASSIS
  22. -INC CCHOLE
  23. -INC CCREEL
  24.       SEGMENT KIVPO(IIMAX)
  25.       SEGMENT KIVLO(IIMAX)
  26.  
  27.       external graco5i
  28.       DATA COEAUG/0.75000000000000/
  29.       SAVE IPASV
  30.       DATA IPASV/0/
  31. C     character*8 zen
  32. C     equivalence (zen,izen)
  33.       call timespv(ittime,oothrd)
  34.       kcour=(ittime(1)+ittime(2))/10
  35.       kcourp=kcour
  36.       kcouri=kcour
  37.       kdiff=0
  38.       kcour=0
  39.       nbchan=1
  40.       nbopit=0
  41.       iposm=0
  42. C     zen='CPU'//char(0)
  43. C     le=4
  44.       lolig=0
  45.       nvaor=0
  46.       nbthro=nbthrs
  47.       ithrd=0
  48.       if (nbthro.gt.1) then
  49.        ithrd=1
  50.        call threadii
  51.       endif
  52.       nbthr=nbthro
  53.       do ith=1,nbthro
  54.        nbop(ith)=0
  55.       enddo
  56. C* en fait la parallelisation apporte peu    car la matrice est trop creuse
  57. C*    nbthr=1
  58.       MACTIT=OOOVAL(1,1)
  59.       prec=xspeti
  60.       nvstrm=0
  61.       MMATRI=MMATRX
  62.       SEGACT,MMATRI*MOD
  63.       MILIG1=IILIGN
  64.       IASLIG= IILIGN
  65.       SEGINI,MILIGN=MILIG1
  66.       IILIGN = MILIGN
  67.       INO=ILIGN(/1)
  68.       MDIA1=IDIAG
  69.       IASDIA=IDIAG
  70.       SEGINI,MDIAG=MDIA1
  71. C------------------------
  72.       IDIAG = MDIAG
  73. C-------------------------
  74.       NBLIG=INO
  75. C* comme la matrice est beaucoup plus creuse qu'en direct, on prend un masdim plus petit
  76.       precc=prec
  77.       INC=DIAG(/1)
  78.       INCC=INC
  79.       MIMIK=IIMIK
  80.       MINCPO=IINCPO
  81.       SEGACT,MINCPO,MIMIK
  82.       IPLUMI=IMIK(/2)*2 +4
  83.       IL2=0
  84.       LTRK=MAX(LPCRAY+1,OOOVAL(1,4))
  85.       IIMAX=((((IJMAX+IPLUMI) +LTRK)/LTRK)+1)*LTRK
  86.       SEGINI KIVPO,KIVLO
  87.       INEG=0
  88.       NBLAG=0
  89.       NENSLX=0
  90.       NVSTOC=0
  91.       NVSTOR=0
  92. C
  93. C
  94. C
  95. C  **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  96. C  **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  97. C
  98. C  ****  LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  99. C
  100.   1   CONTINUE
  101.       IVALMA=IJMAX+IPLUMI
  102.       IL1=IL2+1
  103.       IMIN=IL1
  104.  
  105.  
  106.       DO 2 I=IL1,INO
  107.  
  108.       LLIGN=    ILIGN(I)
  109.       SEGACT /ERR=3/LLIGN
  110.  
  111.  
  112.       NA=  IMMMM(/1)
  113. C*    write (6,*) ' chole ligne noeud inconnues ',i,ipno(i),na
  114.       NBPAR=NA+1
  115.       NVALL=  NJMAX
  116.       LMASQ=NVALL/MASDIM+2
  117.  
  118.       SEGINI /ERR=3/ LIGN
  119.  
  120.  
  121.  
  122.  
  123.  
  124.  
  125. C  recuperer la longueur du segment
  126.  
  127.  
  128.  
  129. *  ilu(0) pas de remplissage
  130. *     lglig = (nvall/na)**1.333333333333333333333
  131.       lglig = (nvall/na)**2
  132.  
  133.  
  134.       NVSTOC=NVSTOC + NVALL
  135.       IVALMA=IVALMA + NVALL
  136.  
  137.  
  138.       nvaor = nvaor + xxva(/1)
  139. C
  140. C  ****  DECOMPACTAGE
  141. C
  142.       IPA=1
  143.       valmi=xgrand
  144.       valma=xpetit
  145.       DO 121 JPA=1,NA
  146.       IVPO(JPA)=IPA
  147.       KPA=IPPO(JPA+1)-IPPO(JPA)
  148.       IPP=IPPO(JPA)
  149.       IPPVV(JPA)=IPA-1
  150.       LPA=LDEB(JPA)
  151.       LPA1=LPA-IPA
  152.       DO 122 MPA=1,KPA
  153.       LL=LINC(MPA+IPP)
  154.       IPLA=LL-LPA1
  155.       if(abs(xxva(mpa+ipp)).gt.xpetit/xzprec) then
  156.         VAL(IPLA)=XXVA(MPA+IPP)
  157.         IMASQ(IPLA/MASDIM+1)=1
  158.         XVAL1 = abs(val(ipla))
  159.         valmi=min(valmi,XVAL1)
  160.         valma=max(valma,XVAL1)
  161.       endif
  162. 122   CONTINUE
  163.       IPAHAU=IPA+IMMMM(JPA)-LPA
  164.       IPA=IPA+IMMMM(JPA)-LPA+1
  165.       IMMM(JPA)=IPNO(LPA)
  166.       IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA)
  167.  121  CONTINUE
  168.  
  169.       if (na.gt.0) then
  170.           IPREL=  IMMMM(1)
  171.           IDERL=  IMMMM(NA)
  172.       endif
  173.       IPPVV(NA+1)=IPA-1
  174. C     SEGSUP LLIGN
  175. C     SEGDES LLIGN
  176.       ILIGN(I)=LIGN
  177.       IF(IIMPI.EQ.1525)  THEN
  178.        WRITE( IOIMP,4987)   I
  179.  4987  FORMAT (' NOEUD NUMERO ',I5)
  180.        LL=VAL(/1)
  181.        WRITE(IOIMP, 4926)(VAL(MPA),MPA=1,LL)
  182. 4926  FORMAT(' VAL ' , 10E11.4)
  183.       ENDIF
  184.       nbthro=min(nbthrs,lglig/400+1)
  185.       if (i+1-il1.ge.nbthro) then
  186.        nbthro=min(nbthrs,i+1-il1)
  187.          il2=i
  188.          GOTO 4
  189.       endif
  190.     2 CONTINUE
  191.       IL2=INO
  192.       GO TO 4
  193.     3 IL2=I-1
  194.       if(lign.ne.0) segsup lign
  195.  
  196.     4 CONTINUE
  197.       nbthro=min(nbthrs,nbthro)
  198.       nbthr=nbthro
  199.  
  200.       IF(IL2.GE.IL1) GO TO 40
  201. C
  202. C ****  APPEL AUX ERREURS  MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  203. C
  204.       ITYP=48
  205.       CALL ERREUR(48)
  206.       if (ithrd.eq.1) call threadis
  207.       RETURN
  208.    40 CONTINUE
  209.       IM=INC
  210.  
  211.  
  212.  
  213.  
  214.       DO 352 IH=IL2 ,IL1,-1
  215.       LIGN=   ILIGN(IH  )
  216.       IL=INC
  217.       DO 354 JH=1,    IMMM(/1)
  218.       IM=MIN(IM,    IMMM(JH))
  219.       IL=MIN(IL,    IMMM(JH))
  220.  354  CONTINUE
  221.       IML=IL
  222.       IMM=ipno(IM)
  223.   352 CONTINUE
  224. C  353 CONTINUE
  225.       LIGN=ILIGN(IL1)
  226.       IL11=IPREL
  227. C
  228. C  ****  BOUCLE  *5*  TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  229. C
  230.       ipos=0
  231.       iper=imin
  232.       ider=imin-1
  233.  
  234.       DO 5 I=IMIN,IL2
  235.       if (ierr.ne.0) then
  236.        if (ithrd.eq.1) call threadis
  237.        return
  238.       endif
  239.       LIG1=ILIGN(I)
  240.       IF(I.LT.IL1)  GO TO 7
  241. C
  242. C  ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  243. C  ******* IPREL  IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  244. C  ******* LE TERME DIAGONAL
  245. C
  246.       LIGN=LIG1
  247.       DO 156 KHG=1,IMMM(/1)
  248.       II=IPREL-1+KHG
  249.       IMMM(KHG)=0
  250.       NN=IPPVV(KHG+1)
  251.       NNM1=IPPVV(KHG)
  252.       N=NN-NNM1
  253.       DIAG(II)=VAL(NN)
  254.       IF(N.EQ.1)  GO TO 8
  255.       NMI=N-II
  256.       IDEP=MAX(IL11,2-NMI)
  257.       KIDEP=IDEP+NMI
  258.       KI1=N-1
  259.       KQ=-NMI
  260.       VAL(NN)=VAL(NN)+
  261.      # GRACO3(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI),
  262.      # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG),
  263.      # PREC,ittr(1),nc)
  264.        imasq(nn/masdim+1)=1
  265.        imasq(n/masdim+1)=1
  266.    8  CONTINUE
  267.       if (diag(ii).ne.0.d0) then
  268.        if (ittr(ii).eq.0) then
  269.         if (abs(val(nn)).lt.abs(diag(ii))*coeaug) then
  270. **      write(6,*) 'augmentation raideur ii ',ii,val(nn),diag(ii)
  271.         val(nn)=diag(ii)*coeaug
  272.         endif
  273.        else
  274.         if (abs(val(nn)).lt.abs(diag(ii))*coeaug) then
  275. **      write(6,*) 'augmentation raideur ii ',ii,val(nn),diag(ii)
  276.         val(nn)=diag(ii)*coeaug
  277.         endif
  278.        endif
  279.       endif
  280.  
  281.       IF(ITTR(II).EQ.0.AND.
  282.      & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-10)  GO TO 12
  283.       IF(ITTR(II).NE.0.AND.
  284.      & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-6)  GO TO 12
  285. C
  286. C  **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  287. C  **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  288. C  **** AU PREALABLE SUR CETTE INCONNUE.
  289. C
  290.       VAL(NN)=DIAG(II)*1d-4
  291.       NENS=NENS+1
  292.       IMMM(KHG)=NENS
  293.       if(ittr(ii).NE.0) NENSLX=NENSLX+1
  294.    12 CONTINUE
  295.        IMASQ(NN/MASDIM+1)=1
  296.       DIAG(II)=VAL(NN)
  297.       IF(DIAG(II).NE.0.d0)  GO TO 41
  298.       KQ1=1+NNM1
  299.       KQN=N+NNM1
  300.       DO 16 LFG=KQ1,KQN
  301.       if (imasq(lfg/masdim+1).le.0) goto 16
  302.       IF(VAL(LFG).NE.0.d0)  GO TO 17
  303.    16 CONTINUE
  304.       DIAG(II)=1.D0
  305.       if (ittr(ii).ne.0) diag(ii)=-1d0
  306.       VAL(NN)=DIAG(II)
  307.       GO TO 41
  308.    17 CONTINUE
  309. C
  310. C  ****  ENVOI ERREUR MATRICE SINGUIERE
  311. C
  312.       ITYP=49
  313.       INTERR(1)=I
  314.       CALL ERREUR(49)
  315.       if (ithrd.eq.1) call threadis
  316.       RETURN
  317. C
  318. C  ****  ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  319. C        ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  320. C
  321.    41 IF(DIAG(II).LT.0.D0)  INEG=INEG+1
  322.       IF(ITTR(II).NE.0) NBLAG=NBLAG+1
  323.       diag(ii)=1.d0/diag(ii)
  324.   156 CONTINUE
  325.       NA=IMMM(/1)
  326. C
  327. C  RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE)
  328. C
  329.       izro=0
  330.       CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1),
  331.      #    NVALL,IPPVV(1),izro ,NA ,xpetit/xzprec,imasq(1),iprel,iderl)
  332. C on recree lign car la compacter en place emiette la memoire
  333.        lmasq=0
  334.        segini,lig1
  335. *  deplacement fait ici maintenant, avec unrolling
  336.       do 300 nbp=1,nbpar-1
  337.        kdif =kivpo(nbp)-kivlo(nbp)
  338.        do     iv=kivlo(nbp),kivlo(nbp+1)-4,4
  339.           lig1.val(iv)=val(iv+kdif )
  340.           lig1.val(iv+1)=val(iv+1+kdif )
  341.           lig1.val(iv+2)=val(iv+2+kdif )
  342.           lig1.val(iv+3)=val(iv+3+kdif )
  343.        enddo
  344.        do     iv1=iv,kivlo(nbp+1)-1
  345.           lig1.val(iv1)=val(iv1+kdif )
  346.        enddo
  347.  300  continue
  348.        do it=1,na
  349.         lig1.immm(it)=immm(it)
  350.         lig1.ippvv(it)=ippvv(it)
  351.        enddo
  352.         lig1.ippvv(na+1)=ippvv(na+1)
  353.         lig1.iml=iml
  354.         lig1.iprel=iprel
  355.         lig1.iderl=iderl
  356.         lig1.iml=iml
  357.         lcara(2,i)=lig1.iprel
  358.         lcara(3,i)=lig1.iderl
  359.         lcara(1,i)=lig1.iml
  360.         segsup lign
  361.        lign =lig1
  362.        ilign(i)=lign
  363.  
  364.       NVSTOR=NVSTOR+NVALL
  365.       nvstrm=max(nvstrm,nvall)
  366.       DO 143 LHG=1,NBPAR
  367.        IVPO(2*LHG-1)=KIVPO(LHG)
  368.        IVPO(2*LHG)=KIVLO(LHG)
  369.  143  CONTINUE
  370.       if (i.gt.1) then
  371.        lig1=ilign(i-1)
  372.        segdes lig1
  373.        iderac=min(iderac,i-2)
  374.       endif
  375. C*    WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',VAL(VAL(/1))
  376. C*    WRITE (6,*) ' LIGNE APRES COMPACTION VAL '
  377. C*    WRITE (6,*) (VAL(IK),IK=1,VAL(/1))
  378. C*    WRITE (6,*) 'KIVPO '
  379. C*    WRITE (6,*) (KIVPO(IK),IK=1,NBPAR)
  380. C*    WRITE (6,*) 'KIVLO '
  381. C*    WRITE (6,*) (KIVLO(IK),IK=1,NBPAR)
  382.  
  383. C  ****  ON TRIANGULARISE LES AUTRES LIGNES
  384. C
  385.       il1=il1+1
  386.       if (il1.gt.il2) goto 5
  387.       LIG1=ILIGN(I)
  388.       lign=ilign(il1)
  389.       IL11=IPREL
  390.       goto 7
  391.   71  continue
  392. C      soit parce qu'on a fini, soit parce qu'on manque de memoire
  393. C      il faut executer les lignes activees puis les desactiver
  394. C      lancer les chole3 et attendre qu'ils soient finis
  395.       if (ipos.ne.0) then
  396. **     write (6,*) ' lancement thread ',iper,ider,il1,il2
  397.        if (iper.gt.ider) then
  398.         write (6,*) ' erreur interne chole '
  399.         write (6,*) ' noeud ',i,' sur ',ino
  400.         write (6,*) ' nombre de termes actuels ',nvstor
  401.         call erreur(5)
  402.        endif
  403.        nbthr=min(nbthr,il2-il1+1)
  404.        do ith=1,nbthr-1
  405.         call threadid(ith,graco5i)
  406.        enddo
  407.         call graco5i(nbthr)
  408.        do ith=nbthr-1,1,-1
  409.         call threadif(ith)
  410.        enddo
  411.       endif
  412.  
  413. C  test ctrlC
  414.       if (ierr.ne.0) goto 9999
  415.       iposm=max(iposm,ipos)
  416.       ipos=0
  417.       iderf=ider-1
  418.       if (ider.ne.il1-1) iderf=ider
  419.       do il=iderf,iper,-1
  420.         lign=ilign(il)
  421.         segdes lign
  422. C       write (6,*) ' desactivation ligne ',il
  423.       enddo
  424.        iderac=min(iderac,iper-1)
  425.       iper=ider+1
  426. C     write (6,*) ' iper ider il1 ',iper,ider,il1
  427.       if (iper.ne.il1) goto 7
  428.  
  429.       goto 5
  430.     7 CONTINUE
  431.       SEGACT/err=71/LIG1
  432.       ipos=ipos+1
  433.       ider=i
  434.       if (i.gt.iderac) iderac=i
  435.       if (i.eq.il1-1)   goto 71
  436.     5 CONTINUE
  437. C     write (6,*) ' il1 il2 apres 5 ',il1,il2
  438. C      DO 11 I=IL1,IL2
  439. C        LIGN=   ILIGN(I)
  440.         SEGDES,LIGN
  441.  
  442.  
  443. C   11 CONTINUE
  444.       nbopt=0
  445.       do ith=1,nbthro
  446.        nbopt=nbopt+nbop(ith)
  447.        nbop(ith)=0
  448.       enddo
  449.       nbopin=nbopt
  450.       nbopit=nbopit+nbopin
  451.       call timespv(ittime,oothrd)
  452.       kcourp=kcour
  453.       kcour=(ittime(1)+ittime(2))/10
  454.       kdiff=kcour-kcourp
  455. C*      write (6,*) ' nb operation temps ',nbopin,kdiff
  456.       iderac=min(iderac,il1-1)
  457.       IF(IL2.LT.INO)  GO TO 1
  458. C  ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  459. C   ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  460. C       INEG=INEG-NBLAG
  461. C     on ne compte pas 2 fois les multiplicateurs qui vont etre
  462. C     elimines lors de la resolution car mode d'ensemble
  463.       INEG=INEG-(NBLAG-NENSLX)
  464.       if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  465.  4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  466.      1' DE LAGRANGE DETECTES')
  467.  
  468.       IF(IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC
  469.  4821 FORMAT( ' NOMBRE DE VALEURS  DANS LE PROFIL',I12)
  470.       IF(IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR
  471.  4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  472.       INTERR(1)=NVSTOR
  473.       reaerr(1)=nvstor/inc**(4./3)
  474.       reaerr(2)=2*nbopit/1D3/max(1,(kcour-kcouri))
  475.       IF (IPASV.EQ.0) CALL ERREUR(-278)
  476.       IPASV=1
  477.       SEGDES,MINCPO
  478.       SEGDES,MIMIK
  479.       SEGDES,MMATRI
  480.       SEGDES,MILIGN
  481.       SEGDES,MDIAG
  482.       MMATRX=MMATRI
  483.       SEGSUP KIVPO,KIVLO
  484.  9999 continue
  485.       if (ithrd.eq.1) call threadis
  486.       RETURN
  487.       END
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  
  510.  
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  
  519.  
  520.  
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  
  527.  
  528.  
  529.  
  530.  
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  
  544.  
  545.  
  546.  

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