Télécharger la plaquette Cast3M


fr en

Télécharger separm.eso

Retour à la liste

Numérotation des lignes : 

separm

  1. C SEPARM    SOURCE    PV090527  23/09/13    21:15:03     11739          
  2.       subroutine separm(mrigid,ri1,ri2,nounil,lagdua,nelim,if)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5.       integer oooval
  6. -INC SMRIGID
  7. -INC SMCOORD
  8. -INC SMELEME
  9.  
  10. -INC PPARAM
  11. -INC CCOPTIO
  12. -INC CCREEL
  13. *  extrait de mrigid dans ri1 la partie des relations pouvant etre
  14. *  traitees comme des dependances (elimination des inconnues)
  15. *  il faut que la premiere inconnue d'une dependance n'apparaisse pas
  16. *  dans les autres dependances
  17. *  juillet 2003 passage aux simples multiplicateurs de Lagrange, mais
  18. *  separm dualise ceux qu'il garde
  19. *-INC CCHAMP
  20.       segment trav1
  21.        character*(lochpo) compp(lcomp1)
  22. C*       character*4 compp(nbincp)
  23.       endsegment
  24.       segment trav2
  25.        integer ielim(nbnoc,nbincp),iautr(nbnoc,nbincp)
  26.        integer icomb(nbnoc,nbincp),ideja(nbnoc,nbincp)
  27.       endsegment
  28.       segment trav3
  29.        integer posinc(nligrp)
  30.       endsegment
  31.       segment trav4
  32.        integer itrv1(nbprl)
  33.        integer itrv2(nbprl)
  34.        real*8 dtrv1(nbprl)
  35.        real*8 dtrv2(nbprl)
  36.       endsegment
  37.       segment icpr(nbpts)
  38.       character*(lochpo) co1
  39.       integer ri1s,ri2s,ri1f
  40. **     write(6,*) ' entree de separm '
  41. **     call ooodmp(0)
  42. **     nbav=oooval(2,1)
  43. *      write (6,*) ' dans separm nouinl ',nounil
  44.       nbdep=0
  45.       nbpiv=0
  47.       nbpvt=0
  48.       nbno=nbpts
  49.       segini icpr
  50.       nbnoc=0
  51.       segact mrigid
  52.       do 760 irig=1,irigel(/2)
  53.         meleme=irigel(1,irig)
  54.         segact meleme
  55. **      write(6,*) '1 itypel dans separm',itypel
  56.          do il=1,num(/2)
  57.          do ip=1,num(/1)
  58.          ipt=num(ip,il)
  59.          if (icpr(ipt).eq.0) then
  60.           nbnoc=nbnoc+1
  61.           icpr(ipt)=nbnoc
  62.          endif
  63.          enddo
  64.          enddo
  65.  760  continue
  66.  
  67. *  remplissage du segment trav1 (liste des noms de composantes)
  68.       lcomp1 = 50
  69.       segini trav1
  70.       nbincp=0
  71.       do 10 irig=1,irigel(/2)
  72. ***    if (irigel(6,irig).ne.0.and.irigel(6,irig).ne.2
  73. ***  &     .and.nounil.eq.0) goto 10
  74.          if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 10
  75.        meleme=irigel(1,irig)
  76.        segact meleme
  77.        if (itypel.ne.22.and.itypel.ne.28) goto 20
  78.        if (irigel(7,irig).ne.0) goto 20
  79.        descr=irigel(3,irig)
  80.        segact  descr
  81.          do 40 nligrp=2,lisinc(/2)
  82. C           print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
  83.            do 50 inc=1,nbincp
  84.              if (compp(inc).eq.lisinc(nligrp)) goto 40
  85.   50       continue
  86.            nbincp=nbincp+1
  87.            if (nbincp.GT.lcomp1) then
  88.              lcomp1 = lcomp1 + 50
  89.              segadj trav1
  90.            endif
  91.            compp(nbincp)=lisinc(nligrp)
  92.   40     continue
  93.        segdes,descr
  94.   20   continue
  95.   10  continue
  96. c*      lcomp1 = nbincp
  97. c*      segadj trav1
  98. *  on construit ri1 et ri2 en notant dans le segment trav2 ce qu'on fait
  99. *     write(6,*) 'separm nbno nbnoc ',nbno,nbnoc
  100.       segini trav2
  101. *  on interdit d'eliminer les inconnues apparaissant dans des conditions unilaterales
  102. *  ainsi qur toutes si on est el elim 0
  103.       do 60 irig=1,irigel(/2)
  104.        meleme=irigel(1,irig)
  105.        if (nelim.ne.0) then
  106.         if(itypel.ne.28) then
  107.         if (irigel(6,irig).eq.0.or.nounil.eq.1) goto 60
  108.         endif
  109.        endif
  110. **      write(6,*) 'itypel dans separm',itypel
  111.        segact meleme
  112. **      if(itypel.eq.28)  write(6,*) 'taille super el ',num(/1)
  113.        if (itypel.ne.22.and.
  114.      >  (itypel.ne.28.or.num(/1).lt.255.or.if.lt.3)) goto 60
  115.        if (irigel(7,irig).ne.0) goto 60
  116.        descr=irigel(3,irig)
  117.        segact  descr
  118.          nld=2
  119.          if(itypel.eq.28) then
  120.          nld=1
  121. **       write(6,*) 'taille super el ',num(/1)
  122.          endif
  123.          do 70 nligrp=nld,lisinc(/2)
  124.            do 80 inc=1,nbincp
  125. C            print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
  126. C     &             ,'  compp(',inc,')=',compp(inc)
  127.              if (compp(inc).eq.lisinc(nligrp)) goto 90
  128.   80       continue
  129. ***      write(ioimp,*) '- ',lisinc(nligrp),' non trouve dans trav1'
  130.            goto 70
  131.   90       continue
  132.            ip=noelep(nligrp)
  133.            do 100 j=1,num(/2)
  134. C            IF(num(ip,j).GT.ielim(/1).OR.inc.GT.ielim(/2).OR.
  135. C     *         num(ip,j).GT.iautr(/1).OR.inc.GT.iautr(/2))THEN
  136. C               print *,'BUG DANS SEPARM : ',
  137. C     &                 num(ip,j),'>',ielim(/1),iautr(/1),
  138. C     &                 inc,'>',ielim(/2),iautr(/2)
  139. C            ENDIF
  140.             ielim(icpr(num(ip,j)),inc)=1
  141.             iautr(icpr(num(ip,j)),inc)=1
  142.  100       continue
  143.   70     continue
  144.        segdes descr
  145.   60   continue
  146.   61   continue
  147. *
  148.        do 700 irig=1,irigel(/2)
  149.         if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 690
  150.         meleme=irigel(1,irig)
  151.         descr=irigel(3,irig)
  152.         xmatri=irigel(4,irig)
  153.         segact meleme
  154.         if (itypel.ne.22) goto 690
  155.         if (irigel(7,irig).ne.0) goto 690
  156.         segact descr
  157. *  remplir posinc pour accelerer les tests sur les composantes
  158.         nligrp=lisinc(/2)
  159.         segini trav3
  160.         do 730 inc=2,nligrp
  161.          do 720 incp=1,nbincp
  162.           if (lisinc(inc).eq.compp(incp)) then
  163.            posinc(inc)=incp
  164.            goto 730
  165.           endif
  166.   720    continue
  167.          call erreur(5)
  168.   730   continue
  169.         do 750 ince=2,lisinc(/2)
  170.          incp=posinc(ince)
  171.          do 740 j=1,num(/2)
  172.           ip=num(noelep(ince),j)
  173.           icomb(icpr(ip),incp)=icomb(icpr(ip),incp)+1
  174.   740    continue
  175.   750   continue
  176.         segsup trav3
  177.   690   continue
  178.   700  continue
  179.       ipass=1
  180.  2000 continue
  181.       segini,ri1=mrigid
  182.       segini,ri2=mrigid
  183. *  quelques preparatifs
  184. *  on  remplit le segment qui sert a trier les relations
  185. *  trier pour attaquer en premier les relations portant sur le moins de noeud
  186.       nbprl=irigel(/2)
  187.       segini trav4
  188.         do 765 irig=1,nbprl        
  189.         descr=irigel(3,irig)
  190.         segact descr
  191.         dtrv1(irig)= lisinc(/2)+0.1
  192.         if(irigel(6,irig).ne.0) dtrv1(irig)=dtrv1(irig)+1000000
  193.         if(irigel(6,irig).eq.2) dtrv1(irig)=dtrv1(irig)+1000000
  194. **      if(irigel(6,irig).eq.0) dtrv1(irig)=dtrv1(irig)+1000000
  195. **      itrv1(irig)=nbprl-irig+1
  196.         itrv1(irig)=      irig   
  197. 765   continue
  198.         call triflo(dtrv1,dtrv2,itrv1,itrv2,nbprl)
  199.         do 200 iri=1,irigel(/2)
  200.         irig=itrv1(iri) 
  201. **      irig=iri   
  202. **      if (ipass.eq.1) irig=itrv1(iri)
  203. *       write(6,*) ' boucle 200 irig iri nbprl ',irig,iri,nbprl
  204.         if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 190
  205.         meleme=irigel(1,irig)
  206.         descr=irigel(3,irig)
  207.         Xmatri=irigel(4,irig)
  208.         segact meleme
  209.         if (itypel.ne.22) goto 190
  210.         if (irigel(7,irig).ne.0) goto 190
  211.         segact descr,Xmatri
  212.         if (abs(re(1,1,1)).gt.1d-30) goto 190
  213. *  remplir posinc pour accelerer les tests sur les composantes
  214.         nligrp=lisinc(/2)
  215.         segini trav3
  216.         do 230 inc=2,nligrp
  217.          do 220 incp=1,nbincp
  218.            if (lisinc(inc).eq.compp(incp)) then
  219.              posinc(inc)=incp
  220.              goto 230
  221.            endif
  222.  220     continue
  223.         call erreur(5)
  224.  230    continue
  225.         nbdepe=0
  226. *  pour arriver à eliminer, il est interessant de faire pivoter les matrices
  227. *  on va le faire en se basant sur le milieu du paquet
  228.         iel=num(/2)/2+1
  229.         incf=0
  230.         do 1000 ince=2,lisinc(/2)
  231.           ipp=num(noelep(ince),iel)
  232.           incpp=posinc(ince)
  233. *  cette inconnue apparait elle dejà dans une autre relation
  234.          if (iautr(icpr(ipp),incpp).eq.1) goto 1000
  235. *  cette matrice touche t elle une autre relation
  236.          do 1240 inc=2,lisinc(/2)
  237.          if (inc.eq.ince) goto 1240
  238.           incp=posinc(inc)
  239.           ip=num(noelep(inc), iel)
  240.           if (ielim(icpr(ip),incp).eq.1) goto 1000
  241. 1240     continue
  242. *  le pivot est il correct
  243.          remax=abs(re(1,1,iel))
  244.          do 1270 ilig=2,re(/2)
  245.            remax=max(remax,abs(re(1,ilig,iel)))
  246. 1270     continue
  247.          remax=remax*0.9
  248.          if (abs(re(1,ince,iel)).lt.remax) goto 1000
  249. *  ok on pivote si il y a lieu
  250. *        write (6,*) ' ipp incpp icomb ',ipp,incpp,icomb(icpr(ipp),incpp)
  251.          if (icomb(icpr(ipp),incpp).eq.1) goto 1010
  252.          if (incf.eq.0) incf=ince
  253.          incf=ince
  254.          goto 1000
  255. 1000    continue
  256.         ince=2
  257.         if (incf.ne.0) ince=incf
  258. 1010    continue
  259.         inm=0
  260.         if (ince.ne.2) then
  261. *  on pivote 2 et ince
  262. *       write (6,*) ' on pivote avec la ligne ',ince
  263.         nbpvt=nbpvt+1
  264.         incp=posinc(ince)
  265.         posinc(ince)=posinc(2)
  266.         posinc(2)=incp
  267.         segini,des1=descr
  268.         segdes descr
  269.         descr=des1
  270.         ri1.irigel(3,irig)=descr
  271.         ri2.irigel(3,irig)=descr
  272. *       write (6,*) noelep(2),noelep(ince),lisinc(2),lisinc(ince)
  273.         co1=lisinc(2)
  274.         lisinc(2)=lisinc(ince)
  275.         lisinc(ince)=co1
  276.         co1=lisdua(2)
  277.         lisdua(2)=lisdua(ince)
  278.         lisdua(ince)=co1
  279.         noe=noelep(2)
  280.         noelep(2)=noelep(ince)
  281.         noelep(ince)=noe
  282.         noe=noeled(2)
  283.         noeled(2)=noeled(ince)
  284.         noeled(ince)=noe
  285.         inm=1
  286.         segini,Xmatr1=Xmatri
  287.         segdes Xmatri
  288.         Xmatri=Xmatr1
  289.         ri1.irigel(4,irig)=Xmatri
  290.         ri2.irigel(4,irig)=Xmatri
  291.         do 1100 im=1,re(/3)
  292.          do 1130 il=1,re(/1)
  293.           ret=re(il,2,im)
  294.           re(il,2,im)=re(il,ince,im)
  295.           re(il,ince,im)=ret
  296.  1130    continue
  297.          do 1160 il=1,re(/2)
  298.           ret=re(2,il,im)
  299.           re(2,il,im)=re(ince,il,im)
  300.           re(ince,il,im)=ret
  301.  1160    continue
  302. *         segdes xmatri
  303.  1100   continue
  304.         endif
  305. *  peut etre une dependance
  306.         segini,ipt1=meleme
  307.         if(inm.eq.0) then
  308.          segini,xmatr1=xmatri
  309.         else
  310.          xmatr1=xmatri  
  311.         endif
  312. *       write(6,*) 'separm xmatr1 ipt1 ipt1.itypel ',
  313. *    >        xmatr1,ipt1,ipt1.itypel
  314.         do 210 j=1,num(/2)
  315.           ip=ipt1.num(noelep(2),j)
  316. *  cette matrice contient t elle des inconnues
  317.           if (noelep(/1).le.1) then
  318.            ipt1.num(1,j)=0
  319.            goto 210
  320.           endif
  321.           incp=posinc(2)
  322.          if (ipass.eq.1.and.icomb(icpr(ip),incp).ne.1) then
  323.             ipt1.num(1,j)=0
  325.             goto 210
  326.          endif
  327.          if (iautr(icpr(ip),incp).eq.1) then
  328.            ipt1.num(1,j)=0
  330.            goto 210
  331.          endif
  332. *  cette matrice touche t elle une autre relation
  333.          do 240 inc=2,lisinc(/2)
  334.           incp=posinc(inc)
  335.           ip=ipt1.num(noelep(inc),j)
  336.           if (ielim(icpr(ip),incp).eq.1) then
  337.             ipt1.num(1,j)=0
  339.             goto 210
  340.           endif
  341.  240      continue
  342. *  le pivot est il correct
  343.           remax=xpetit*1d4
  344.           do 270 ilig=2,re(/2)
  345.             remax=max(remax,abs(re(1,ilig,j)))
  346.  270      continue
  347.           remax=remax*1d-2
  348.           if (abs(re(1,2,j)).le.remax) then
  349.             ipt1.num(1,j)=0
  350.             nbpiv=nbpiv+1
  351.             goto 210
  352.           endif
  353. * Y a t'il deux fois la meme inconnue?
  354.                ideux=0
  355.           do inc=2,lisinc(/2)
  356.          if (ideja(icpr(ipt1.num(noelep(inc),j)),posinc(inc)).eq.1) then
  357.                ideux=inc
  358.             else
  359.               ideja(icpr(ipt1.num(noelep(inc),j)),posinc(inc))=1
  360.             endif
  361.           enddo
  362.           do inc=2,lisinc(/2)
  363.               ideja(icpr(ipt1.num(noelep(inc),j)),posinc(inc))=0
  364.           enddo
  365.           if (ideux.ne.0) then
  366.               moterr(1:4)=lisinc(ideux)
  367.               interr(1)=ipt1.num(noelep(ideux),j)
  368.               call erreur(-361)
  369.               ipt1.num(1,j)=0
  370.                goto 210
  371.           endif
  372. * la matrice est elle augmentee
  373.           if (abs(re(1,1,j)).gt.1d-5) then
  374.             ipt1.num(1,j)=0
  375.             goto 210
  376.           endif
  377. *          segdes xmatri
  378. *  ok on la prend en dependance
  379.           nbdepe=nbdepe+1
  380.           nbdep=nbdep+1
  381. *         write (6,*) ' inconnu eliminee ',ipt1.num(noelep(2),j),
  382. *    >      lisinc(2)
  383.           ielim(icpr(ipt1.num(noelep(2),j)),posinc(2))=1
  384.           do 280 inc=2,lisinc(/2)
  385.             iautr(icpr(ipt1.num(noelep(inc),j)),posinc(inc))=1
  386.  280      continue
  387.  210    continue
  388.         segsup trav3
  389. * il faut maintenant tasser ipt1 et creer ipt2
  390.         nbnn=num(/1)
  391.         nbelem=num(/2)-nbdepe
  392. **      if(nbelem.eq.0) write(6,*) ' separm nbelem ',nbelem
  393. **      write(6,*) ' nombre d elements crees pour ipt2 ',nbelem
  394.         nbsous=0
  395.         nbref=0
  396.         segini ipt2
  397.         ipt2.itypel=itypel
  398.         nelrig=nbelem
  399.         nligrd=re(/1)
  400.         nligrp=re(/2)
  401.         segini xmatr2
  402.         xmatr2.symre=symre
  403.         idec=0
  404.         do 300 j=1,num(/2)
  405.          if (ipt1.num(1,j).eq.0) then
  406.            idec=idec+1
  407.            do 310 i=1,num(/1)
  408.             ipt2.num(i,idec)=num(i,j)
  409.  310       continue
  410.             do io=1,nligrp
  411.                  do iu=1,nligrd
  412.                    xmatr2.re(iu,io,idec)=re(iu,io,j)
  413.                  enddo
  414.             enddo
  415. *           imatr2.imattt(idec)=imattt(j)
  416.          else
  417.            do 320 i=1,num(/1)
  418.             ipt1.num(i,j-idec)=num(i,j)
  419.  320       continue
  420.             do io=1,nligrp
  421.                  do iu=1,nligrd
  422.                    xmatr1.re(iu,io,j-idec)=re(iu,io,j)
  423.                  enddo
  424.             enddo
  425. *           imatr1.imattt(j-idec)=imattt(j)
  426.          endif
  427.  300    continue
  428.  
  429.         nbelem=nbdepe
  430.         segadj ipt1
  431.         nelrig=nbelem
  432.         segadj xmatr1
  433.         ri1.irigel(1,irig)=ipt1
  434.         ri2.irigel(1,irig)=ipt2
  435.         ri1.irigel(4,irig)=xmatr1
  436.         ri2.irigel(4,irig)=xmatr2
  437.         segdes xmatr1,xmatr2
  438.         goto 200
  439.  190    continue
  440. *  pas une dependance
  441.         ri1.irigel(1,irig)=0
  442.  200  continue
  443.       segsup trav4
  444. *  compression de ri1
  445.       idec=0
  446.       do 500 irig=1,irigel(/2)
  447.         meleme=ri1.irigel(1,irig)
  448.         if (meleme.eq.0) then
  449.           idec=idec+1
  450.         else
  451.           segact meleme
  452.           if (num(/2).eq.0) then
  453.             idec=idec+1
  454.             segsup meleme
  455.             xmatri=ri1.irigel(4,irig)
  456.             segsup xmatri
  457.           else
  458.             do 510 ir=1,irigel(/1)
  459.              ri1.irigel(ir,irig-idec)=ri1.irigel(ir,irig)
  460.  510        continue
  461.             ri1.coerig(irig-idec)=ri1.coerig(irig)
  462.           endif
  463.         endif
  464.  500  continue
  465.       nrige=max(irigel(/1),8)
  466.       nrigel=irigel(/2)-idec
  467. *     write (6,*) ' dimension de ri1 ',nrige,nrigel
  468.       segadj ri1
  469.       do irig=1,nrigel
  470.         ri1.irigel(8,irig)=1
  471.       enddo
  472. *  compression de ri2
  473.       idec=0
  474.       do 600 irig=1,irigel(/2)
  475.         meleme=ri2.irigel(1,irig)
  476.         if (meleme.eq.0) then
  477.           idec=idec+1
  478.         else
  479.           segact meleme
  480.           if (num(/2).eq.0) then
  481.             idec=idec+1
  482.             segsup meleme
  483.             xmatri=ri2.irigel(4,irig)
  484.             segsup xmatri
  485.           else
  486.             do 610 ir=1,irigel(/1)
  487.              ri2.irigel(ir,irig-idec)=ri2.irigel(ir,irig)
  488.  610        continue
  489.             ri2.coerig(irig-idec)=ri2.coerig(irig)
  490.           endif
  491.         endif
  492.  600  continue
  493.       nrige=max(irigel(/1),8)
  494.       nrigel=irigel(/2)-idec
  495. *     write (6,*) ' dimension de ri2 ',nrige,nrigel
  496.       segadj ri2
  497.       do irig=1,nrigel
  498.         ri2.irigel(8,irig)=0
  499.       enddo
  500. *  On va voir si on ne peut pas faire pivoter ri2
  501.       if (ipass.eq.2) goto 2001
  502.       ri1s=ri1
  503.       ri2s=ri2
  504.       mrigis=mrigid
  505.       mrigid=ri2
  506.       nbpiv=0
  508.       ipass=ipass+1
  509.       goto 2000
  510.  2001 continue
  511.       mrigid=mrigis
  512.       call fusrig(ri1s,ri1,ri1f)
  513.       ri1=ri1f
  514.       segsup trav1,trav2,icpr
  515.       if (iimpi.ne.0) then
  516.       write (6,*) ' nombre de relations éliminées ',nbdep
  517.       write (6,*) ' nombre de relations gardées à cause du pivot ',
  518.      >            nbpiv
  519.       write (6,*) ' nombre de relations gardées car non indépendantes ',
  521.       write (6,*) ' nombre de paquets pivotés ',
  522.      >            nbpvt
  523.       endif
  524. *  dualisation des multiplicateurs de Lagrange dans ri2
  525. *  si on a des conditions unilatérales, on ne dualise pas, ce sera fait
  526. *  dans le resou de unilater
  527.        iunil=0
  528.        segact ri2
  529.        do ir=1,ri2.irigel(/2)
  530.        if (ri2.irigel(6,ir).ne.0) iunil=1
  531.        enddo
  532.        if ((iunil.eq.0.or.nounil.eq.1).and.lagdua.ge.0)
  533.      $      call dbblx(ri2,lagdua)
  534. *      write(6,*) ' dans separm lagdua ',lagdua
  535. *      write (6,*) ' matrice ri1 '
  536. *      call prrigi(ri1)
  537. *      write (6,*) ' matrice ri2 '
  538. *      call prrigi(ri2)
  539. **     write(6,*) ' sortie de separm '
  540. **     call ooodmp(0)
  541. **    nbap=oooval(2,1)
  542. **    write(6,*) 'nb segmts dans separm avant apres ',nbav,nbap
  543.  
  544.       end
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  

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