C GAUNEW    SOURCE    CHAT      09/10/09    21:18:31     6519           
      SUBROUTINE GAUNEW
C
C     Determination de la direction de descente pour un probleme
C        de moindre carre par la methode de Gauss-Newton
C
C     Possiblilite d'introduire un parametre de viscosite (facultatif)
C
C   TAB2=GANE TAB1 ('AMOR' FLOT1);
C   
C   CHPO1 RIGI1=GANE TAB1 'MATR' ('AMOR' FLOT1);
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCOPTIO
-INC SMRIGID
-INC SMCHPOI
-INC SMELEME
-INC CCHAMP
-INC SMLENTI
-INC SMLREEL
-INC SMCOORD
C
      CHARACTER*8 TYPOBJ,CHARR0
      CHARACTER*30 CHARRE
      REAL*8 X0,X1
      LOGICAL L0,L1
*
      PARAMETER(NOPTIO=2)
      CHARACTER*4 OPTION(NOPTIO)
      LOGICAL LMATR
      DATA OPTION/'AMOR','MATR'/
C
C     LECTURE DES OPTIONS
C
      LMATR=.FALSE.
 1    CALL LIRMOT(OPTION,NOPTIO,IRETOU,0)
      IF(IRETOU.EQ.0)GOTO 9
      GOTO(2,3),IRETOU
C
C     OPTION AMOR: On introduit un amortissement
C
 2    CALL LIRREE(AMOR,1,IRETOU)
      IF(IERR.NE.0)RETURN
      GOTO 1
C
C     OPTION MATR: Osort avec la matrice et le second membre
C
 3    LMATR=.TRUE.
      GOTO 1
C
 9    CONTINUE
C
C     LECTURE DE LA TABLE
C
      CALL LIROBJ('TABLE',MTABLE,1,IRETOU)
      IF(IERR.NE.0)RETURN
C
C     INDEX TAB1.SOUSTYPE
C
      CALL ACCTAB(MTABLE,'MOT',I0,X0,'SOUSTYPE',L0,IPO0,
     $                   'MOT',I1,X1,CHARRE,L1,IPO1)
      IF(IERR.NE.0) THEN
        IRET=1
        RETURN
      ENDIF
      IF(CHARRE(1:7).NE.'VECTEUR')THEN
        CALL ERREUR(648)
        IRET=1
        RETURN
      ENDIF
C
C     CHPOINT INDEX 0 (valeur de la fonction)
C     
      CALL ACCTAB(MTABLE,'ENTIER',0,X0,'SOUSTYPE',L0,IPO0,
     $                   'CHPOINT',I1,X1,CHARRE,L1,MCHPOI)
      IF(IERR.NE.0) THEN
        IRET=1
        RETURN
      ENDIF    
C
C     VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
C     SANS SERIE DE FOURIER
C
      SEGACT,MCHPOI
      IF(IPCHP(/1).NE.1)THEN
        WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
        MOTERR(1:8)='CHPOINT '
        INTERR(1)=MCHPOI
        CALL ERREUR(110)
        GOTO 9999
      ENDIF
      IF(IFOPOI.EQ.1)THEN
        MOTERR(1:8)='GANE    '
        CALL ERREUR(333)
        GOTO 9999
      ENDIF
C
C     ACTIVATIONS DIVERSES pour la fonction
C
      MSOUPO=IPCHP(1)
      SEGACT,MSOUPO
      MELEME=IGEOC
      MPOVAL=IPOVAL
      SEGACT,MELEME,MPOVAL
      N1=VPOCHA(/1)
      NC1=VPOCHA(/2)
C
C     ON PASSE AUX AUTRES INDEX ENTIERS DE TYPE CHPO CONSECUTIFS DE LA TABLE
C
      JG=0
      SEGINI,MLENTI
10    CONTINUE
        TYPOBJ=' '
        JG=JG+1
        CALL ACCTAB(MTABLE,'ENTIER',JG,X0,CHARR0,L0,IPO0,
     $                       TYPOBJ,I1,X1,CHARRE,L1,IPO1)
        IF(TYPOBJ.EQ.'CHPOINT')THEN
          SEGADJ,MLENTI
          LECT(JG)=IPO1
          GOTO 10  
        ENDIF 
      CONTINUE
C
C     ... ET ON VERIFIE QUE TOUT BAIGNE
C  
      JG=LECT(/1)
      IF(JG.EQ.0)THEN
        WRITE(IOIMP,*)'GAuss NEwton: The gradient is empty'
        CALL ERREUR(314)
        GOTO 9998
      ENDIF
C
      DO IE1=1,JG
        MCHPO1=LECT(IE1)
C
C     --> VERIFICATION DU CARACTERE ELEMENTAIRE DU CHPOINT,
C         SANS SERIE DE FOURIER
C
        SEGACT,MCHPO1
        IF(MCHPO1.IPCHP(/1).NE.1)THEN
          WRITE(IOIMP,*)'GAuss NEwton: work only with elementary CHPO'
          MOTERR(1:8)='CHPOINT '
          INTERR(1)=MCHPOI
          CALL ERREUR(110)
          SEGDES,MCHPO1
          GOTO 9998
        ENDIF
        IF(MCHPO1.IFOPOI.EQ.1)THEN
          MOTERR(1:8)='GANE    '
          CALL ERREUR(333)
          SEGDES,MCHPO1
          GOTO 9998
        ENDIF        
C
C     --> MEME COMPOSANTES, MEME MAILLAGE
C
        MSOUP1=MCHPO1.IPCHP(1)
        SEGACT,MSOUP1
C
        IF (NC1.NE.MSOUP1.NOHARM(/1))THEN
          WRITE(IOIMP,*)
     >    'GAuss NEwton: all CHPO should have the same nb of components'
          CALL ERREUR(665)
          SEGDES,MSOUP1,MCHPO1
          GOTO 9998         
        ENDIF
C
        DO IE2=1,NC1
          IF(NOCOMP(IE2).NE.MSOUP1.NOCOMP(IE2))THEN
          WRITE(IOIMP,*)
     >    'GAuss NEwton: all CHPO should have the same components',
     >    ' in the same order'
            CALL ERREUR(488)
            SEGDES,MSOUP1,MCHPO1
            GOTO 9998
          ENDIF
        ENDDO
C
        IPT1=MSOUP1.IGEOC
        SEGACT,IPT1
        IF(N1.NE.IPT1.ICOLOR(/1))THEN
          WRITE(IOIMP,*)
     >    'GAuss NEwton: all CHPO should have the same number',
     >    ' of points'
          MOTERR(1:8) ='MAILLAGE'
          MOTERR(9:16)='NOEUD    '
          CALL ERREUR(403)
          SEGDES,MSOUP1,MCHPO1,IPT1
          GOTO 9998
        ENDIF
C
        DO IE2=1,N1
          IF(NUM(1,IE2).NE.IPT1.NUM(1,IE2))THEN
            WRITE(IOIMP,*)
     >      'GAuss NEwton: all CHPO should have the same',
     >      ' points in the same order'
            CALL ERREUR(403)
            SEGDES,MSOUP1,MCHPO1,IPT1
            GOTO 9998
          ENDIF
        ENDDO
C
        LECT(IE1)=MSOUP1.IPOVAL
        SEGDES,MSOUP1,MCHPO1,IPT1
C
      ENDDO
C
C     ON PREPARE LA MATRICE ET LE SECOND MEMBRE
C
C     1) NOEUDS DU SUPER ELEMENT ET DU CHPO
C
      NBPTT=XCOOR(/1)/(IDIM+1)
      NBPTS=NBPTT+JG
      SEGADJ,MCOORD
C
C     2) SUPER ELEMENT
C
      NBSOUS=0
      NBELEM=1
      NBNN=JG
      NBREF=0
      SEGINI,IPT1
      IPT1.ITYPEL=28
      DO IE1=1,JG
        IPT1.NUM(IE1,1)=NBPTT+ IE1
      ENDDO
      SEGDES,IPT1
C
C     4) DECRIPTEUR POUR LA RIGIDITE
C
      NLIGRP=JG
      NLIGRD=JG
      SEGINI,DESCR
      DO IE1=1,JG
        LISINC(IE1)=NOMDD(11)
        LISDUA(IE1)=NOMDU(11)
        NOELEP(IE1)=IE1
        NOELED(IE1)=IE1
      ENDDO
      SEGDES,DESCR
C
C     5) CONTENU DE LA RIGIDITE
C
      nelrig=1
      SEGINI,XMATRI
      DO IE1=1,NLIGRP
        MPOVA1=LECT(IE1)
        SEGACT,MPOVA1
        DO IE2=1,IE1
          MPOVA2=LECT(IE2)
          SEGACT,MPOVA2
C
          XVAL=0.D0
          DO IE3=1,NC1
            DO IE4=1,N1
              XVAL=XVAL+MPOVA1.VPOCHA(IE4,IE3)*MPOVA2.VPOCHA(IE4,IE3)
            ENDDO
          ENDDO
          RE(IE1,IE2,1)=XVAL
          RE(IE2,IE1,1)=XVAL
C
          SEGDES,MPOVA2
        ENDDO
      ENDDO
C
C     6) CORRECTION EVENTUELLE DE RIGIDITE SUR LA DIAGONALE
C
      IF(AMOR.NE.0.D0)THEN
        DO IE1=1,NLIGRP
          RE(IE1,IE1,1)=RE(IE1,IE1,1)+AMOR
        ENDDO
      ENDIF
      SEGDES,XMATRI
C
C     7) IMATRI
C
*      NELRIG=1
*      SEGINI,IMATRI
*      IMATTT(1)=XMATRI
*      SEGDES,IMATRI
C
C     8) CHAPEAU MRIGID DE LA RIGIDITE
C
      NRIGEL=1
      NRIGE=8
      SEGINI,MRIGID
      MTYMAT='RIGIGITE'
      COERIG(1)=1.D0
      IRIGEL(1,1)=IPT1
      IRIGEL(2,1)=0
      IRIGEL(3,1)=DESCR
      IRIGEL(4,1)=xMATRI
      IRIGEL(5,1)=0
      IRIGEL(6,1)=0
      IRIGEL(7,1)=0
      IRIGEL(8,1)=0
      ICHOLE=0
      IMGEO1=0
      IMGEO2=0
      IFORIG=IFOPOI
      ISUPEQ=0
      SEGDES,MRIGID
C
C     9) SUPPORT DU CHPO
C
      NBSOUS=0
      NBELEM=JG
      NBNN=1
      NBREF=0
      SEGINI,IPT1
      IPT1.ITYPEL=1
      DO IE1=1,JG
        IPT1.NUM(1,IE1)=NBPTT+ IE1
      ENDDO
      SEGDES,IPT1
C
C     10) VALEUR DU CHPO
C
      N=JG
      NC=1
      SEGINI,MPOVA2
C
      DO IE1=1,N
        MPOVA1=LECT(IE1)
        SEGACT,MPOVA1
C
        XVAL=0.D0
        DO IE2=1,NC1
          DO IE3=1,N1
            XVAL=XVAL-VPOCHA(IE3,IE2)*MPOVA1.VPOCHA(IE3,IE2)
          ENDDO
        ENDDO
        MPOVA2.VPOCHA(IE1,1)=XVAL
C
        SEGDES,MPOVA1
      ENDDO
C
      SEGDES,MPOVA2
C
C     11) MSOUPO
C
      SEGINI,MSOUP1
      MSOUP1.NOCOMP(1)=NOMDU(11)
      MSOUP1.IGEOC=IPT1
      MSOUP1.IPOVAL=MPOVA2
      SEGDES,MSOUP1
C
C     12) MCHPOI
C
      NAT=2
      NSOUPO=1
      SEGINI,MCHPO1
      MCHPO1.JATTRI(1)=JATTRI(1)
      MCHPO1.IFOPOI=IFOPOI
      MCHPO1.IPCHP(1)=MSOUP1
      SEGDES,MCHPO1
C
C     RETOUR A GIBIANE
C
      CALL ECROBJ('RIGIDITE',MRIGID)
      CALL ECROBJ('CHPOINT ',MCHPO1)
C
C     ON SORT SI LAMATR
C
      IF(LMATR)GOTO 9998
C
C     DESACTIVATIONS
C
      SEGDES,MELEME,MPOVAL,MSOUPO
      SEGSUP,MLENTI
      SEGDES,MCHPOI
C
C     ON APPELLE RESO
C
      CALL RESOU
C
C     ... ET ON REGARDE LE RESULTAT
C
      CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU)
      IF(IERR.NE.0)RETURN
C
      SEGACT,MCHPOI
      MSOUPO=IPCHP(1)
      SEGACT,MSOUPO
      MELEME=IGEOC
      MPOVAL=IPOVAL
      SEGACT,MELEME,MPOVAL
C
      SEGINI,MLENTI,MLENT1
      DO IE1=1,JG
        LECT(IE1)=NUM(1,IE1)
        MLENT1.LECT(IE1)=IE1
      ENDDO
      CALL GENOR2(LECT(1),MLENT1.LECT(1),JG)
C
      SEGINI,MLREEL
      DO IE1=1,JG
        PROG(IE1)=VPOCHA(MLENT1.LECT(IE1),1)
      ENDDO
      CALL VECTAB(MLREEL,JG,MTABLE)
C
      CALL ECROBJ('TABLE   ',MTABLE)
C
C     DESTRUCTION
C
      SEGSUP,MLENTI,MLENT1,MLREEL
C
      SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
C
      MCHPOI=MCHPO1
      SEGACT,MCHPOI
      MSOUPO=IPCHP(1)
      SEGACT,MSOUPO
      MELEME=IGEOC
      MPOVAL=IPOVAL
      SEGSUP,MELEME,MPOVAL,MSOUPO,MCHPOI
C
      SEGACT,MRIGID
      MELEME=IRIGEL(1,1)
      DESCR=IRIGEL(3,1)
      xMATRI=IRIGEL(4,1)
*      SEGACT,IMATRI
*      XMATRI=IMATTT(1)
      SEGSUP,MELEME,DESCR,XMATRI,MRIGID
C
      NBPTS=NBPTT
      SEGADJ,MCOORD
C   
      RETURN
C
C     ERREURS
C
 9998 SEGDES,MELEME,MPOVAL,MSOUPO
      SEGSUP,MLENTI
 9999 SEGDES,MCHPOI
      RETURN
      END