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    MB234859  25/01/03    21:15:30     12105          
  2.       SUBROUTINE SEPARM(MRIGID,RI1,RI2,NOUNIL,BDBLX,NELIM,IF)
  3. C----------------------------------------------------------------------
  4. C  Distinguer dans le MRIGID les rigidites elementaires pouvant etre 
  5. C  eliminees et celles devant etre conservees.
  6. C
  7. C Entrees :
  8. C ---------
  9. C  mrigid : pointeur sur la rigidite totale
  10. C  nelim  : ne rien eliminer lorsque vaut 0
  11. C  bdblx  : est vrai si les mult de Lagrange doivent etre dedoubles
  12. C     if  : numero de la passe d'elimination si elimination recursive
  13. C
  14. C Sorties :
  15. C ---------
  16. C     ri1 : pointeur sur les rigidites elementaires a eliminer, i.e. 
  17. C           pouvant etre traitee comme des dependances
  18. C     ri2 : pointeur sur les rigidites elementaires a conserver
  19. C
  20. C  juillet 2003 passage aux simples multiplicateurs de Lagrange, mais
  21. C  separm dualise ceux qu'il garde
  22. C----------------------------------------------------------------------
  23.       IMPLICIT INTEGER(I-N)
  24.       IMPLICIT REAL*8 (A-H,O-Z)
  25. C
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCREEL
  29. *-INC CCHAMP
  30. -INC SMRIGID
  31. -INC SMCOORD
  32. -INC SMELEME
  33. C
  34.       INTEGER OOOVAL
  35.       LOGICAL bdblx
  36. C
  37.       segment trav1
  38.        character*(lochpo) compp(lcomp1)
  39.       endsegment
  40.       segment trav2
  41.        integer ielim(nbnoc,nbincp),iautr(nbnoc,nbincp)
  42.        integer icomb(nbnoc,nbincp),ideja(nbnoc,nbincp)
  43.       endsegment
  44. C     Remplir posinc pour accelerer les tests sur les composantes
  45.       segment trav3
  46.        integer posinc(nligrp)
  47.       endsegment
  48.       segment trav4
  49.        integer itrv1(nbprl)
  50.        integer itrv2(nbprl)
  51.        real*8 dtrv1(nbprl)
  52.        real*8 dtrv2(nbprl)
  53.       endsegment
  54.       segment trav5
  55.        integer idata(2,nligrp)
  56.       endsegment
  57.       segment icpr(nbpts)
  58.       segment itemp(ntemp)
  59.       character*(lochpo) co1
  60.       character*(8) typmat
  61.       integer ri1p,ri1s,ri1f,ri2f
  62. C
  63. **     write(ioimp,*) ' entree de separm '
  64. **     call ooodmp(0)
  65. **     nbav=oooval(2,1)
  66. *      write (ioimp,*) ' dans separm nouinl ',nounil
  67. CC     write (ioimp,*) ' matrice mrigid'
  68. CC     call prrigi(mrigid,0)
  69. CC     segact,mrigid
  70.       nbdep=0
  71.       nbpiv=0
  73.       nbpvt=0
  74.       iunil=0
  75.       typmat='TEMPORAI'
  76. C     -----------------------------------------------------------------
  77. C     Distinguer ce qui peut etre elimine (ri4) et le reste (ri3)
  78. C     -----------------------------------------------------------------
  79.       if (nelim.ne.0) then
  80.         nbno=nbpts
  81.         segini icpr
  82.       endif
  83.       nbnoc=0
  84.       segact mrigid
  85.       nbrig=irigel(/2)
  86.       ntemp=nbrig
  87.       segini,itemp
  88.       nrige3=0
  89.       do 760 irig=1,nbrig
  90.         if (irigel(6,irig).ne.0) iunil=1
  91.         meleme=irigel(1,irig)
  92.         segact meleme
  93.         if ((itypel.ne.22.and.itypel.ne.28).or.
  94.      &      (irigel(7,irig).ne.0).or.(nelim.eq.0)) then
  95.           itemp(irig)=1
  96.           nrige3=nrige3+1
  97.         endif
  98. C
  99.         if (nelim.eq.0) goto 760
  100. C
  101. **      write(ioimp,*) '1 itypel dans separm',itypel
  102.         do il=1,num(/2)
  103.           do ip=1,num(/1)
  104.             ipt=num(ip,il)
  105.             if (icpr(ipt).eq.0) then
  106.               nbnoc=nbnoc+1
  107.               icpr(ipt)=nbnoc
  108.             endif
  109.           enddo
  110.         enddo
  111.  760  continue
  112. C
  113.       nrigel=nrige3
  114.       segini,ri3
  115.       ri3.iforig=iforig
  116.       ri3.mtymat=mtymat
  117.       nrige4=nbrig-nrige3
  118.       nrigel=nrige4
  119.       segini,ri4
  120.       ri4.iforig=iforig
  121.       ri4.mtymat=typmat
  122.       i3=0
  123.       i4=0
  124.       do 555 jj=1,nbrig
  125.         if (itemp(jj).eq.1) then
  126.           i3=i3+1
  127.           ii=i3
  128.           ri5=ri3
  129.         else
  130.           i4=i4+1
  131.           ii=i4
  132.           ri5=ri4
  133.         endif
  134.         ri5.coerig(ii)=coerig(jj)
  135.         do kk=1,8
  136.           ri5.irigel(kk,ii)=irigel(kk,jj)
  137.         enddo
  138.  555  continue
  139.       segsup,itemp
  140.       if (nelim.eq.0) then
  141.         ri1=ri4
  142.         ri2=ri3
  143.         GOTO 2002
  144.       endif
  145. C     -----------------------------------------------------------------
  146. C     Recuperer les noms de composantes (segment trav1)
  147. C     -----------------------------------------------------------------
  148.       lcomp1 = 50
  149.       segini trav1
  150.       nbincp=0
  151.       do 10 irig=1,nrige4
  152.         descr=ri4.irigel(3,irig)
  153.         segact descr
  154.         do 40 nligrp=2,lisinc(/2)
  155. C         print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
  156.           do 50 inc=1,nbincp
  157.             if (compp(inc).eq.lisinc(nligrp)) goto 40
  158.   50      continue
  159.           nbincp=nbincp+1
  160.           if (nbincp.GT.lcomp1) then
  161.             lcomp1 = lcomp1 + 50
  162.             segadj trav1
  163.           endif
  164.           compp(nbincp)=lisinc(nligrp)
  165.   40    continue
  166.   10  continue
  167. C     -----------------------------------------------------------------
  168. C     Quelles inconnues peuvent etre eliminees?
  169. C     -----------------------------------------------------------------
  170. C     Interdit d'eliminer les inconnues apparaissant dans des conditions 
  171. C     unilaterales lorsque nounil = 0.
  172.       segini,trav2
  173.       do 60 irig=1,nrige4
  174.         meleme=ri4.irigel(1,irig)
  175.         if(itypel.ne.28) then
  176.           if (ri4.irigel(6,irig).eq.0.or.nounil.eq.1) goto 60
  177.         endif
  178. C
  179. C       Super-elements : si trop gros ce n'est pas avantageux 
  180. C       d'eliminer. Dans ce cas, idnetifier ses noeuds et ddl pour ne 
  181. C       pas les eliminer.
  182.         if (itypel.ne.22.and.
  183.      >     (itypel.ne.28.or.num(/1).lt.255.or.if.lt.3)) goto 60
  184. C
  185.         descr=ri4.irigel(3,irig)
  186.         nld=2
  187.         if(itypel.eq.28) nld=1
  188. C
  189.         do 70 nligrp=nld,lisinc(/2)
  190.           do 80 inc=1,nbincp
  191. C            print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
  192. C     &             ,'  compp(',inc,')=',compp(inc)
  193.             if (compp(inc).eq.lisinc(nligrp)) goto 90
  194.   80      continue
  195. ***       write(ioimp,*) '- ',lisinc(nligrp),' non trouve dans trav1'
  196.           goto 70
  197.   90      continue
  198. C
  199. C         La composante lisinc(nligrp) ne doit pas etre eliminee pour
  200. C         les noeuds num(ip,j)
  201.           ip=noelep(nligrp)
  202.           do 100 j=1,num(/2)
  203. C            IF(num(ip,j).GT.ielim(/1).OR.inc.GT.ielim(/2).OR.
  204. C     *         num(ip,j).GT.iautr(/1).OR.inc.GT.iautr(/2))THEN
  205. C               print *,'BUG DANS SEPARM : ',
  206. C     &                 num(ip,j),'>',ielim(/1),iautr(/1),
  207. C     &                 inc,'>',ielim(/2),iautr(/2)
  208. C            ENDIF
  209.             ielim(icpr(num(ip,j)),inc)=1
  210.             iautr(icpr(num(ip,j)),inc)=1
  211.  100      continue
  212.   70    continue
  213.   60  continue
  214. C     -----------------------------------------------------------------
  215. C     Nombre de conditions associees a chaque ddl de chaque noeud
  216. C     -----------------------------------------------------------------
  217.       do 700 irig=1,nrige4
  218.         if (ri4.irigel(6,irig).ne.0.and.nounil.eq.0) goto 700
  219.         meleme=ri4.irigel(1,irig)
  220.         if (itypel.ne.22) goto 700
  221.         descr=ri4.irigel(3,irig)
  222. C
  223.         nligrp=lisinc(/2)
  224.         segini trav3
  225.         do 730 inc=2,nligrp
  226.           do 720 incp=1,nbincp
  227.             if (lisinc(inc).eq.compp(incp)) then
  228.               posinc(inc)=incp
  229.               goto 730
  230.             endif
  231.   720     continue
  232.           call erreur(5)
  233.   730   continue
  234. C
  235.         do 750 ince=2,lisinc(/2)
  236.           incp=posinc(ince)
  237.           do 740 j=1,num(/2)
  238.             ip=num(noelep(ince),j)
  239.             icomb(icpr(ip),incp)=icomb(icpr(ip),incp)+1
  240.   740     continue
  241.   750   continue
  242.         segsup trav3
  243.   700 continue
  244.       ipass=1
  245. C
  246. C     -----------------------------------------------------------------
  247.  2000 CONTINUE
  248.       nrigel=0
  249.       segini,ri1
  250.       ri1.iforig=iforig
  251.       ri1.mtymat=typmat
  252.       segini,ri2=ri4
  253.       ri2.mtymat=mtymat
  254. C
  255. C     Trier pour attaquer en premier les relations portant sur le moins
  256. C     d'inconnues
  257.       nbprl=ri4.irigel(/2)
  258.       segini trav4
  259.       do 765 irig=1,nbprl        
  260.         descr=ri4.irigel(3,irig)
  261.         dtrv1(irig)= lisinc(/2)+(irig/(nbprl+1.))
  262.         if(irigel(6,irig).ne.0) dtrv1(irig)=dtrv1(irig)+1000000
  263.         if(irigel(6,irig).eq.2) dtrv1(irig)=dtrv1(irig)+1000000
  264. **      if(irigel(6,irig).eq.0) dtrv1(irig)=dtrv1(irig)+1000000
  265. **      itrv1(irig)=nbprl-irig+1
  266.         itrv1(irig)=      irig   
  267. 765   continue
  268.       call triflo(dtrv1,dtrv2,itrv1,itrv2,nbprl)
  269. C
  270.       do 200 iri=1,nbprl
  271.         irig=itrv1(iri) 
  272. **      irig=iri   
  273. **      if (ipass.eq.1) irig=itrv1(iri)
  274.         if (ri4.irigel(6,irig).ne.0.and.nounil.eq.0) goto 190
  275.         meleme=ri4.irigel(1,irig)
  276.         if (itypel.ne.22) goto 190
  277.         Xmatri=ri4.irigel(4,irig)
  278.         segact Xmatri
  279.         if (abs(re(1,1,1)).gt.1d-30) goto 190
  280.         descr=ri4.irigel(3,irig)
  281. C
  282.         nligrp=lisinc(/2)
  283.         segini,trav3,trav5
  284.         do 230 inc=2,nligrp
  285.          do 220 incp=1,nbincp
  286.            if (lisinc(inc).eq.compp(incp)) then
  287.              posinc(inc)=incp
  288.              goto 230
  289.            endif
  290.  220     continue
  291.         call erreur(5)
  292.  230    continue
  293. C
  294.         segini,ipt8=meleme
  295.         nbdepe=0
  296.         do 210 j=1,num(/2)
  297. C
  298. C         Cette matrice elementaire contient elle des inconnues
  299.           if (noelep(/1).le.1) then
  300.             ipt8.num(1,j)=0
  301.             goto 210
  302.           endif
  303. C
  304. C         La matrice est elle augmentee
  305.           if (abs(re(1,1,j)).gt.1d-5) then
  306.             ipt8.num(1,j)=0
  307.             goto 210
  308.           endif
  309. C
  310. C         Cette matrice elementaire contient-elle une inco. eliminee
  311.           remax=abs(re(1,1,j))
  312.           do 240 inc=2,nligrp
  313.             incpp=posinc(inc)
  314.             ipp=ipt8.num(noelep(inc),j)
  315.             if (ielim(icpr(ipp),incpp).eq.1) then
  316.               ipt8.num(1,j)=0
  318.               goto 210
  319.             endif
  320.             remax=max(remax,abs(re(1,inc,j)))
  321.  240      continue
  322. C
  323. C         Choix du pivot
  324.           incf=0
  325.           remaz=remax*0.9
  326.           do 250 inc=2,nligrp
  327.             incpp=posinc(inc)
  328.             ipp=ipt8.num(noelep(inc),j)
  329.             if (iautr(icpr(ipp),incpp).eq.1) goto 250
  330.             if (abs(re(1,inc,j)).lt.remaz) goto 250
  331. *           if (icomb(icpr(ipp),incpp).eq.1) then
  332.               incf=inc
  333.               goto 260
  334. *           endif
  335.  250      continue
  336.  260      continue
  337.           ince=incf
  338.           if (incf.eq.0) ince=2
  339. C
  340. C         Traitement de l'inconnue incp du noeud ip
  341.           incp=posinc(ince)
  342.           ip=ipt8.num(noelep(ince),j)
  343. C
  344. C         Le pivot est il correct
  345.           remax=remax*1d-2
  346.           if (abs(re(1,ince,j)).le.remax) then
  347.             ipt8.num(1,j)=0
  348.             nbpiv=nbpiv+1
  349.             goto 210
  350.           endif
  351. C
  352. C         Cette inconnue apparait-elle dans d'autres CL
  353.           if (ipass.eq.1.and.icomb(icpr(ip),incp).ne.1) then
  354.             ipt8.num(1,j)=0
  356.             goto 210
  357.           endif
  358. C
  359. C         Cette inconnue apparait-elle dans une dependance
  360.           if (iautr(icpr(ip),incp).eq.1) then
  361.             ipt8.num(1,j)=0
  363.             goto 210
  364.           endif
  365. C
  366. C         Cette inconnue apparait-elle deux fois dans la relation
  367.           ideux=0
  368.           do inc=2,nligrp
  369.          if (ideja(icpr(ipt8.num(noelep(inc),j)),posinc(inc)).eq.1) then
  370.               ideux=inc
  371.             else
  372.               ideja(icpr(ipt8.num(noelep(inc),j)),posinc(inc))=1
  373.             endif
  374.           enddo
  375.           do inc=2,nligrp
  376.             ideja(icpr(ipt8.num(noelep(inc),j)),posinc(inc))=0
  377.           enddo
  378.           if (ideux.ne.0) then
  379.             moterr(1:4)=lisinc(ideux)
  380.             interr(1)=ipt8.num(noelep(ideux),j)
  381.             call erreur(-361)
  382.             ipt8.num(1,j)=0
  383.             goto 210
  384.           endif
  385. C
  386. C         Nouvelle dependance
  387.           nbdepe=nbdepe+1
  388.           nbdep=nbdep+1
  389. C         write (ioimp,*) 'Elimination noeud,inco,posi',ip,lisinc(ince),ince
  390. C         Reperer l'inconnue eliminee et le noeud 
  391.           ipt8.num(1,j)=ince
  392.           idata(1,ince)=idata(1,ince)+1
  393.           ielim(icpr(ip),incp)=1
  394.           do 280 inc=2,nligrp
  395.             iautr(icpr(ipt8.num(noelep(inc),j)),posinc(inc))=1
  396.  280      continue
  397. C
  398.  210    continue
  399.         segsup trav3
  400. C
  401. C       Creation de ri1 et ri2 pour scinder la sous-matrice irig
  402. C       ---------------------------------------------------------------
  403. C       Dimensions communes a RI1 et RI2
  404.         nbnn=num(/1)
  405.         nbsous=0
  406.         nbref=0
  407.         nligrd=re(/1)
  408.         nligrp=re(/2)
  409. C
  410. C       Creation de RI2 : rigidites elementaires a conserver
  411.         nbelem=num(/2)-nbdepe
  412.         if (nbelem.gt.0) then
  413.           segini,ipt2
  414.           ipt2.itypel=itypel
  415.           nelrig=nbelem
  416.           segini,xmatr2
  417.           xmatr2.symre=symre
  418.           ri2.irigel(1,irig)=ipt2
  419.           ri2.irigel(3,irig)=descr
  420.           ri2.irigel(4,irig)=xmatr2
  421.           do iii=5,8
  422.             ri2.irigel(iii,irig)=ri4.irigel(iii,irig)
  423.           ENDDO
  424.           ri2.coerig(irig)=ri4.coerig(irig)
  425.         else
  426.           ri2.irigel(1,irig)=0
  427.           xmatr2=0
  428.         endif
  429. C
  430. C       Creation de RI1 : rigidites elementaires a eliminer. 
  431. C       L'inconnue a eliminer doit etre en 2eme position (apres LX)
  432. C       Si l'inconnue a eliminer n'est pas celle apparaissant en   
  433. C       deuxieme position il faut pivoter la rigidite elementaire 
  434. C       et le descripteur associe
  435.         ri1p=ri1
  436.         nrigel=0
  437.         do kkk=1,nligrp
  438.           if (idata(1,kkk).gt.0) nrigel=nrigel+1
  439.         enddo
  440.         nbpvt=nbpvt+nrigel
  441.         segini,ri1
  442.         ri1.iforig=iforig
  443.         ri1.mtymat=typmat
  444.         icpt=0
  445.         do 25 jjj=1,nligrp
  446.           nelrig=idata(1,jjj)
  447. C         Nombre d'elements elimines en choisissant la jjjeme inconnue
  448.           if (nelrig.eq.0) goto 25
  449.           nbelem=nelrig
  450.           segini,ipt1
  451.           ipt1.itypel=itypel
  452.           segini,xmatr1
  453.           xmatr1.symre=symre
  454.           des1=descr
  455.           if (jjj.ne.2) then
  456.             segini,des1=descr
  457.             co1=des1.lisinc(2)
  458.             des1.lisinc(2)=des1.lisinc(jjj)
  459.             des1.lisinc(jjj)=co1
  460.             co1=des1.lisdua(2)
  461.             des1.lisdua(2)=des1.lisdua(jjj)
  462.             des1.lisdua(jjj)=co1
  463.             noe=des1.noelep(2)
  464.             des1.noelep(2)=des1.noelep(jjj)
  465.             des1.noelep(jjj)=noe
  466.             noe=des1.noeled(2)
  467.             des1.noeled(2)=des1.noeled(jjj)
  468.             des1.noeled(jjj)=noe
  469.           endif
  470. C
  471.           icpt=icpt+1
  472.           ri1.irigel(1,icpt)=ipt1
  473.           ri1.irigel(3,icpt)=des1
  474.           ri1.irigel(4,icpt)=xmatr1
  475.           do iii=5,7
  476.             ri1.irigel(iii,icpt)=ri4.irigel(iii,irig)
  477.           enddo
  478.           ri1.irigel(8,icpt)=1
  479.           ri1.coerig(icpt)=ri4.coerig(irig)
  480. C
  481. C         Compteur des elements et identifiant de la rigidite elem
  482.           idata(1,jjj)=1
  483.           idata(2,jjj)=icpt
  484.  25     continue
  485. C
  486. C       Remplir ri1 et ri2 en creant les xmatri/descr/maillages
  487. C       ---------------------------------------------------------------
  488.         idec=0
  489.         do 300 j=1,num(/2)
  490.           if (ipt8.num(1,j).eq.0) then
  491.             idec=idec+1
  492.             do 310 i=1,num(/1)
  493.              ipt2.num(i,idec)=num(i,j)
  494.  310        continue
  495.             do io=1,nligrp
  496.               do iu=1,nligrd
  497.                 xmatr2.re(iu,io,idec)=re(iu,io,j)
  498.               enddo
  499.             enddo
  500.           else
  501.             ince=ipt8.num(1,j)
  502.             iriel=idata(2,ince)
  503.             ipt1=ri1.irigel(1,iriel)
  504.             xmatr1=ri1.irigel(4,iriel)
  505.             ielt=idata(1,ince)
  506.             idata(1,ince)=idata(1,ince)+1
  507.             do 320 i=1,num(/1)
  508.               ipt1.num(i,ielt)=num(i,j)
  509.  320        continue
  510.             do io=1,nligrp
  511.               do iu=1,nligrd
  512.                 xmatr1.re(iu,io,ielt)=re(iu,io,j)
  513.               enddo
  514.             enddo
  515. C
  516. C           Pivoter les lignes/colonnes ince et 2
  517.             if (ince.ne.2) then
  518.               do 1130 il=1,nligrd
  519.                 ret=xmatr1.re(il,2,ielt)
  520.                 xmatr1.re(il,2,ielt)=xmatr1.re(il,ince,ielt)
  521.                 xmatr1.re(il,ince,ielt)=ret
  522.  1130         continue
  523.               do 1160 il=1,nligrp
  524.                 ret=xmatr1.re(2,il,ielt)
  525.                 xmatr1.re(2,il,ielt)=xmatr1.re(ince,il,ielt)
  526.                 xmatr1.re(ince,il,ielt)=ret
  527.  1160         continue
  528.             endif
  529. C
  530.           endif
  531.  300    continue
  532. C       ---------------------------------------------------------------
  533. C
  534.         call fusrig(ri1p,ri1,ri1f)
  535.         ri1=ri1f
  536.         if (xmatr2.ne.0) segdes,xmatr2
  537.         segsup,trav5
  538.         goto 200
  539. C
  540.  190    continue
  541. C
  542. C       Sous-matrice irig a conserver integralement
  543.         ri2.irigel(1,irig)=ri4.irigel(1,irig)
  544.         ri2.irigel(3,irig)=ri4.irigel(3,irig)
  545.         ri2.irigel(4,irig)=ri4.irigel(4,irig)
  546.         do iii=5,7
  547.           ri2.irigel(iii,irig)=ri4.irigel(iii,irig)
  548.         enddo
  549.         ri2.irigel(8,irig)=0
  550.         ri2.coerig(irig)=ri4.coerig(irig)
  551. C
  552.  200  continue
  553. C     -----------------------------------------------------------------
  554.       segsup trav4
  555. C
  556. C     Compression de ri2
  557.       idec=0
  558.       do 600 irig=1,nbprl
  559.         meleme=ri2.irigel(1,irig)
  560.         if (meleme.eq.0) then
  561.           idec=idec+1
  562.         else
  563.           do 610 ir=1,8
  564.            ri2.irigel(ir,irig-idec)=ri2.irigel(ir,irig)
  565.  610      continue
  566.           ri2.coerig(irig-idec)=ri2.coerig(irig)
  567.         endif
  568.  600  continue
  569.       nrigel=nbprl-idec
  570. C     write (ioimp,*) ' dimension de ri2 ',nrigel
  571.       segadj,ri2
  572. C
  573. C     On va voir si on ne peut pas faire pivoter ri2
  574.       if (ipass.eq.2) goto 2001
  575.       ri1s=ri1
  576.       ri4=ri2
  577.       nbpiv=0
  579.       ipass=ipass+1
  580.       goto 2000
  581. C
  582.  2001 CONTINUE
  583. C     -----------------------------------------------------------------
  584. C
  585. C     RI1 : rigidites elementaires a eliminer
  586.       call fusrig(ri1s,ri1,ri1f)
  587.       ri1=ri1f
  588. C     RI2 : rigidites elementaires a conserver
  589.       call fusrig(ri3,ri2,ri2f)
  590.       ri2=ri2f
  591. C
  592.       segsup trav1,trav2,icpr
  593. C
  594.       if (iimpi.ne.0) then
  595.       write(ioimp,*)'nombre de relations eliminees',nbdep
  596.       write(ioimp,*)'nombre de relations gardees a cause du pivot',nbpiv
  597.       write(ioimp,*)'nombre de relations gardees car non independantes',
  599.       write(ioimp,*)'nombre de paquets pivotes',nbpvt
  600.       endif
  601. C
  602.  2002 CONTINUE
  603. C
  604. C     -----------------------------------------------------------------
  605. C     Dualisation des multiplicateurs de Lagrange
  606. C     -----------------------------------------------------------------
  607. C     si on a des conditions unilaterales, on ne dualise pas, ce sera fait
  608. C     dans le resou de unilater
  609.       if ((iunil.eq.0.or.nounil.eq.1).and.bdblx) call dbblx(ri2)
  610. C
  611. CC    write (ioimp,*) ' matrice ri1 '
  612. CC    call prrigi(ri1,0)
  613. CC    write (ioimp,*) ' matrice ri2 '
  614. CC    call prrigi(ri2,0)
  615. CC    segact,ri1,ri2
  616. CC    write(ioimp,*) ' sortie de separm '
  617. **    call ooodmp(0)
  618. **    nbap=oooval(2,1)
  619. **    write(ioimp,*) 'nb segmts dans separm avant apres ',nbav,nbap
  620.       END
  621.  
  622.  
  623.  

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