C SEPARM    SOURCE    FANDEUR   13/02/07    21:15:14     7694           
      subroutine separm(mrigid,ri1,ri2,nounil,lagdua)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC SMRIGID
-INC SMCOORD
-INC SMELEME
-INC CCOPTIO
*  extrait de mrigid dans ri1 la partie des relations pouvant etre 
*  traitees comme des dependances (elimination des inconnues)
*  il faut que la premiere inconnue d'une dependance n'apparaisse pas 
*  dans les autres dependances
*  juillet 2003 passage aux simples multiplicateurs de Lagrange, mais 
*  separm dualise ceux qu'il garde
*-INC CCHAMP
      segment trav1
       character*4 compp(lcomp1)
C*       character*4 compp(nbincp)
      endsegment
      segment trav2
       integer ielim(nbno,nbincp),iautr(nbno,nbincp),icomb(nbno,nbincp)
      endsegment
      segment trav3
       integer posinc(nligrp)
      endsegment
      character*4 co1
      integer ri1s,ri2s,ri1f

*      write (6,*) ' dans separm nouinl ',nounil
      nbdep=0
      nbpiv=0
      nbinc=0
      nbpvt=0
      nbno=xcoor(/1)/(idim+1)
      segact mrigid
*  remplissage du segment trav1 (liste des noms de composantes)
      lcomp1 = 50
      segini trav1
      nbincp=0
      do 10 irig=1,irigel(/2)
       if (irigel(6,irig).ne.0.and.irigel(6,irig).ne.2
     &     .and.nounil.eq.0) goto 10
       meleme=irigel(1,irig)
       segact meleme
       if (itypel.ne.22) goto 20
       descr=irigel(3,irig)
       segact  descr
         do 40 nligrp=2,lisinc(/2)
C           print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
           do 50 inc=1,nbincp
             if (compp(inc).eq.lisinc(nligrp)) goto 40
  50       continue
           nbincp=nbincp+1
           if (nbincp.GT.lcomp1) then
             lcomp1 = lcomp1 + 50
             segadj trav1
           endif
           compp(nbincp)=lisinc(nligrp)
  40     continue
       segdes,descr
  20   continue
       segdes,meleme
  10  continue
c*      lcomp1 = nbincp
c*      segadj trav1
*  on construit ri1 et ri2 en notant dans le segment trav2 ce qu'on fait
      segini trav2
*  on interdit d'eliminer les inconnues apparaissant dans le frottement
      do 60 irig=1,irigel(/2)
       if (irigel(6,irig).eq.0.or.nounil.eq.1) goto 60
       meleme=irigel(1,irig)
       segact meleme
       if (itypel.ne.22) goto 60
       descr=irigel(3,irig)
       segact  descr
         do 70 nligrp=2,lisinc(/2)
           do 80 inc=1,nbincp
C            print *,'lisinc(',nligrp,'/',lisinc(/2),')=',lisinc(nligrp)
C     &             ,'  compp(',inc,')=',compp(inc)
             if (compp(inc).eq.lisinc(nligrp)) goto 90
  80       continue
c*      write(ioimp,*) '- ',lisinc(nligrp),' non trouve dans trav1'
           goto 70
  90       continue
           ip=noelep(nligrp)
           do 100 j=1,num(/2)
C            IF(num(ip,j).GT.ielim(/1).OR.inc.GT.ielim(/2).OR.
C     *         num(ip,j).GT.iautr(/1).OR.inc.GT.iautr(/2))THEN
C               print *,'BUG DANS SEPARM : ',
C     &                 num(ip,j),'>',ielim(/1),iautr(/1),
C     &                 inc,'>',ielim(/2),iautr(/2)
C            ENDIF
            ielim(num(ip,j),inc)=1
            iautr(num(ip,j),inc)=1
 100       continue
  70     continue
       segdes descr,meleme
  60   continue
       do 700 irig=1,irigel(/2)
        if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 690
        meleme=irigel(1,irig)
        descr=irigel(3,irig)
        xmatri=irigel(4,irig)
        segact meleme
        if (itypel.ne.22) goto 690
        segact descr
*  remplir posinc pour accelerer les tests sur les composantes
        nligrp=lisinc(/2)
        segini trav3
        do 730 inc=2,nligrp
         do 720 incp=1,nbincp
          if (lisinc(inc).eq.compp(incp)) then
           posinc(inc)=incp
           goto 730
          endif
  720    continue
         call erreur(5)
  730   continue
        do 750 ince=2,lisinc(/2)
         incp=posinc(ince)
         do 740 j=1,num(/2)
          ip=num(noelep(ince),j)
          icomb(ip,incp)=icomb(ip,incp)+1
  740    continue
  750   continue
        segsup trav3
  690   continue
  700  continue
      ipass=1
 2000 continue
      segini,ri1=mrigid
      segini,ri2=mrigid
      do 200 irig=1,irigel(/2)
        if (irigel(6,irig).ne.0.and.nounil.eq.0) goto 190
        meleme=irigel(1,irig)
        descr=irigel(3,irig)
        Xmatri=irigel(4,irig)
        segact meleme
        if (itypel.ne.22) goto 190
        segact descr,Xmatri
*  remplir posinc pour accelerer les tests sur les composantes
        nligrp=lisinc(/2)
        segini trav3
        do 230 inc=2,nligrp
         do 220 incp=1,nbincp
           if (lisinc(inc).eq.compp(incp)) then
             posinc(inc)=incp
             goto 230
           endif
 220     continue
        call erreur(5)
 230    continue
        nbdepe=0
*  pour arriver à eliminer, il est interessant de faire pivoter les matrices
*  on va le faire en se basant sur le milieu du paquet
        iel=num(/2)/2+1
        incf=0
        do 1000 ince=2,lisinc(/2)
          ipp=num(noelep(ince),iel)
          incpp=posinc(ince)
*  cette inconnue apparait elle dejà dans une autre relation
         if (iautr(ipp,incpp).eq.1) goto 1000
*  cette matrice touche t elle une autre relation
         do 1240 inc=2,lisinc(/2)
         if (inc.eq.ince) goto 1240
          incp=posinc(inc)
          ip=num(noelep(inc), iel)
          if (ielim(ip,incp).eq.1) goto 1000
1240     continue
*  le pivot est il correct
*         xmatri=imattt(iel)
*         segact xmatri
         remax=abs(re(ince,1,iel))
         do 1270 ilig=2,re(/1)
           remax=max(remax,abs(re(ilig,1,iel)))
1270     continue
         remax=remax*1e-2
         if (abs(re(ince,1,iel)).lt.remax) goto 1000
*         segdes xmatri
*  ok on pivote si il y a lieu
*        write (6,*) ' ipp incpp icomb ',ipp,incpp,icomb(ipp,incpp)
         if (icomb(ipp,incpp).eq.1) goto 1010
         if (incf.eq.0) incf=ince  
         incf=ince
         goto 1000
1000    continue
        ince=2
        if (incf.ne.0) ince=incf
1010    continue
        if (ince.ne.2) then
*  on pivote 2 et ince
*       write (6,*) ' on pivote avec la ligne ',ince
        nbpvt=nbpvt+1
        incp=posinc(ince)
        posinc(ince)=posinc(2)
        posinc(2)=incp
        segini,des1=descr
        segdes descr
        descr=des1
        ri1.irigel(3,irig)=descr
        ri2.irigel(3,irig)=descr
*       write (6,*) noelep(2),noelep(ince),lisinc(2),lisinc(ince)
        co1=lisinc(2)
        lisinc(2)=lisinc(ince)
        lisinc(ince)=co1
        co1=lisdua(2)
        lisdua(2)=lisdua(ince)
        lisdua(ince)=co1
        noe=noelep(2)
        noelep(2)=noelep(ince)
        noelep(ince)=noe
        noe=noeled(2)
        noeled(2)=noeled(ince)
        noeled(ince)=noe
        segini,Xmatr1=Xmatri
        segdes Xmatri
        Xmatri=Xmatr1
        ri1.irigel(4,irig)=Xmatri
        ri2.irigel(4,irig)=Xmatri
        do 1100 im=1,re(/3)
*        xmatri=imattt(im)
*         segini,xmatr1=xmatri
*         segdes xmatri
*         xmatri=xmatr1
*         imattt(im)=xmatri
         do 1130 il=1,re(/1)
          ret=re(il,2,im)
          re(il,2,im)=re(il,ince,im)
          re(il,ince,im)=ret
 1130    continue
         do 1160 il=1,re(/2)
          ret=re(2,il,im)
          re(2,il,im)=re(ince,il,im)
          re(ince,il,im)=ret
 1160    continue
*         segdes xmatri
 1100   continue
        endif
*  peut etre une dependance
        segini,ipt1=meleme
        segini,xmatr1=xmatri
        do 210 j=1,num(/2)
          ip=ipt1.num(noelep(2),j)
*  cette matrice contient t elle des inconnues
          if (noelep(/1).le.1) then
           ipt1.num(1,j)=0
           goto 210
          endif
          incp=posinc(2)
         if (ipass.eq.1.and.icomb(ipp,incpp).ne.1) then
            ipt1.num(1,j)=0
            nbinc=nbinc+1
            goto 210
         endif
         if (iautr(ip,incp).eq.1) then
           ipt1.num(1,j)=0
           nbinc=nbinc+1
           goto 210
         endif
*  cette matrice touche t elle une autre relation
         do 240 inc=2,lisinc(/2)
          incp=posinc(inc)
          ip=ipt1.num(noelep(inc),j)
          if (ielim(ip,incp).eq.1) then
            ipt1.num(1,j)=0
            nbinc=nbinc+1
            goto 210
          endif
 240      continue
*  le pivot est il correct
*          xmatri=imattt(j)
*          segact xmatri
          remax=abs(re(2,1,j))
          do 270 ilig=2,re(/1)
            remax=max(remax,abs(re(ilig,1,j)))
 270      continue
          remax=remax*1e-2
          if (abs(re(2,1,j)).lt.remax) then
            ipt1.num(1,j)=0
            nbpiv=nbpiv+1
            goto 210
          endif
*          segdes xmatri
*  ok on la prend en dependance
          nbdepe=nbdepe+1
          nbdep=nbdep+1
*         write (6,*) ' inconnu eliminee ',ipt1.num(noelep(2),j),
*    >      lisinc(2)
          ielim(ipt1.num(noelep(2),j),posinc(2))=1
          do 280 inc=2,lisinc(/2)
            iautr(ipt1.num(noelep(inc),j),posinc(inc))=1
 280      continue
 210    continue
        segsup trav3
* il faut maintenant tasser ipt1 et creer ipt2
        nbnn=num(/1)
        nbelem=num(/2)-nbdepe
*       write(6,*) ' nombre d elements crees pour ipt2 ',nbelem
        nbsous=0
        nbref=0
        segini ipt2
        ipt2.itypel=itypel
        nelrig=nbelem
        nligrd=re(/1)
        nligrp=re(/2)
        segini xmatr2
        idec=0
        do 300 j=1,num(/2)
         if (ipt1.num(1,j).eq.0) then
           idec=idec+1
           do 310 i=1,num(/1)
            ipt2.num(i,idec)=num(i,j)
 310       continue
            do io=1,nligrp
                 do iu=1,nligrd
                   xmatr2.re(iu,io,idec)=re(iu,io,j)
                 enddo
            enddo
*           imatr2.imattt(idec)=imattt(j)
         else
           do 320 i=1,num(/1)
            ipt1.num(i,j-idec)=num(i,j)
 320       continue
            do io=1,nligrp
                 do iu=1,nligrd
                   xmatr1.re(iu,io,j-idec)=re(iu,io,j)
                 enddo
            enddo 
*           imatr1.imattt(j-idec)=imattt(j)
         endif
 300    continue
        nbelem=nbdepe
        segadj ipt1
        nelrig=nbelem
        segadj xmatr1
        ri1.irigel(1,irig)=ipt1
        ri2.irigel(1,irig)=ipt2
        ri1.irigel(4,irig)=xmatr1
        ri2.irigel(4,irig)=xmatr2
        segdes ipt1,ipt2,meleme,xmatr1,xmatr2
        goto 200
 190    continue
*  pas une dependance
        ri1.irigel(1,irig)=0
 200  continue
*  compression de ri1
      idec=0
      do 500 irig=1,irigel(/2)
        meleme=ri1.irigel(1,irig)
        if (meleme.eq.0) then
          idec=idec+1
        else
          segact meleme
          if (num(/2).eq.0) then
            idec=idec+1
          else
            do 510 ir=1,irigel(/1)
             ri1.irigel(ir,irig-idec)=ri1.irigel(ir,irig)
 510        continue
            ri1.coerig(irig-idec)=ri1.coerig(irig)
          endif
          segdes meleme
        endif
 500  continue
      nrige=max(irigel(/1),8)
      nrigel=irigel(/2)-idec
*     write (6,*) ' dimension de ri1 ',nrige,nrigel
      segadj ri1
      do irig=1,nrigel
        ri1.irigel(8,irig)=1
      enddo
*  compression de ri2
      idec=0
      do 600 irig=1,irigel(/2)
        meleme=ri2.irigel(1,irig)
        if (meleme.eq.0) then
          idec=idec+1
        else
          segact meleme
          if (num(/2).eq.0) then
            idec=idec+1
          else
            do 610 ir=1,irigel(/1)
             ri2.irigel(ir,irig-idec)=ri2.irigel(ir,irig)
 610        continue
            ri2.coerig(irig-idec)=ri2.coerig(irig)
          endif
          segdes meleme
        endif
 600  continue
      nrige=max(irigel(/1),8)
      nrigel=irigel(/2)-idec
*     write (6,*) ' dimension de ri2 ',nrige,nrigel
      segadj ri2
      do irig=1,nrigel
        ri2.irigel(8,irig)=0
      enddo
*  On va voir si on ne peut pas faire pivoter ri2
      if (ipass.eq.2) goto 2001
      ri1s=ri1
      ri2s=ri2
      mrigis=mrigid
      mrigid=ri2
      nbpiv=0
      nbinc=0
      ipass=2
      goto 2000
 2001 continue 
      mrigid=mrigis
      call fusrig(ri1s,ri1,ri1f)
      ri1=ri1f
      segsup trav1,trav2
      if (iimpi.ne.0) then
      write (6,*) ' nombre de relations éliminées ',nbdep
      write (6,*) ' nombre de relations gardées à cause du pivot ',
     >            nbpiv
      write (6,*) ' nombre de relations gardées car non indépendantes ',
     >            nbinc
      write (6,*) ' nombre de paquets pivotés ',
     >            nbpvt
      endif
*  dualisation des multiplicateurs de Lagrange dans ri2   
*  si on a des conditions unilatérales, on ne dualise pas, ce sera fait
*  dans le resou de unilater
       iunil=0
       segact ri2
       do ir=1,ri2.irigel(/2)
       if (ri2.irigel(6,ir).ne.0) iunil=1
       enddo
       if ((iunil.eq.0.or.nounil.eq.1).and.lagdua.ge.0)
     $      call dbblx(ri2,lagdua)
*      write (6,*) ' matrice ri1 '
*      call prrigi(ri1)
*      write (6,*) ' matrice ri2 '
*      call prrigi(ri2)

      return
      end