C MONDES    SOURCE    PV        11/05/03    21:17:44     6959           
      SUBROUTINE MONDES(MMATRX,MVECTX,NOEN)
C
C  **** EXECUTE LA SOLUTION X DE  (Lt D L) X=F
C
CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
CMB
CMB   Plutot la solution de L.D.Lt ou L.D.Mt (cas non symétrique)
CMB   Elle devrait dons s'appeller DESMON et non MONDES.
CMB
CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (LPCRAY=10000)
      INTEGER OOOVAL

C ... Variable décrivant l'EXIStence du Triangle Supérieur différent
C     de l'inférieur (cas non symétrique) ...
      LOGICAL EXISTS

-INC SMMATRI
-INC SMELEME
-INC SMVECTD
-INC CCOPTIO

      SEGMENT ITTRV(MULRE)
      SEGMENT,ITRAA(NENS)
      SEGMENT/BID/(BIDON(IIMAX+10)*D)

      if (noen.eq.0) then
       nbelem=0
       nbnn=1
       nbsous=0
       nbref=0
       segini ipt8
       ipt8.itypel=1
      endif
      CALL OOOMRU(1)
      IVALMA=0
      IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX
 1000 FORMAT('  SUBROUTINE  MONDES : POINTEUR DE LA MATRICE=',I5,
     1                          '    POINTEUR DU VECTEUR=',I5)
C
C  ****  ACTIVATION DES SEGMENTS
C
      MMATRI=MMATRX
      SEGACT,MMATRI*MOD
*  BEC00SE OPTIMISEUR
      ITRAA=MMATRI
      IF(NENS.NE.0) SEGINI,ITRAA
      NE1=NENS
      IAA=0

      EXISTS=.FALSE.
      IF(IILIGS.NE.0) THEN
         MILIGN=IILIGS
         SEGACT MILIGN
         IF(ILIGN(/1).GT.0) EXISTS=.TRUE.
      ENDIF

      MILIGN=IILIGN
      SEGACT,MILIGN
      INC=IPNO(/1)

      MVECTD=MVECTX
      SEGACT MVECTD*MOD

      MDNOR=IDNORM
      SEGACT MDNOR
      IF(IDNORD.GT.0) THEN
         MDNO1=IDNORD
         SEGACT MDNO1
      ELSE
         MDNO1=MDNOR
      ENDIF

      MDIAG=IDIAG
      SEGACT,MDIAG

C ... MULRE = nombre de seconds membres ...
      MULRE = VECTBB(/1)/INC
C ... MUINC ne servira que comme une borne sur les indices de boucles ...
      MUINC = ( MULRE-1)*INC
      SEGINI ITTRV

C ... Multiplication du second membre par les DNOR(.) (à cause du
C     conditionnement des matrices de blocages) ...
      DNORMA=0.D0
      DO 450 K=0,MUINC,INC
         DO 45 I=1,INC
            VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I)
            DNORMA=MAX(DNORMA,ABS(VECTBB(I+K)))
 45      CONTINUE
 450  CONTINUE
      DNORMA=DNORMA*1D-18

C ... Détection du premier terme important du second membre ...
      II=0
      DO 451 K=0,MUINC,INC
         II=II+1
         DO 452 I=1,INC
            IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452
C       ... Le numéro du noeud concerné par le premier terme important va à ITTRV ...
            ITTRV(II)=IPNO(I)
            GO TO 451
  452    CONTINUE
  451 CONTINUE

C ... NNOE = nombre de noeuds concernés ...
      NNOE=ILIGN(/1)
      IDEP=NNOE
C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ...
      DO 453 I=1,MULRE
         IF(ITTRV(I).LT.IDEP) THEN
            IF(ITTRV(I).NE.0) THEN
               IDEP=ITTRV(I)
            ELSE
               IDEP=1
               GO TO 4530
            ENDIF
         ENDIF
 453  CONTINUE
 4530 CONTINUE
      SEGDES MDNOR
C
C  ****  DESCENTE: ON RESOU L*C=B.  EN FAIT ON STOCKE C DANS B
C
      LTRK=MAX(1+LPCRAY,OOOVAL(1,4))
      IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK
*  test si la place disponible permet de tout mettre en memoire
      xplds=oooval(1,1)-oooval(2,3)
      if (real(nnoe)*ijmax.lt.xplds) iimax=0
      SEGINI BID
C ... IDEP : on commence par ce noeud car le second membre est nul pour
C     tous les précédents ...
      DO 610 IVS=IDEP,NNOE
         LIGN=ILIGN(IVS)
         SEGACT /ERR=50/ LIGN
         IVALMA=IVALMA+VAL(/1)
*pvpv         IF (IVALMA.GT.NGMAXY) GOTO 50
         NA=IMMM(/1)
         IPRELL=IPREL
         DO 611 IA=1,NA
C       ... Si l'inconnue présente un mode d'ensemble ...
            IF(IMMM(IA).NE.0) THEN
C          ... On incrémente le compteur et
               IAA=IAA+1
C          ... on met le N° de la ligne dans ITRAA ...
               ITRAA(IAA)=IPRELL+IA-1
            ENDIF
  611    CONTINUE
C    ... IFIB sera vu dans MONDE2 par le COMMON/CMOND2/ ...
         IFIB=IVPO(/1)
         CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
     >      NA,IPREL,MULRE,INC,IVS,LLOM,IFIB)
 610  CONTINUE
C ... Si on n'a pas eu de pb de mémoire on passe par là ...
      SEGSUP BID
C ... ILMAX vaut le dernier noeud qui a pu être traité ...
      ILMAX=NNOE
C ... On va à la division par la diagonale ...
      GOTO 55

  50  CONTINUE
C ... Si on est là, c'est à cause des pb de mémoire ...
      SEGSUP BID
C ... IVS est le N° du LIGN qui n'a pas pu être traité ...
C ... ILMAX = N° du dernier traité ...
      ILMAX=IVS-1
C ... IILMAX = N° du premier à traiter ...
      IILMAX=IVS
      DO 54 IVS=IILMAX,NNOE
  58     CONTINUE
         LIGN=ILIGN(IVS)
C    ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ...
         SEGACT /ERR=56/ LIGN
         GOTO 57
  56     CONTINUE
*  EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE
C    ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX,
         LIGN=ILIGN(ILMAX)
         SEGDES LIGN*(NOMOD,MRU)
C    ... puis ce dernier est décrémenté ...
         ILMAX=ILMAX-1
         IF (ILMAX.EQ.1) CALL ERREUR(5)
C    ... et on recommence ...
         GOTO 58
C    ... On est là si tout s'est bien passé avec SEGACT ...
  57     CONTINUE
         NA=IMMM(/1)
         IPRELL=IPREL
         DO 612 IA=1,NA
            IF(IMMM(IA).EQ.0) GOTO 612
            IAA=IAA+1
            ITRAA(IAA)=IPRELL+IA-1
  612    CONTINUE
         IFIB = IVPO(/1)
         CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
     >      NA,IPREL,MULRE,INC,IVS,LLOM,IFIB)
         SEGDES,LIGN*(NOMOD,MRU)
   54 CONTINUE
   55 CONTINUE
C
C ... À cet endroit la descente est finie. L'état des LIGN est le suivant :
C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ...
C
C  ****  DIVISION PAR LE TERME DIAGONAL  ****
C
      DO 120 K=0,MUINC,INC
         DO 12 I=1,INC
            J=I+K
            VECTBB(J)=VECTBB(J)/DIAG(I)
   12    CONTINUE
  120 CONTINUE
      SEGDES,MDIAG

C
C  ****  MONTEE  ****
C
C ... Si la matrice est asymétrique ...
      IF(EXISTS) THEN

C    ... On commence par désactiver les LIGN qui sont restés actifs ...
         DO 2000 I=IDEP,ILMAX
            LIGN=ILIGN(I)
            SEGDES,LIGN*(NOMOD,MRU)
 2000    CONTINUE

C    ... Puis on désactive MILIGN ...
         SEGDES,MILIGN

C    ... On change de MILIGN ...
         MILIGN=IILIGS
         SEGACT,MILIGN

C    ... Et on active des LIGN pointés par ce dernier ...
         ILMAX=IDEP-1
         DO 2100 I=IDEP,NNOE
            LIGN=ILIGN(I)
            SEGACT /ERR=2110/ LIGN
            ILMAX = I
 2100    CONTINUE
C    ... Et on passe à la montée ...
         GOTO 3000

C    ... Si on n'a même pas activé le premier, alors ADIOS ...
 2110    IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5)

 3000    CONTINUE
      ENDIF


      J=NNOE
C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ...
  70  CONTINUE
         IF (J.EQ.ILMAX) GOTO 72
         LIGN=ILIGN(J)
         SEGACT /ERR=68/ LIGN
         GO TO 67
C    ... Si on a des pb avec activation, on désactive et on diminue ILMAX ...
 68      CONTINUE
         LIGN = ILIGN( ILMAX)
         SEGDES LIGN*(NOMOD,MRU)
         ILMAX=ILMAX-1
         IF(ILMAX.EQ.1) CALL ERREUR (5)
C    ... puis on recommence la tentative ...
         GO TO 70
C    ... On a réussi à activer ...
  67     CONTINUE
         NA=IMMM(/1)
         IPRELL=IPREL
         IFIB=IVPO(/1)
         CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
     >     NA,IPREL,MULRE,INC,IFIB)
         SEGDES,LIGN*(NOMOD,MRU)
         J = J-1
         GO TO 70
  72  CONTINUE
CC  FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!!
CC  FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1

C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ...
      DO 13 J=ILMAX,IDEP,-1
         LIGN=ILIGN(J)
         NA=IMMM(/1)
         IPRELL=IPREL
         IFIB=IVPO(/1)
         CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
     >     NA,IPREL,MULRE,INC,IFIB)
         SEGDES,LIGN*(NOMOD,MRU)
   13 CONTINUE

C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ...
      DO 113 J=IDEP-1,1,-1
         LIGN=ILIGN(J)
         SEGACT LIGN
         NA=IMMM(/1)
         IPRELL=IPREL
         DO 1131 ILM=1,NA
            IF(IMMM(ILM).EQ.0)  GO TO 1131
            IAA=IAA+1
            ITRAA(IAA)=IPRELL+ILM-1
 1131    CONTINUE
         IFIB=IVPO(/1)
         CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
     >     NA,IPREL,MULRE,INC,IFIB)
         SEGDES,LIGN*(NOMOD,MRU)
  113 CONTINUE

C
C  ****  ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE
C
C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ...
      NENS=IAA
      NBEN=0
      IF(NENS.EQ.0)  GO TO 26

      MINCPO=IINCPO
      MIMIK=IIMIK
      SEGACT,MINCPO,MIMIK
      MELEME=IGEOMA
      SEGACT,MELEME
      NNOE=INCPO(/2)
      IINC1=INCPO(/1)

C ... Boucle sur les seconds membres ...
      DO 300 KI=0,MUINC,INC
C    ... XMA sera le maximum des val. abs. des termes (multiplicateurs 
C        compris) du résultat partiel (avant la multiplication par MDNOR)
C        N° KI+1 multiplié par 1.e-10 ...
         XMA=0.
         DO 30 KK=1,INC
            I=KK+KI
            XMA=MAX(XMA,ABS(VECTBB(I)))
   30    CONTINUE
         XMA=XMA*1.E-10

C    ... Boucle sur les modes d'ensemble ...
         DO 21 IA=1,NENS
            I1=ITRAA(IA)
            J=IPNO(I1)
            DO 22 K=1,IINC1
               IF(INCPO(K,J).EQ.I1) GO TO 23
   22       CONTINUE
               CALL ERREUR(5)
               RETURN

   23       CONTINUE
C       ... Si ce n'est pas un multiplicateur, le déplacement doit être
C           < à XMA, sinon ERREUR ...
            IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA    ) then
              vectbb(i1+ki)=0.d0
              GO TO 20
            endif
C       ... Si c'est un multiplicateur, le multiplicateur doit être
C           < à 1.e+7*XMA, sinon ERREUR ...
            if (ittr(i1).ne.0) then
             i2=ittr(i1)
             if (abs(vectbb(i1+ki)+vectbb(i2+ki)).le.xma    ) then
              vectbb(i1+ki)=0.d0
              vectbb(i1+ki)=0.d0
              goto 20
             endif
            endif

*           IF(ITTR(I1).NE.0.AND.ABS(VECTBB(I1+KI)).LE.XMA*1E4) then
*             vectbb(i1+ki)=0.d0
*             GO TO 20
*           endif
C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
C ... Si option NOEN on accumule les pts dans un meleme
            IF(NOEN.EQ.0) THEN
             nbelem=ipt8.num(/2)+1
             segadj ipt8
             ipt8.num(1,nbelem)=num(1,j)
            ENDIF
               NBEN=NBEN+1
               MOTERR(1:4)=IMIK(K)
               INTERR(1)=NUM(1,J)
               IF(NOEN.EQ.0) THEN
                 GO TO 21
               ENDIF
               ITY=149
               CALL ERREUR(ITY)
               RETURN

   20       CONTINUE
C       ... Messages d'information ...
            JJK=NUM(1,J)
            KIK=KI/INC +1

            IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1) then
              if (imik(k).ne.'LX  ')
     >                WRITE(IOIMP,41) JJK,IMIK(K)
            endif

            IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1)
     &                WRITE(IOIMP,42) JJK,IMIK(K)

            IF(ITTR(I1).NE.1 .AND. MULRE.NE.1)
     &                WRITE(IOIMP,410)KIK,JJK,IMIK(K)

            IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
     &                WRITE(IOIMP,420)KIK,JJK,IMIK(K)

   21    CONTINUE

  300 CONTINUE

   41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
     *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
     * 'LA SUSDITE')
   42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
     * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
     * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE ')
  410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
     *,I6,' INCONNUE ',
     *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
     * 'LA SUSDITE')
  420 FORMAT(' VECTEUR NUMERO ',I3,/,
     * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
     * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
     * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE ')

      SEGDES,MELEME
      SEGDES,MINCPO
      SEGDES,MIMIK
   26 CONTINUE
*  indiquer si besoin le nombre de modes d'ensemble excites
*  et le maillage des noeuds freinés
      IF (NOEN.EQ.0) then
*      write (6,*) ' mondes nben ',nben
       CALL ECRENT(NBEN)
       call ecrobj('MAILLAGE',ipt8)
       segdes ipt8
      endif
      MDNOR=IDNORM
      SEGACT,MDNOR
      DO 350 K=0,MUINC,INC
         DO 35 I = 1, INC
            VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
   35    CONTINUE
*  on force l'egalite des multiplicateurs de lagrange
         DO 36 I = 1, INC
            if (ITTR(I).ne.0) then
*            write (6,*) ' dans mondes ',i,ittr(i)
             if (vectbb(i+k).eq.0.d0.or.vectbb(ittr(i)+k).eq.0.d0) then
*      write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k)
             vectbb(i+k)=0.d0
             vectbb(ittr(i)+k)=0.d0
             goto 36
             endif
             vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))/2
             vectbb(ittr(i)+k)=vectbb(i+k)
            endif
   36    CONTINUE
  350 CONTINUE
      SEGDES,MDNOR
C ... On ne désactive MDNO1 que dans le cas où il est
C     différent de MDNOR ...
      IF(IDNORD.GT.0) THEN
         SEGDES,MDNO1
      ENDIF
      SEGDES MMATRI
      SEGDES,MILIGN
      SEGDES,MVECTD
      SEGSUP ITTRV
      IF(NE1.NE.0) SEGSUP,ITRAA

      IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
 1001 FORMAT('  SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)

      CALL OOOMRU(0)
      RETURN
      END
C