C RESOU     SOURCE    PV090527  25/10/15    21:15:03     12383          
          
      SUBROUTINE RESOU
C----------------------------------------------------------------------
C **** CET OPERATEUR SERT A RESOUDRE UN SYSTEME D EQUATIONS LINEAIRES
C ****   CHPOINT = RESOU RIGIDITE CHPOINT
C----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
-INC PPARAM
-INC CCOPTIO
-INC SMRIGID
-INC SMCOORD
-INC SMTEXTE
-INC SMTABLE
-INC SMCHPOI
-INC SMELEME
-INC SMLCHPO
-INC CCREEL
C
      PARAMETER(ZERO=0.D0)
      SEGMENT IDEMEM(0)
      segment ideme0(idemem(/1),30)
      segment ideme1(idemem(/1),30)
      segment idnote(0)
C
      CHARACTER*4  LISM(9)
      CHARACTER*72 CHARRE
      REAL*8 XVA
      LOGICAL ILOG,ILUG,casfimp,notok,bdblx
      DATA LISM/'NOID','NOUN','ENSE','GRAD','CHOL','STAB','ELIM',
     >'NOST','SOUC'/
      DATA ILOG/.FALSE./
C
      XVA=REAL(0.D0)
      IOB=0
C
*  sauvegarde norinc
*     write(6,*) 'norinc norind ',norinc,norind
      norins=norinc
      iverif=1
      ipt8=0
      iunil=0
*  le defaut est de faire une passe d'elimination
      nelim=30
*  experimentalement 2 passes est mieux
***   nelim=2
      IMTVID=0
      NOUNIL=0
      NOID=0
      NOEN=1
      IGRADJ = 0
      ICHSKI = 0
      INSYM  = 0
      KIKI=0
      KSYMET = 0
      IPSHPO = 0
      ISTAB=0
      ISOUCI = 0
      ifochs = -99
      NDEMEM=0
C----------------------------------------------------------------------
C     LECTURE DES ARGUMENTS DE L'OPERATEUR RESOU
C----------------------------------------------------------------------
C     LECTURE DES OPTIONS/MOTS-CLES
  5   CONTINUE
      CALL LIRMOT(LISM,9,KIKI,0)
      IF (KIKI.EQ.1) NOID=1
      IF (KIKI.EQ.2) NOUNIL=1
      IF (KIKI.EQ.3) NOEN=0
*     IF (KIKI.EQ.4) IGRADJ = 1
*     IF (KIKI.EQ.5) ICHSKI = 1
*     IF (KIKI.EQ.4.OR.KIKI.EQ.5) KSYMET = 1
      IF (KIKI.EQ.6) ISTAB=1
      IF (KIKI.EQ.7) then
       call lirent(nelim,1,iretou)
       nelim=min(30,max(0,nelim))
      endif
      IF (KIKI.EQ.8) ISTAB=0
      IF (KIKI.EQ.9) ISOUCI=1
      IF (KIKI.NE.0) GOTO 5
C
      if (noid.eq.1) iverif=0
      IF(NUCROU.EQ.0) THEN
        ICHSKI=1
      ELSEIF(NUCROU.EQ.1) THEN
        IGRADJ=1
        KSYMET=1
      ENDIF
*      WRITE(6,*) ' nucrou', nucrou
*      IF ( IGRADJ + ICHSKI .EQ. 0 ) ICHSKI = 1
C
C     LECTURE DE LA RIGIDITE
      CALL LIROBJ('RIGIDITE',IPOIRI,1,IRETOU)
      IF(IERR.NE.0) GO TO 5000
      IPRIGO=IPOIRI
C
C     LECTURE DE LA PRECISION
      PREC=REAL(xspeti)
      CALL LIRREE(PREC,0,IRETOU)
      IF(IERR.NE.0) GO TO 5000
C
C     REMPLISSAGE DU 2ND MEMBRE IDEMEM(**) A PARTIR DE
C     ... CHPOINT
      SEGINI IDEMEM
   1  CONTINUE
      CALL LIROBJ('CHPOINT ',ISECO,0,IRETOU)
      IF(IRETOU.NE.0) THEN
        IDEMEM(**)=ISECO
        NDEMEM=NDEMEM+1
        mchpoi=iseco
        segact mchpoi
        if (mchpoi.jattri(/1).ge.1) then
          notok=(mchpoi.jattri(1).ne.2)
        else
          notok=.true.
        endif
        if (notok) then
          INTERR=mchpoi
          CALL ERREUR(-387)
        endif
        GOTO 1
      ENDIF
C
C     ... LISTCHPO
      CALL LIROBJ('LISTCHPO',ISECO,0,IRETOU)
      IF(IRETOU.NE.0) THEN
        mlchpo=ISECO
        segact mlchpo
        NDEMEM = ichpoi(/1)
        do iu = 1,NDEMEM
          idemem(**) = ichpoi(iu)
*         write(6,*) ' extension idemem 2 ',idemem(/1)
        enddo
        segdes mlchpo
        n1 = ndemem
        segini mlchpo
        ipshpo = mlchpo
      ENDIF
      IF (IERR.NE.0) RETURN
C
C     ... TABLE DE SOUS-TYPE LIAISONS_STATIQUES
      CALL LIRTAB('LIAISONS_STATIQUES',ITBAS,0,IRET)
      IF (ITBAS.EQ.0) GOTO 90
C
      IF (IIMPI.EQ.333) THEN
        WRITE(IOIMP,*) 'on a lu la table des conditions aux limites'
      ENDIF
      mtab1 = itbas
      segact mtab1
      ima = mtab1.mlotab - 1
      segini idnote
      im = 0
      segdes mtab1
 80   CONTINUE
      im = im + 1
      itmod = 0
      ichp0 = 0
      if (im.gt.ima) then
        IF (NDEMEM.GT.0) GOTO 90
*       pas de champs de force
        CALL ERREUR(1)
        return
      endif
      CALL ACCTAB(ITBAS,'ENTIER',IM,0.d0,' ',.true.,IP0,
     &                 'TABLE',I1,X1,CHARRE,.true.,ITMOD)
      if (ierr.ne.0) return
c     table itmod trouvee --> on recupere la force
      if (itmod.gt.0) then
        CALL ACCTAB(ITMOD,'MOT',0,0.d0,'FORCE',.true.,IP0,
     &                  'CHPOINT',I1,X1,CHARRE,.true.,ICHP0)
        if (ierr.ne.0) return
        if (ichp0.gt.0) then
          idemem(**) = ichp0
          NDEMEM=NDEMEM+1
*       write(6,*) ' extension idemem 3 ',idemem(/1)
          idnote(**) = im
        else
          call erreur(1)
          return
        endif
c       on cree le point repere ici
        CALL CREPO1 (ZERO, ZERO, ZERO, IPOIN)
        CALL ECCTAB(ITMOD,'MOT',0,0.0D0,'POINT_REPERE',.TRUE.,0,
     &        'POINT',0,0.0D0,' ',.TRUE.,IPOIN)
      endif
      GOTO 80
      IF (IERR.NE.0) RETURN
C----------------------------------------------------------------------
C     DEBUT DU TRAVAIL
 90   CONTINUE
      SEGINI IDEME0,IDEME1
C
C     Verifier qu'il n'y a pas de blocage en double
***   call verlag(ipoiri)
      if (ierr.ne.0) return
*     y a t il des matrices de relations non unilaterales
      mrigid=ipoiri
C     call prrigi(ipoiri,1)
      segact mrigid
      ifochs=iforig
      idepe=0
*     write(ioimp,*) 'dans resou mrigid iforig ',mrigid,iforig
C
      nbr = irigel(/2)
      if (jrcond.ne.0) nelim=30
      do 1000 irig = 1,nbr
       meleme=irigel(1,irig)
       segact meleme
       if ((irigel(6,irig).eq.0.or.nounil.eq.1).and.itypel.eq.22)
     >    idepe=idepe+num(/2)
       if (irigel(6,irig).ne.0) iunil=1
       if (irigel(6,irig).eq.2) nelim=30
       if (irigel(7,irig).ne.0) then
         insym=1
         ichski=0
       endif
1000  continue
C
C     Elimination recursive des conditions aux limites
*  on la fait en gradient conjugue ou en appel de unilater
      if (igradj.eq.1.or.(iunil.eq.1.and.nounil.eq.0)) nelim=30
      nfois=nelim-1
      bdblx=.false.
      imult=1
      icond=idepe
      icondi=(icond*10)/9+1
      if=0
      do ifois=1,nfois
        if(imult.ne.0.and.icond.ne.0.and.(icond*10)/9.lt.icondi.and.
     >      (icondi-icond.gt.0.or.igradj.eq.1)) then
          icondi=icond
          if=if+1
          if(ierr.ne.0) return
          call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
     >                nounil,bdblx,icond,imult,if,imtvid,nelim)
C         write(ioimp,*) ' passe ',if,' condition ',icond,ifois
          if(ierr.ne.0) return
          mrigid=mrigic
        endif
      enddo
C
C     S'il reste des conditions : dedoubler les mult de Lagrange restants
C     -> nouvel appel pour creer lagdua et adapter les seconds membres
      if (iunil.eq.0.or.nounil.eq.1) then
        if (icond.ne.0) then
          if=if+1
          bdblx=.true.
          if(ierr.ne.0) return
          call resouc(mrigid,mrigic,idemem,ideme0,ideme1,
     >                nounil,bdblx,icond,imult,if,imtvid,nelim)
*         write(ioimp,*) ' passe ','finale',' condition ',icond
          if(ierr.ne.0) return
          mrigid=mrigic
        endif
      endif
*     write (ioimp,*) 'nombre de passes',if,mrigid.imlag
      if (idepe.ne.0) noid=1
C
C     Rigidite reduite aux inconnues independantes
      ipoiri=mrigid
*     call prrigi(ipoiri,1)
C-------------------------------------------------------
*
* Si au moins une des matrices n'est pas symétrique, on passera
* par le solveur non-symétrique LDMT.
*
      SEGACT MRIGID*MOD
      NRG = IRIGEL(/1)
      NBR = IRIGEL(/2)
C ... Ceci peut arriver si par exemple on extrait la partie
C     symétrique d'une matrice purement antisymétrique ...
*     IF(NBR.EQ.0) THEN
*        SEGDES MRIGID
*        CALL ERREUR(727)
*        RETURN
*     ENDIF
C ... Mais avant on va tester si la normalisation des variables
C     primales et duales a été demandée - ceci entraîne la perte
C     de la symétrie ...
      IF(NORINC.GT.0  .AND. NORIND.GT.0) THEN
        IF(KSYMET.EQ.1) THEN
          CALL ERREUR(19)
          SEGDES,MRIGID
          RETURN
        ENDIF
        INSYM = 1
        IGRADJ = 0
        ICHSKI = 0
        GOTO 15
      ENDIF

C    ... On teste si la matrice contient des matrices non symétriques ...
C    ... Si OUI, ce n'est pas la peine de faire les autres tests ...
      DO  9 IN = 1,NBR
        IF(IRIGEL(7,IN).GT.0) THEN
C        ... On vérifie si l'utilisateur n'a pas demandé explicitement
C            la résolution par Choleski ou gradient conjugué,
C            si OUI on râle puis on s'en va !!! ...
          IF(KSYMET.EQ.1) THEN
            CALL ERREUR(1126)
            SEGDES,MRIGID
            RETURN
          ENDIF
          IF(NORINC.NE.0.AND.NORIND.EQ.0) THEN
*  on supprime la normalisation au lieu de faire une erreur
            norinc=0
**          CALL ERREUR(760)
**          SEGDES,MRIGID
**          RETURN
          ENDIF
          INSYM = 1
          IGRADJ = 0
          ICHSKI = 0
          GOTO 15
        ENDIF
  9   CONTINUE
 15   CONTINUE
C
C  ****  ON S'ASSURE QU'IL N'Y A PAS D'APPUIS UNILATERAUX
C
      IF (IUNIL.EQ.0.OR.NOUNIL.EQ.1) GOTO 30
C
C  **** EXISTENCE DES APPUIS UNILATERAUX
C  **** SI ON EST DEJA PASSE DANS LES APPUIS UNILATERAUX
C       ISUPEQ POINTE SUR UNE TABLE CONTENANT LES DONNEES A PASSER
C       A LA PROCEDURE UNILATER
C
      ISUPLO=ISUPEQ
      IF (ISUPLO.NE.0) GOTO 27
C
C     Distinguer ce qui est unilateral (RI2) du reste (RI1)
      NNOR=0
      DO 22 I=1,IRIGEL(/2)
        IF(IRIGEL(6,I).EQ.0) NNOR=NNOR+1
  22  CONTINUE
      IF(NNOR.EQ.0)  THEN
        CALL ERREUR(312)
        GOTO 5000
      ENDIF
C
      NRIGE=IRIGEL(/1)
      NRIGEL=NNOR
      SEGINI RI1
      RI1.IFORIG = IFORIG
c*      RI1.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
      RI1.MTYMAT = '        '
      NRIGEL=IRIGEL(/2)-NNOR
      SEGINI RI2
      RI2.IFORIG = IFORIG
c*      RI2.MTYMAT = MTYMAT <- type TEMPORAI(RE) plantage severe
      RI2.MTYMAT = '        '
      II1=0
      II2=0
      DO 23 I=1,IRIGEL(/2)
        IF(IRIGEL(6,I).NE.0) THEN
          RI3=RI2
          II2=II2+1
          II=II2
        ELSE
          RI3=RI1
          II1=II1+1
          II=II1
        ENDIF
        DO 24 J=1,NRIGE
          RI3.IRIGEL(J,II) = IRIGEL(J,I)
   24   CONTINUE
        RI3.COERIG(II)=COERIG(I)
   23 CONTINUE
C
C     Creation de la table a transmettre a UNILATER (stockee dans la
C     raideur initiale)
      CALL CRTABL(MTABLE)
      ISUPEQ=MTABLE
      MRIGID=IPRIGO
      SEGACT MRIGID*MOD
      ISUPEQ=MTABLE
      CALL ECCTAB(MTABLE,'ENTIER  ',1  ,XVA,' ',ILOG,IOB,
     $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI1)
      CALL ECCTAB(MTABLE,'ENTIER  ',2  ,XVA,' ',ILOG,IOB,
     $                   'RIGIDITE',IOB,XVA,' ',ILOG,RI2)
      CALL ECCTAB(MTABLE,'ENTIER  ',3  ,XVA,' ',ILOG,IOB,
     $                   'LOGIQUE ',IOB,XVA,' ',ILOG,IOB)

      ISUPLO=MTABLE
      SEGDES RI1,RI2,MTABLE
C
C     Recuperer la table a transmettre a UNILATER
  27  CONTINUE
      MTABLE=ISUPLO
      SEGACT MTABLE
      ILUG=(INSYM.EQ.1)
C
      CALL ECCTAB(MTABLE,'MOT     ',4  ,XVA,'NSYM',ILOG,IOB,
     $                   'LOGIQUE ',IOB,XVA,' '   ,ILUG,IOB)
      if(idepe.ne.0) then
*  on passe les ideme* a mrem sous forme de listchpo
        n1=if
        segini mlchpo,mlchp1
        do i=1,if
          mlchpo.ichpoi(i)=ideme0(1,i)
          mlchp1.ichpoi(i)=ideme1(1,i)
        enddo
        CALL ECCTAB(MTABLE,'ENTIER  ',10 ,XVA,' ',ILOG,IOB,
     $                     'LISTCHPO',IOB,XVA,' ',ILOG,mlchpo)
        CALL ECCTAB(MTABLE,'ENTIER  ',11 ,XVA,' ',ILOG,IOB,
     $                     'LISTCHPO',IOB,XVA,' ',ILOG,mlchp1)
*  pour mrem on met la derniere raideur condensee. Elle contient les pointeurs pour remonter
        CALL ECCTAB(MTABLE,'ENTIER  ',50 ,XVA,' ',ILOG,IOB,
     $                     'RIGIDITE',IOB,XVA,' ',ILOG,ipoiri)
      endif
      SEGDES MRIGID
      DO 26 I=NDEMEM,1,-1
        ISECO=IDEMEM(I)
        CALL ACTOBJ ('CHPOINT ',ISECO,1)
        CALL ECROBJ ('CHPOINT ',ISECO)
  26  CONTINUE
      SEGSUP IDEMEM
      CALL ECROBJ ('TABLE   ',ISUPLO)
      SEGINI MTEXTE
      LTT=8
      MTEXT(1:LTT) ='UNILATER'
      NCART=8
      SEGDES MTEXTE
      CALL ECROBJ('TEXTE',MTEXTE)
      mrigid=iprigo
      segdes mrigid
      GOTO 5000
C
C ... On arrive ici dans le cas où il n'y a pas d'appuis unilatéraux ...
   30 CONTINUE
* il se peut que le dernier chp soit du frottement
* on l'enleve car il ne sert a rien si on n'appele pas unilater
      if (NDEMEM.gt.1.and.idepe.ne.0) then
        mchpoi=ideme0(NDEMEM,if)
        segact MCHPOI
        if (mtypoi.eq.'LX  ') NDEMEM=NDEMEM-1
      endif
C
*     write(6,*) ' ichski, igradj,insym ',ichski, igradj,insym
*     write (6,*) ' imtvid ',imtvid
C
C     Cas matrice vide
      if (imtvid.eq.1) then
***     write(6,*) ' attention matrice vide. Système surcontraint '
        if (nounil.eq.0) call erreur(-364)
*
        nsoupo=0
        nat=0
        segact idemem*mod

        do i=1,idemem(/1)
          segini mchpoi
          mchpoi.ifopoi = ifochs
          idemem(i)=mchpoi
        enddo
        if (noen.eq.0) then
          nat=2
          nsoupo=0
          segini mchpo4
          call ecrobj('CHPOINT ',mchpo4)
          call ecrent(0)
        endif
      else
*       write(6,*) ' appel resou1 --  idemem(1)'
*       segact idemem
*       idesec= idemem(1)
*       call ecchpo(idesec,0)
*       write(6,*) ' appel resou1 -- ipoiri'
*       call prrigi ( ipoiri,1)
*       write(6,*) ' ichski  insym  ', ichski, insym
        IF(ICHSKI.EQ.1) CALL RESOU1(IPOIRI,IDEMEM,NOID,NOEN,PREC,
     >          ISTAB,ISOUCI)
        IF(IGRADJ.EQ.1) CALL GRACO0(IPOIRI,IDEMEM,NOID,NOEN)
        IF(INSYM .EQ.1) CALL LDMT  (IPOIRI,IDEMEM,NOID,NOEN,PREC,ISOUCI)
        IF(IERR.NE.0) GO TO 5000
      ENDIF
C
C     -----------------------------------------------------------------
C     Reintroduire les inconnues eliminees
      do 2010 ifois=1,30
        segact mrigid
        mrigid=jrsup
        if (mrigid.eq.0) goto 2011
        if(ierr.ne.0) GOTO 5000
        call resour(idemem,ideme0,ideme1,mrigid,if,ipt8,isouci,iverif)
        if(ierr.ne.0) GOTO 5000
        if=if-1
 2010 continue
 2011 continue
C
C     -----------------------------------------------------------------
C     ECRITURE DES OBJETS RESULTATS
      SEGACT IDEMEM
      do 3 i=1,NDEMEM
       il = NDEMEM + 1 - i
       iret=idemem(il)
C
       mchpoi=iret
       segact mchpoi*mod
       mchpoi.jattri(1)=1
C
      if (itbas.ne.0) then
         ilo = idnote(il)
         CALL ACCTAB(ITBAS,'ENTIER',ILO,0.d0,' ',.true.,IP0,
     &                   'TABLE',I1,X1,CHARRE,.true.,ITMOD)
         if (ierr.ne.0) GOTO 5000

         CALL ECCTAB(ITMOD,'MOT',0,0.D0,'DEFORMEE',
     &        .TRUE.,0,'CHPOINT',0,0.D0,' ',.TRUE.,IRET)

         if (i.EQ.NDEMEM) then
           segdes mtab1
           segsup idnote
           CALL ECROBJ ('TABLE   ',itbas)
        endif
      else if (ipshpo.gt.0) then
         mlchpo = ipshpo
         ichpoi(il) = iret
         if (i.EQ.NDEMEM) then
          mlchpo = ipshpo
          CALL ACTOBJ ('LISTCHPO ',ipshpo,1)
          CALL ECROBJ ('LISTCHPO ',ipshpo)
        endif
      else
         CALL ACTOBJ ('CHPOINT ',IRET,1)
         CALL ECROBJ ('CHPOINT ',IRET)
      endif

   3  CONTINUE
      SEGSUP IDEMEM
C
 5000 CONTINUE
      norinc=norins
      END