Télécharger la plaquette Cast3M


fr en

Télécharger scnd2.eso

Retour à la liste

Numérotation des lignes : 

scnd2

  1. C SCND2     SOURCE    PV090527  26/04/30    21:16:24     12529          
  2. *
  3.       subroutine scnd2(mrigid,ri1,ri4)
  4. *
  5. *    entree
  6. *         mrigid rigidité  sans sous matrices de dependances
  7. *         ri1 rigidite  de dependance
  8. *
  9. *    sortie
  10. *         ri4     rigidite condensee les noeuds dependants sont
  11. *                 elimines
  12. *
  13. *    on suppose que la premiere inconnue de la raideur de dependance
  14. *    (celle qu'on elimine) n'apparait que dans cette dependance
  15. *    Cela va permettre de constituer les tableau :
  16. *    noeu-inconnue (primale ou duale) ==> numéro raideur de dépendance
  17. *    au moment de la création des matrices, on va s'efforcer de factoriser
  18. *    les descripteurs
  19. *
  20.       IMPLICIT INTEGER(I-N)
  21.       IMPLICIT REAL*8 (A-H,O-Z)
  22. *
  23. -INC SMRIGID
  24. -INC SMELEME
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC SMCOORD
  29. -INC CCHAMP
  30. -INC CCREEL
  31.       segment icpr(nbno)
  32.       segment trav1
  33.       character*4 compd(nbincl),compdd(nbincl)
  34.       endsegment
  35.       segment trav2
  36.       integer irigd(nbnoc,nbincd),nelrigd(nbnoc,nbincd)
  37.       endsegment
  38.       segment trav3
  39.       integer cor1p(nligrp),cor1d(nligrd)
  40.       integer imatt(nelrig)
  41.       endsegment
  42.       segment trav4
  43.       integer lisdes(klisd)
  44.       endsegment
  45.       segment trav5
  46.       real*8 xmatmpv(nligpv)
  47.       integer imatpv(nligpv)
  48.       endsegment
  49.       segment imatno(lmatno)
  50.       segment idesno(ldesno)
  51.       segment igeomno(lgeomno)
  52.       segment coerno(lcoerno)
  53.       segment irg7no(lrg7no)
  54.       segment irg6no(lrg6no)
  55.       segment irg5no(lrg5no)
  56.       segment irg2no(lrg2no)
  57.       logical dupli,descar,dupli4
  58.       character*4 co1,co1d
  59.       segment svu
  60.       logical vu(lrvu)
  61.       endsegment
  62.       integer lvu,symr
  63. *
  64. *      write (6,*) ' raideur initiale '
  65. *      call prrigi(mrigid,0)
  66. *      write (6,*) ' raideur dependance '
  67. *      call prrigi(ri1,1)
  68.  
  69. *  meme valeur de test que dans chole
  70.       xzref=5.d-12        
  71.       lrvu=100
  72.       lmatno=1000
  73.       ldesno=1000
  74.       lgeomno=1000
  75.       lcoerno=1000
  76.       lrg7no=1000
  77.       lrg6no=1000
  78.       lrg5no=1000
  79.       lrg2no=1000
  80.       kmatno=0
  81.       kdesno=0
  82.       kgeomno=0
  83.       kcoerno=0
  84.       krg7no=0
  85.       krg6no=0
  86.       krg5no=0
  87.       krg2no=0
  88.       segini imatno,idesno,igeomno,coerno,irg7no,irg5no,svu
  89.       lvu=0
  90.       segini irg2no,irg6no
  91.       nbdeg=0
  92. *  remplissage du segment trav1 (liste des noms de composantes)
  93.       nbincd=0
  94.       nbincl=16
  95.       segini trav1
  96.       nbno=nbpts
  97.       segact ri1
  98.       do 10 irig=1,ri1.irigel(/2)
  99.          descr=ri1.irigel(3,irig)
  100.          segact descr
  101.          do 40 nligrd=1,lisdua(/2)
  102.             do 50 inc=1,nbincd
  103.                if (compd(inc).eq.lisdua(nligrd)) goto 40
  104.  50         continue
  105.             nbincd=nbincd+1
  106.             if(nbincd.gt.nbincl) then
  107.              nbincl=2*nbincd
  108.              segadj trav1
  109.             endif
  110.             compd(nbincd)=lisdua(nligrd)
  111.             compdd(nbincd)=' '
  112. **          write(6,*) 'inc-1 compdd ',inc,compd(inc),compdd(inc)
  113.  40      continue
  114.          ldua=lisdua(/2)
  115.          do 45 nligrd=1,lisinc(/2)
  116.             do 55 inc=1,nbincd
  117.                if (compd(inc).eq.lisinc(nligrd)) goto 46
  118.  55         continue
  119.             nbincd=nbincd+1
  120.             inc=nbincd
  121.             if(nbincd.gt.nbincl) then
  122.              nbincl=2*nbincd
  123.              segadj trav1
  124.             endif
  125.             compd(nbincd)=lisinc(nligrd)
  126.  46      continue
  127.  
  128. ** si la matrice n'est pas carre, avec quoi on rempli?
  129.             compdd(inc)=lisdua(min(ldua,nligrd))
  130. **          write(6,*) 'inc-2 compdd ',inc,compd(inc),compdd(inc)
  131.  45      continue
  132.  10   continue
  133. *      write (6,*) ' scnd2 - 1 ',(compd(i),i=1,nbincd)
  134.       segact mrigid
  135.       do 110 irig=1,irigel(/2)
  136.          descr=irigel(3,irig)
  137.          segact descr
  138. **       write(6,*) 'nligrp,nligrd ',lisinc(/2),lisdua(/2)
  139.          ldua=lisdua(/2)
  140.          do 140 nligrd=1,lisinc(/2)
  141.             do 150 inc=1,nbincd
  142.                if (compd(inc).eq.lisinc(nligrd)) goto 141
  143.  150        continue
  144.             nbincd=nbincd+1
  145.             inc=nbincd
  146.             if(nbincd.gt.nbincl) then
  147.              nbincl=2*nbincd
  148.              segadj trav1
  149.             endif
  150.             compd(nbincd)=lisinc(nligrd)
  151.  141     continue
  152. *  au  cas ou la primale n'est pas dans nomdd, on sauve la duale
  153. **            write(6,*) 'inc lisdua descr',inc,lisdua(nligrd),descr
  154. ** si la matrice n'est pas carre, avec quoi on rempli la duale?
  155.             compdd(inc)=lisdua(min(ldua,nligrd))
  156. **          write(6,*) 'inc-3 compdd ',inc,compd(inc),compdd(inc)
  157.  140     continue
  158.  110  continue
  159. *      write (6,*) ' scnd2 - 2 ',(compd(i),i=1,nbincd)
  160. *  au tour de trav2 (inconnue => raideur de dépendance)
  161. *  mais d'abord icpr
  162.       segini icpr
  163.       nbnoc=0
  164.       do irig=1,ri1.irigel(/2)
  165.          meleme=ri1.irigel(1,irig)
  166.          descr=ri1.irigel(3,irig)
  167.          segact meleme,descr
  168.          do  j=1,num(/2)
  169.             ip=num(noeled(1),j)
  170.             if(icpr(ip).eq.0) then
  171.               nbnoc=nbnoc+1
  172.               icpr(ip)=nbnoc
  173.             endif
  174.          enddo
  175.       enddo
  176.       do irig=1,irigel(/2)
  177.          meleme=irigel(1,irig)
  178.          segact meleme
  179.          do  j=1,num(/2)
  180.          do  i=1,num(/1)
  181.             ip=num(i,j) 
  182.             if(icpr(ip).eq.0) then
  183.               nbnoc=nbnoc+1
  184.               icpr(ip)=nbnoc
  185.             endif
  186.          enddo
  187.          enddo
  188.       enddo
  189.       segini trav2
  190.       do 200 irig=1,ri1.irigel(/2)
  191.          meleme=ri1.irigel(1,irig)
  192.          descr=ri1.irigel(3,irig)
  193.          xmatri=ri1.irigel(4,irig)
  194.          segact xmatri
  195.          do 250 incd=1,nbincd
  196.             if (compd(incd).eq.lisdua(1)) goto 260
  197.  250     continue
  198.          call erreur(5)
  199.  260     continue
  200.          do 240 i=1,num(/2)
  201.             ip=num(noeled(1),i)
  202.             irigd(icpr(ip),incd)=irig
  203.             nelrigd(icpr(ip),incd)=i
  204.  240     continue
  205.  200  continue
  206. *  rechercher les noms d'inconnues duales
  207.       do 300 i=1,compd(/2)
  208.          do 310 j=1,lnomdd
  209.             if (compd(i).eq.nomdd(j)) goto 320
  210.  310     continue
  211. **          write (6,*) ' pas de nom dual trouvé ',compd(i),compdd(i)
  212.          goto 300
  213.  320     continue
  214.             compdd(i)=nomdu(j)
  215. **        write (6,*) ' correspondance ', compd(i), compdd(i)
  216.  300  continue
  217. *        write (6,*) ' scnd2 - 3 ',compdd(/2),(compdd(i),i=1,compdd(/2))
  218. *
  219. *  maintenant, on balaye les raideurs. Pour chacune d'entre elle, on
  220. *  balaye les inconnues pour trouver une dépendance. Si oui on construit
  221. *  le résultat de l'opération
  222.       segini,ri4=mrigid
  223. *     write(6,*) 'ri4 mrigid ',ri4,mrigid
  224.       ri4.ichole=0
  225.       ri4.imgeo1=0
  226.       ri4.imgeo2=0
  227.       ri4.isupeq=0
  228. *  boucle raideur en cours de condensation
  229. *  on range dans des2 le descripteur precedent pour le reutiliser si possible
  230.       des2=0
  231.       isup=0
  232.       klisd=100
  233.       llisd=0
  234.       segini trav4
  235. *     call prrigi(ri4,1)
  236.       do 1000 irig=1,ri4.irigel(/2)
  237. *         write(ioimp,*) 'Balayage raideur ',irig
  238. ** on ne dupliquera meleme et xmatri que si necessaire
  239.          meleme=ri4.irigel(1,irig)
  240.          ipt4=meleme
  241.          dupli4=.false.
  242.          ri4.irigel(1,irig)=ipt4
  243.          des4=ri4.irigel(3,irig)
  244.          segact des4
  245.          Xmatri=ri4.irigel(4,irig)
  246.          segact xmatri
  247. *  au cas ou symre ait ete oublie
  248.          if (symre.ne.ri4.irigel(7,irig)) then
  249.            write(6,*) 'incoherence symre irigel(7) pour ',ri4,xmatri
  250.            segact xmatri*mod
  251.            symre = max(symre,ri4.irigel(7,irig))
  252.            ri4.irigel(7,irig)=symre
  253.          endif
  254.          Xmatr4=Xmatri
  255. *        write(6,*) 'scnd2 xmatri xmatr4 dims ',xmatri,xmatr4,
  256. *    >     xmatr4.re(/1),xmatr4.re(/2),xmatr4.re(/3)
  257.          nelrig=xmatr4.re(/3)
  258.          ri4.irigel(4,irig)=Xmatr4
  259. * On regarde le descripteur. S'il est symétrique, on pourra
  260. * faire le travail en une passe sinon deux.
  261. *   Correspondance lisinc, lisdua, numéro incd dans trav1
  262.          nligrp=des4.lisinc(/2)
  263.          nligrd=des4.lisdua(/2)
  264.          segini trav3
  265.          nligpv=nligrp*2+256
  266.          segini trav5
  267.          do 410 iligrp=1,nligrp
  268.             do 420 incd=1,nbincd
  269.                if (compd(incd).eq.des4.lisinc(iligrp)) then
  270.                   cor1p(iligrp)=incd
  271.                   goto 410
  272.                endif
  273.  420        continue
  274. **          write(6,*) 'cor1p pas de correspondance'
  275. **          call erreur(5)
  276.  410     continue
  277.          do 510 iligrd=1,nligrd
  278.             do 520 incd=1,nbincd
  279. **             write(6,*) 'compdd ',compdd(incd),des4.lisdua(iligrd)
  280.                if (compdd(incd).eq.des4.lisdua(iligrd)) then
  281.                   cor1d(iligrd)=incd
  282.                   goto 510
  283.                endif
  284.  520        continue
  285. **          write(6,*) 'cor1d pas de correspondance descar',descar
  286. **          call erreur(5)
  287.  510     continue
  288. *         write(ioimp,*) 'nligrp=',nligrp
  289. *         write(ioimp,*) 'nligrd=',nligrd
  290.          descar=.false.
  291.          if (nligrp.ne.nligrd) goto 600
  292.          do iligr=1,nligrp
  293.             if (cor1p(iligr).ne.cor1d(iligr)) goto 600
  294.             if (des4.noelep(iligr).ne.des4.noeled(iligr)) goto 600
  295.          enddo
  296.          if (xmatr4.symre.eq.0) descar=.true.
  297.  600     continue
  298. *         WRITE(IOIMP,194)(des4.NOELEP(J),des4.LISINC(J),cor1p(j),J=1
  299. *     $        ,nligrp)
  300. *         WRITE(IOIMP,194)(des4.NOELED(J),des4.LISDUA(J),cor1d(j),J=1
  301. *     $        ,nligrd)
  302. *  194     FORMAT( I6,9X,A4,9X,I6)
  303. *         if (descar) then
  304. *            write(ioimp,*) 'le descripteur est carre'
  305. *         else
  306. *            write(ioimp,*) 'le descripteur n''est pas carre'
  307. *         endif
  308.          npass=2
  309.          idsup=0
  310.          do 1040 i=1,ipt4.num(/2)
  311.             dupli=.false.
  312.             renorm=0.d0
  313. *    On fait deux passes : une pour la primale et une pour la duale
  314.             do 1045 ipass=1,npass
  315.                if (ipass.eq.1) then
  316.                   nligr4=des4.lisinc(/2)
  317.                else
  318.                   nligr4=des4.lisdua(/2)
  319.                endif
  320.                do 1010 nligri=1,nligr4
  321.                   if (ipass.eq.1) then
  322.                      incd=cor1p(nligri)
  323.                   else
  324.                      incd=cor1d(nligri)
  325.                   endif
  326. ***               write(6,*) 'ipass incd ',ipass,incd
  327.                   if (incd.eq.0) goto 1010
  328.                   if (ipass.eq.1) then
  329.                      ip=ipt4.num(des4.noelep(nligri),i)
  330.                   else
  331.                      ip=ipt4.num(des4.noeled(nligri),i)
  332.                   endif
  333. *  y a t'il une dependance concernant ce noeud
  334. *        write(6,*) ' descar',descar
  335. ***                 write(6,*) 'ipass ip icpr incd nelrig',ipass,ip,
  336. **   >              icpr(ip),incd,nelrigd(icpr(ip),incd)
  337.                   if (nelrigd(icpr(ip),incd).eq.0) then
  338.                     goto 1010
  339.                   endif
  340. *  ici il faut eliminer l'inconnue
  341. *  on commence par dupliquer la la matrice si ce n'est pas deja fait
  342.                   if (.not.dupli) then
  343.                      ird=irigel(/1)
  344.                      krg7no=krg7no+1
  345.                      if (krg7no.gt.irg7no(/1)) then
  346.                        lrg7no=irg7no(/1)+1000
  347.                        segadj irg7no
  348.                      endif
  349.                      if (ird.ge.7) then
  350.                         irg7no(krg7no)=ri4.irigel(7,irig)
  351.                      else
  352.                         irg7no(krg7no)=0
  353.                      endif
  354.                      krg6no=krg6no+1
  355.                      if (krg6no.gt.irg6no(/1)) then
  356.                        lrg6no=irg6no(/1)+1000
  357.                        segadj irg6no
  358.                      endif
  359.                      if (ird.ge.6) then
  360.                         irg6no(krg6no)=ri4.irigel(6,irig)
  361.                      else
  362.                         irg6no(krg6no)=0
  363.                      endif
  364.                      krg5no=krg5no+1
  365.                      if (krg5no.gt.irg5no(/1)) then
  366.                        lrg5no=irg5no(/1)+1000
  367.                        segadj irg5no
  368.                      endif
  369.                      if (ird.ge.5) then
  370.                         irg5no(krg5no)=ri4.irigel(5,irig)
  371.                      else
  372.                         irg5no(krg5no)=0
  373.                      endif
  374.                      krg2no=krg2no+1
  375.                      if (krg2no.gt.irg2no(/1)) then
  376.                        lrg2no=irg2no(/1)+1000
  377.                        segadj irg2no
  378.                      endif
  379.                      irg2no(krg2no)=ri4.irigel(2,irig)
  380.                      kdesno=kdesno+1
  381.                      if (kdesno.gt.idesno(/1)) then
  382.                        ldesno=idesno(/1)+1000
  383.                        segadj idesno
  384.                      endif
  385.                      nbnn=ipt4.num(/1)
  386.                      nbelem=1
  387.                      nbref=0
  388.                      nbsous=0
  389.                      segini meleme
  390.                      itypel=28
  391. *  si c'est une matrice de relation, il faut garder l'attribut
  392. **                   if(ipt4.itypel.eq.49) write(6,*)
  393. **   >                  'ipt4 dans scnd2 ',ipt4
  394.                      if (ipt4.itypel.eq.49) itypel=49
  395.                      if (ipt4.itypel.eq.22) itypel=22
  396.                      do i1=1,nbnn
  397.                         num(i1,1)=ipt4.num(i1,i)
  398.                      enddo
  399.                      kgeomno=kgeomno+1
  400.                      if (kgeomno.gt.igeomno(/1)) then
  401.                        lgeomno=igeomno(/1)+1000
  402.                        segadj igeomno
  403.                      endif
  404.                      igeomno(kgeomno)=meleme
  405.                      kligrp=xmatr4.re(/2)
  406.                      kligrd=xmatr4.re(/1)
  407.                   irig1=irigd(icpr(ip),incd)
  408.                   des1=ri1.irigel(3,irig1)
  409.                   segact des1
  410. *  on agrandit xmatri du nb d'inc. primales presentent dans la relation
  411.                   nligrp1=max(des1.lisinc(/2),0)
  412.                      nligrp=kligrp+nligrp1
  413.                      nligrd=kligrd+nligrp1
  414.                      nelrig=1
  415.                      rigrel=0
  416.                      segini xmatri
  417.                      xmatri.symre=xmatr4.symre
  418. *  re est plus grand que xmatr4.re donc pas d'appel a scndc possible
  419.                      renorm=0.d0
  420.                      ii=1
  421.                      if(itypel.eq.22) ii=2
  422.                      if(itypel.eq.49) ii=3
  423.                      do io=1,kligrp
  424.                         do iu=1,kligrd
  425.                            retmp=xmatr4.re(iu,io,i)
  426.                            re(iu,io,1)=retmp
  427.                            if (io.ge.ii.or.iu.ge.ii) renorm=renorm+
  428.      >                                                abs(retmp)
  429.                         enddo
  430.                      enddo
  431.                      renorm=max(xpetit,renorm*xzref)
  432.                      segini,descr
  433.                      idesno(kdesno)=descr
  434.                      do io=1,des4.noelep(/1)
  435.                         noelep(io)=des4.noelep(io)
  436.                         lisinc(io)=des4.lisinc(io)
  437.                      enddo
  438.                      do io=1,des4.noeled(/1)
  439.                         noeled(io)=des4.noeled(io)
  440.                         lisdua(io)=des4.lisdua(io)
  441.                      enddo
  442.  
  443.                      kmatno=kmatno+1
  444.                      if (kmatno.gt.imatno(/1))then
  445.                        lmatno=imatno(/1)+1000
  446.                        segadj imatno
  447.                      endif
  448.                      imatno(kmatno)=xmatri
  449.                      kcoerno=kcoerno+1
  450.                      if (kcoerno.gt.coerno(/1))then
  451.                        lcoerno=coerno(/1)+1000
  452.                        segadj coerno
  453.                      endif
  454.                      coerno(kcoerno)=ri4.coerig(irig)
  455.                      dupli=.true.
  456.                      imatt(i)=1
  457.                   endif
  458. *  on elimine dans la matrice dupliquee
  459.                   irig1=irigd(icpr(ip),incd)
  460.                   i1=nelrigd(icpr(ip),incd)
  461.                   ipt1=ri1.irigel(1,irig1)
  462.                   des1=ri1.irigel(3,irig1)
  463.                   xmatr1=ri1.irigel(4,irig1)
  464.                   segact ipt1,des1,xmatr1
  465. *          nligrd1 =des1.lisdua(/2)
  466. *   nligrd1=1 est l'inconnue à éliminer (de la matrice de dependance)
  467.                   nligrd1=1
  468.                   nligrp1=des1.lisinc(/2)
  469. *          write (6,*) ' scnd2  nligrp1  nligrd1 i1 ',nligrp1,nligrd1,i1,
  470. *     > ip,incp
  471. *  nlig1 : inconnue primale à remettre dans la matrice (comme primale)
  472.                   ipvm=0
  473.                   do 1050 nlig1=1,nligrp1
  474.                      xmatmp = xmatr1.re(1,nlig1,i1)
  475.                      ip1=ipt1.num(des1.noelep(nlig1),i1)
  476.                      co1=des1.lisinc(nlig1)
  477.                      if (ipass.eq.2.or.descar) then
  478.                         do  incd=1,nbincd
  479.                            if (compd(incd).eq.co1) then
  480.                               co1d=compdd(incd)
  481.                               goto 1055
  482.                            endif
  483.                         enddo
  484.                         write (6,*) ' scnd2 2 pas de correspondance '
  485.      $                       ,co1
  486.                         call erreur(5)
  487.                         return
  488.  1055                   continue
  489.                      endif
  490. *           write (6,*) ' ip1 co1 nligrpi ',ip1,co1,nligrpi
  491. *  le ddl est il deja dans la raideur
  492.                      nbno=0
  493.                      if (ipass.eq.1) then
  494.                         do 1060 nligr=1,kligrp
  495.                            if (noelep(nligr).eq.0) goto 1060
  496.                            if (num(noelep(nligr),1).ne.ip1) goto 1060
  497. * on note le numéro local du noeud pour le reutiliser
  498.                            nbno=noelep(nligr)
  499.                            if (lisinc(nligr).ne.co1) goto 1060
  500.                            goto 1070
  501.  1060                   continue
  502.                      else
  503.                         do 2060 nligr=1,kligrd
  504.                            if (noeled(nligr).eq.0) goto 2060
  505.                            if (num(noeled(nligr),1).ne.ip1) goto 2060
  506. * on note le numéro local du noeud pour le reutiliser
  507.                            nbno=noeled(nligr)
  508.                            if (lisdua(nligr).ne.co1d) goto 2060
  509.                            goto 1070
  510.  2060                   continue
  511.                      endif
  512. *  on n'a pas trouve le ddl, on le rajoute
  513. *            write (6,*) ' ajout de ddl'
  514.                      if (descar) then
  515.                         nligr =kligrp+1
  516.                         kligrp=kligrp+1
  517.                         kligrd=kligrd+1
  518.                      else
  519.                         if (ipass.eq.1) then
  520.                            nligr =kligrp+1
  521.                            kligrp=kligrp+1
  522.                            kligrd=kligrd
  523.                            xmatri.symre=2
  524.                         else
  525.                            nligr =kligrd+1
  526.                            kligrp=kligrp
  527.                            kligrd=kligrd+1
  528.                            xmatri.symre=2
  529.                         endif
  530.                      endif
  531.                      nligrp=lisinc(/2)
  532.                      nligrd=lisdua(/2)
  533.                      if (kligrp.gt.nligrp.or.kligrd.gt.nligrd) then
  534.                        nligrd=kligrd+max(1,nligrp1)
  535.                        nligrp=kligrp+max(1,nligrp1)
  536. *                      write(6,*) ' premier adj xmatri '
  537.                        rigrel=0
  538.                        segadj xmatri
  539.                        segadj descr
  540.                      endif
  541.                      if (nbno.eq.0) then
  542.                         nbsous=0
  543.                         nbelem=1
  544.                         nbnn=num(/1)+1
  545.                         nbref=0
  546.                         segadj meleme
  547.                         nbno=nbnn
  548.                      endif
  549.                      num(nbno,1)=ip1
  550.                      if (ipass.eq.1.or.descar) then
  551.                         noelep(kligrp)=nbno
  552.                         lisinc(kligrp)=co1
  553.                      endif
  554.                      if (ipass.eq.2.or.descar) then
  555.                         lisdua(kligrd)=co1d
  556.                         noeled(kligrd)=nbno
  557.                      endif
  558.  1070                continue
  559. *  on met a jour la raideur
  560.                      nligrp=kligrp
  561.                      nligrd=kligrd
  562. *            segact xmatr1
  563.                         if(nligr.gt.nligpv) then
  564.                           nligpv=nligr*2
  565.                           segadj trav5
  566. **                        write(6,*) 'agrandissement trav5 ',nligpv
  567.                         endif
  568.                      if (ipass.eq.2) then
  569.                         ipvm=ipvm+1
  570.                         imatpv(ipvm)=nligr
  571.                         xmatmpv(ipvm)=xmatmp
  572.                      endif
  573.                      if (ipass.eq.1) then
  574. *  unrolling sur 2
  575.                         do j=1,nligrd-1,2
  576.                            retmp1=re(j,nligr,1)+re(j,nligri,1)
  577.      $                          *xmatmp
  578.                            re(j,nligr,1)=retmp1
  579.                            retmp2=re(j+1,nligr,1)+re(j+1,nligri,1)
  580.      $                          *xmatmp
  581.                            re(j+1,nligr,1)=retmp2
  582.                         enddo
  583.                         if (j.eq.nligrd) then
  584.                            retmp1=re(j,nligr,1)+re(j,nligri,1)
  585.      $                          *xmatmp
  586.                            re(j,nligr,1)=retmp1
  587.                         endif
  588.                      endif
  589. *    la mise a jour est sorti de la boucle pour ppouvoir l'effectuer dans l'ordre efficace pour le cache
  590. ***                  if (ipass.eq.2) then
  591. ***                     do j=1,nligrp
  592. ***                        re(nligr,j,1)=re(nligr,j,1)+re(nligri,j,1)
  593. ***  $                          *xmatmp
  594. ***                   enddo
  595. ***                endif
  596. *              segdes xmatr1,des1
  597.  1050             continue
  598. *  mise a jour transferee ici
  599.                  if(ipass.eq.2) then
  600. *  unrolling sur 2
  601.                   do j=1,nligrp
  602.                    retmp=re(nligri,j,1)
  603.                    do ipv=1,ipvm-1,2    
  604.                     nligr1=imatpv(ipv)
  605.                     re(nligr1,j,1)=re(nligr1,j,1)+retmp*xmatmpv(ipv)
  606.                     nligr2=imatpv(ipv+1)
  607.                     re(nligr2,j,1)=re(nligr2,j,1)+retmp*xmatmpv(ipv+1)
  608.                    enddo
  609.                    if (ipv.eq.ipvm) then
  610.                        nligr=imatpv(ipv)
  611.                        re(nligr,j,1)=re(nligr,j,1)+retmp*xmatmpv(ipv)
  612.                    endif
  613.                   enddo
  614.                  endif
  615. *  on note les inconnues à supprimer (primale et duale)
  616. *           if (ipass.eq.2.or.descar)
  617. *    >      write (6,*) ' il faut eliminer ',nligri,lisdua(nligri)
  618. *      write (6,*) ' lisinc contient ',(lisinc(iou),iou=1,lisinc(/2))
  619.                   if (ipass.eq.1.or.descar) noelep(nligri)=0
  620.                   if (ipass.eq.2.or.descar) noeled(nligri)=0
  621.  1010          continue
  622. *
  623. *            write(ioimp,*) 'Avant suppression des inconnues flaguees'
  624. *          WRITE(IOIMP,194)(NOELEP(J),LISINC(J),J=1,LISINC(/2))
  625. *          WRITE(IOIMP,194)(NOELED(J),LISDUA(J),J=1,LISDUA(/2))
  626. *  194     FORMAT( I6,9X,A4)
  627. *          WRITE(IOIMP,195) ((RE(L,K,1),K=1,RE(/2)),L=1,RE(/1))
  628. * 195      FORMAT(1X,6E12.5)
  629. *  on supprime les inconnues flaguées
  630.  1045       continue
  631.             if (dupli) then
  632.                nligrd=kligrd
  633.                nligrp=kligrp
  634.                idecp=0
  635.                do 1080 nligr=1,nligrp
  636.                   if (noelep(nligr).eq.0) then
  637.                      idecp=idecp+1
  638.                   elseif (idecp.ne.0) then
  639.                      lisinc(nligr-idecp)=lisinc(nligr)
  640.                      noelep(nligr-idecp)=noelep(nligr)
  641.                      do nligd=1,nligrd
  642.                         if (noeled(nligd).ne.0)
  643.      >                  re(nligd,nligr-idecp,1)=re(nligd,nligr,1)
  644.                      enddo
  645.                   endif
  646.  1080          continue
  647. *pv  on imbrique la boucle 1093 dans ce sens pour minimiser l'utilisation du cache
  648.                nligrp=nligrp-idecp
  649.              do 1093 nligr=1,nligrp
  650.                idecd=0
  651.                    do nligd=1,nligrd
  652.                 if (noeled(nligd).eq.0) then
  653.                    idecd=idecd+1
  654.                 elseif (idecd.ne.0) then
  655.                       re(nligd-idecd,nligr,1)=re(nligd,nligr,1)
  656.                 endif
  657.                    enddo
  658. 1093          continue
  659.              idecd=0
  660.              do 1092 nligd=1,nligrd
  661.                 if (noeled(nligd).eq.0) then
  662.                    idecd=idecd+1
  663.                 elseif (idecd.ne.0) then
  664.                    lisdua(nligd-idecd)=lisdua(nligd)
  665.                    noeled(nligd-idecd)=noeled(nligd)
  666.                 endif
  667. 1092         continue
  668.                nligrd=nligrd-idecd
  669. *            write (6,*) ' scnd2 avant apres   ',lisinc(/2),nligrp
  670.           if(descar.and..true.) then
  671. *  symetriser xmatri
  672.           ninc=32
  673.           do ib=1,nligrd,ninc
  674.             do jr=1,nligrd 
  675.               do ir=ib,min(jr-1,ib+ninc-1)
  676.                 re1=re(ir,jr,1)
  677.                 re2=re(jr,ir,1)
  678.                    rem=(re1+re2)/2.d0
  679.                    re(ir,jr,1)=rem
  680.                    re(jr,ir,1)=rem
  681.               enddo
  682.             enddo
  683.           enddo
  684.           endif
  685.  
  686. *  on elimine les matrices ne contenant que les multiplicateurs de Lagrange
  687. *  ou les autres termes non significatifs  
  688. *  la verif du second membre nul est faite dans vechpo
  689. *
  690. **             if(itypel.eq.49) write(6,*) 'scnd2 49 ',meleme
  691.                if (itypel.eq.49.or.itypel.eq.22) then
  692. ***          write (6,*)  (lisinc(iou),iou=1,lisinc(/2))
  693.                    iid=2
  694.                    if (itypel.eq.49) iid=3
  695. *  suppression des termes negligeables de la relation
  696.                    idec=0
  697. *  on utilise une valeur absolu car la relation initiale est supposee normalisee a 1
  698.                    do ii=iid,nligrd
  699.                      if (abs(re(ii,1,1)).lt.renorm.and.
  700.      >                   abs(re(1,ii,1)).lt.renorm) then
  701.                       idec=idec+1
  702.                      elseif(idec.ne.0) then
  703.                         re(ii-idec,1,1)=re(ii,1,1)
  704.                         if(iid.eq.3) re(ii-idec,2,1)=re(ii,2,1)
  705.                         re(1,ii-idec,1)=re(1,ii,1)
  706.                         if(iid.eq.3) re(2,ii-idec,1)=re(2,ii,1)
  707.                         lisinc(ii-idec)=lisinc(ii)
  708.                         lisdua(ii-idec)=lisdua(ii)
  709.                         noelep(ii-idec)=noelep(ii)
  710.                         noeled(ii-idec)=noeled(ii)
  711.                      endif
  712.                    enddo
  713.                    nligrp=nligrp-idec
  714.                    nligrd=nligrd-idec
  715.                    if(nligrd.lt.iid) then
  716. **           write (6,*) 'scnd2 matrice supprimee LX',
  717. **   >       nligrd,itypel,meleme
  718.                      nligrd=0
  719.                      nligrp=0
  720.                      nbdeg=nbdeg+1
  721.                    endif
  722.                endif
  723.                if (re(/1).ne.nligrd.or.re(/2).ne.nligrp) then
  724. **             write(6,*) ' ajustement final xmatri ',re(/1),re(/2),
  725. **   >         nligrd,nligrp
  726.                rigrel=0
  727.                segadj xmatri
  728.                else
  729. **              write(6,*) ' pas d ajustement xmatri '
  730.                endif
  731.                if (nligrp.ne.nligrd.or..not.descar) then
  732.                  xmatri.symre=2
  733.                else
  734.                endif
  735.                segdes xmatri
  736.                segadj descr
  737. * la matrice est fabriquée, on elimine maintenant les noeuds
  738. * non references
  739.                lvu=num(/1)
  740.                if (lvu.gt.lrvu) then
  741.                 lrvu=lvu+100
  742.                 segadj svu
  743.                endif
  744.                do in=1,lvu
  745.                   vu(in)=.false.
  746.                enddo
  747.                do in=1,nligrp
  748.                   vu(noelep(in))=.true.
  749.                enddo
  750.                if (.not.descar) then
  751.                   do in=1,nligrd
  752.                      vu(noeled(in))=.true.
  753.                   enddo
  754.                endif
  755.                idec=0
  756.                do 1095 in=1,num(/1)
  757.                   if (.not.vu(in)) then
  758. *  noeud a eliminer
  759.                      idec=idec+1
  760.                   elseif (idec.ne.0) then
  761.                      do no=1,nligrp
  762.                         if (noelep(no).eq.in) then
  763.                            noelep(no)=noelep(no)-idec
  764.                            if (descar) noeled(no)=noeled(no)-idec
  765.                         endif
  766.                      enddo
  767.                      if (.not.descar) then
  768.                         do no=1,nligrd
  769.                            if (noeled(no).eq.in) then
  770.                               noeled(no)=noeled(no)-idec
  771.                            endif
  772.                         enddo
  773.                      endif
  774.                      num(in-idec,1)=num(in,1)
  775.                   endif
  776.  1095          continue
  777.                nbnn=num(/1)-idec
  778.                nbelem=1
  779.                nbsous=0
  780.                nbref=0
  781.                segadj meleme
  782.                if (nbnn.eq.0) imatno(kmatno)=0
  783. *  recyclage eventuel du descripteur
  784.                do 4032 ides=llisd,1,-1
  785.                des2=lisdes(ides)
  786.                segact des2
  787.          if (nligrd.ne.des2.noeled(/1)) goto 4030
  788.          if (nligrp.ne.des2.noelep(/1)) goto 4030
  789.          do id1=1,nligrd
  790.            if (noeled(id1).ne.des2.noeled(id1)) goto 4030
  791.            if (lisdua(id1).ne.des2.lisdua(id1)) goto 4030
  792.          enddo
  793.          do id1=1,nligrp
  794.            if (noelep(id1).ne.des2.noelep(id1)) goto 4030
  795.            if (lisinc(id1).ne.des2.lisinc(id1)) goto 4030
  796.          enddo
  797.          segsup descr
  798.          idsup=idsup+1
  799.          descr=des2
  800.          idesno(kdesno)=descr
  801.          goto 4031
  802.  4030    continue
  803.  4032    continue
  804.                 llisd=llisd+1
  805.                 if (llisd.gt.klisd) then
  806.                   klisd=llisd+100
  807.                   segadj trav4
  808.                 endif
  809.                 lisdes(llisd)=descr
  810. *               write(6,*) ' nouveau descripteur ',nligrp,nligrd
  811. **              segdes des2
  812.  4031    continue
  813. **             endif
  814.                if (nligrp.ne.noelep(/1).or.nligrd.ne.noeled(/1))
  815.      >             segadj descr
  816.             endif
  817.  1040    continue
  818. *  compression de ipt4 et imatr4
  819.          idec=0
  820.          ipt5=ipt4
  821.          xmatr5=xmatr4
  822.          nligrp=xmatr5.re(/2)
  823.          nligrd=xmatr5.re(/1)
  824.          nelrig=xmatr5.re(/3)
  825.          nbnn=ipt5.num(/1)
  826.          nbelem=ipt5.num(/2)
  827.          idec=0
  828.          do i=1,xmatr4.re(/3)
  829.            if (imatt(i).ne.0) idec=idec+1
  830.          enddo
  831. *  on duplique systematiquement pour pouvoir detruire dans resouc
  832.          if (idec.ne.0.or..true.) then
  833.                  nbelem=nbelem-idec
  834.                  nbsous=0
  835.                  nbref=ipt5.lisref(/1)
  836.                  segini,ipt4
  837.                  ipt4.itypel=ipt5.itypel
  838.                  do ir=1,ipt4.lisref(/1)
  839.                   ipt4.lisref(ir)=ipt5.lisref(ir)
  840.                  enddo
  841.                  ri4.irigel(1,irig)=ipt4
  842.                  nelrig=nelrig-idec
  843.                  segini,xmatr4
  844.                  ri4.irigel(4,irig)=Xmatr4
  845.                  xmatr4.symre=xmatr5.symre
  846.                  idec=0
  847.          endif
  848.          do 2000 i=1,xmatr5.re(/3)
  849.             if(imatt(i).ne.0) then
  850.                idec=idec+1
  851.             elseif (xmatr4.ne.xmatr5) then
  852. **             do io=1,nligrp
  853. **                do iu=1,nligrd
  854. **                   xmatr4.re(iu,io,i-idec)=xmatr5.re(iu,io,i)
  855. **                enddo
  856. **             enddo
  857.                      call scndc(nligrp*nligrd,xmatr5.re(1,1,i),
  858.      >               xmatr4.re(1,1,i-idec))
  859.                do in=1,nbnn
  860.                    ipt4.num(in,i-idec)=ipt5.num(in,i)
  861.                enddo
  862.             endif
  863.  2000    continue
  864.          segdes xmatr4
  865. *        write (6,*) ' nelrig dans scnd2 ',nelrig
  866.          if (nelrig.eq.0) ri4.irigel(4,irig)=0
  867.          segsup trav3
  868.          segsup trav5
  869. **       segdes descr
  870.  1000 continue
  871.  
  872. *
  873. *  compression des xmatri suivant les descripteurs
  874. *
  875.        isup=0
  876.        id=0
  877.        is=1
  878.        nid=1
  879.        nrigf=0
  880.  4100  id=nid
  881.        if (id.gt.kdesno) goto 4120
  882.        xmatri=imatno(id)
  883.        if (xmatri.eq.0) then
  884. *        write(6,*) ' xmatri 0',id
  885.          isup=isup+1
  886.          nid=id+1
  887.          goto 4100
  888.        endif
  889.        meleme=igeomno(id)
  890.        segact meleme
  891.        des4=idesno(id)
  892.        nid=0
  893.        nelrig=num(/2) 
  894. *      if (num(/2).ne.1) write(6,*) 'nelrig-1',nelrig
  895.        nelrii=nelrig
  896.        nrigf=nrigf+1
  897.        lis=id
  898.        do is=id+1,kdesno
  899.          if (idesno(is).ne.des4) goto 4110
  900.          if (imatno(is).eq.0) goto 4110
  901.          if (irg2no(is).ne.irg2no(id)) goto 4110
  902.          if (irg5no(is).ne.irg5no(id)) goto 4110
  903.          if (irg6no(is).ne.irg6no(id)) goto 4110
  904.          if (irg7no(is).ne.irg7no(id)) goto 4110
  905.          lis=is
  906.          meleme=igeomno(is)
  907.          segact meleme
  908.          nelrig=nelrig+num(/2)
  909. *      if (num(/2).ne.1) write(6,*) 'nelrig-2',num(/2)
  910.          goto 4112
  911.  4110   continue
  912.         if (nid.eq.0) nid=is
  913.  4112   continue
  914.        enddo
  915.        if (nid.eq.0) nid=is
  916.        if (nelrig.eq.nelrii) goto 4100
  917.        xmatri=imatno(id)
  918.        segact xmatri
  919.        symr=symre
  920.        nligrd=re(/1)
  921.        nligrp=re(/2)
  922.        rigrel=0
  923.        segini xmatri
  924.        symre=symr
  925.        meleme=igeomno(id)
  926.        segact meleme
  927.        nbnn=num(/1)
  928.        nbelem=nelrig
  929.        nbsous=0
  930.        nbref=0
  931.        itype=itypel
  932.        segini meleme
  933.        itypel=itype
  934.        nelrig=0
  935.        do is=id,lis
  936.         coeref=1.d0
  937.         if (is.ne.id) then
  938.          if (idesno(is).ne.des4) goto 4111
  939.          if (imatno(is).eq.0) goto 4111
  940.          if (irg2no(is).ne.irg2no(id)) goto 4111
  941.          if (irg5no(is).ne.irg5no(id)) goto 4111
  942.          if (irg6no(is).ne.irg6no(id)) goto 4111
  943.          if (irg7no(is).ne.irg7no(id)) goto 4111
  944.          coeref=coerno(is)/coerno(id)
  945.         endif
  946.          xmatr1=imatno(is)
  947.          ipt1=igeomno(is)
  948.          segact xmatr1,ipt1
  949.          symre=max(symre,xmatr1.symre)
  950.          do nel=1,xmatr1.re(/3)
  951.           nelrig=nelrig+1
  952.           do irp=1,nligrp
  953.           do ird=1,nligrd
  954.            re(ird,irp,nelrig)=xmatr1.re(ird,irp,nel)*coeref
  955.           enddo
  956.           enddo
  957.           do ip=1,nbnn
  958.            num(ip,nelrig)=ipt1.num(ip,nel)
  959.           enddo
  960.           imatno(is)=0
  961.          enddo
  962.         segsup xmatr1,ipt1
  963.         isup=isup+1
  964.  4111   continue
  965.        enddo
  966.        imatno(id)=xmatri
  967.        igeomno(id)=meleme
  968.        segdes xmatri
  969.  
  970.        goto 4100
  971.  4120  continue
  972. *      write(6,*) ' suppression de raideurs et descr ',isup,idsup
  973. *      write(6,*) ' nb desc          ',llisd
  974. **     write(6,*) ' lisdes ',(lisdes(il),il=1,lisdes(/1))
  975.          segsup trav4
  976. *  compression de irigel
  977.       idec=0
  978.  
  979.       do 2100 irig=1,ri4.irigel(/2)
  980.          if (ri4.irigel(4,irig).eq.0.or.abs(coerig(irig)).lt.xpetit) 
  981.      >      then
  982.             idec=idec+1
  983.          elseif (idec.ne.0) then
  984.             do in=1,ri4.irigel(/1)
  985.                ri4.irigel(in,irig-idec)=ri4.irigel(in,irig)
  986.             enddo
  987.             ri4.coerig(irig-idec)=ri4.coerig(irig)
  988.          endif
  989.  2100 continue
  990.       nrigeo=ri4.irigel(/2)-idec
  991. *  il ne reste plus qu'a ranger dans ri4 les nouvelles raideurs engendrées
  992.       nrigel=nrigf+nrigeo
  993.       segadj ri4
  994.       ir=nrigeo
  995.       do 3000 irig=1,kmatno
  996.          if (imatno(irig).eq.0) goto 3000
  997.          if (abs(coerno(irig)).lt.xpetit) goto 3000
  998.          ir=ir+1
  999.          ri4.coerig(ir)=coerno(irig)
  1000.          ri4.irigel(1,ir)=igeomno(irig)
  1001.          ri4.irigel(3,ir)=idesno(irig)
  1002.          ri4.irigel(7,ir)=irg7no(irig)
  1003.          ri4.irigel(5,ir)=irg5no(irig)
  1004.          ri4.irigel(2,ir)=irg2no(irig)
  1005.          ri4.irigel(6,ir)=irg6no(irig)
  1006.          ri4.irigel(4,ir)=imatno(irig)
  1007.  3000 continue
  1008. *     write(6,*) ' ir nrigel ',ir,nrigel,nrigeo
  1009. *  il peut y avoir des imatri nul a enlever
  1010.       if (ir.lt.nrigel) then
  1011.        nrigel=ir
  1012.        segadj ri4
  1013.       endif
  1014.       if (ir.ne.nrigel) call erreur(5)
  1015. *     write (6,*) ' raideur resultante '
  1016. *       call prrigi(ri4,0)
  1017.       segsup trav1,trav2,imatno,idesno,igeomno,coerno,irg7no,icpr
  1018.      $     ,irg5no,irg2no,irg6no,svu
  1019. *  desactivation generale
  1020.       segact mrigid
  1021.       do irig=1,irigel(/2)
  1022.          xmatri=irigel(4,irig)
  1023.          meleme=irigel(1,irig)
  1024.          descr=irigel(3,irig)
  1025.          segdes xmatri,descr
  1026.       enddo
  1027.       segact ri4
  1028.       do irig=1,ri4.irigel(/2)
  1029.          xmatri=ri4.irigel(4,irig)
  1030.          meleme=ri4.irigel(1,irig)
  1031.          descr=ri4.irigel(3,irig)
  1032.          segdes xmatri,descr
  1033.       enddo
  1034.       segact ri1
  1035.       do irig=1,ri1.irigel(/2)
  1036.          xmatri=ri1.irigel(4,irig)
  1037.          meleme=ri1.irigel(1,irig)
  1038.          descr=ri1.irigel(3,irig)
  1039.          segdes xmatri,descr
  1040.       enddo
  1041.       if(iimpi.ne.0.and.nbdeg.ne.0)
  1042.      >  write (6,*) ' nombre de relations degenerees'//
  1043.      >  ' apres elimination ',nbdeg
  1044.       end
  1045.  
  1046.  
  1047.  
  1048.  
  1049.  
  1050.  
  1051.  
  1052.  
  1053.  
  1054.  
  1055.  
  1056.  
  1057.  
  1058.  
  1059.  
  1060.  
  1061.  
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070.  
  1071.  
  1072.  
  1073.  
  1074.  
  1075.  
  1076.  
  1077.  
  1078.  
  1079.  
  1080.  
  1081.  

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