Sources Esope | Cast3M

levmar
C LEVMAR    SOURCE    KICH      16/07/12    21:15:04     9027           
      SUBROUTINE LEVMAR
*
* moteur pour algorithme levenberg-marquardt
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMLENTI
-INC SMLREEL
      SEGMENT MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL,MLREE8.MLREEL,
     +        MLREE9.MLREEL,MLRE11.MLREEL,MLRE12.MLREEL,MLRE13.MLREEL
 
      PARAMETER (NTYP=11,NPRO = 10)
      CHARACTER*4 LISTYP(NTYP)
      integer iprocm(NPRO,2)
      real*8 xmordo(NPRO)
      logical dpro
      DATA LISTYP / 'ABSC', 'ORDO', 'PARA', 'TORD', 'DYDA',
     &              'SIGM', 'PROC', 'ALPH', 'BETA', 'CHI2',
     &              'PARO' /
 
      data iprocm/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
      SAVE iprocm
      data xmordo/1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0,1.D0/
      SAVE xmordo
c
      dpro =.true.
      lapro = 0
c
      jg = NTYP
      segini mlenti
c
        call LIROBJ('LISTREEL',ip1,0,iret)
        IF (IERR.NE.0) GOTO 900
       if(iret.eq.1) then
        lect(4) = ip1
        call LIROBJ('LISTREEL',ip1,1,iret)
        IF (IERR.NE.0) GOTO 900
        lect(5) = ip1
       endif
c
      idon = 0
 
 10   continue
      CALL LIRMOT(LISTYP,NTYP,iplac,0)
      IF (IERR.NE.0) GOTO 900
      if (iplac.gt.0) then
        if (iplac.ne.7) then
           call LIROBJ('LISTREEL',ip1,1,iret)
        else if (iplac.eq.7) then
           call LIROBJ('PROCEDUR',ip1,1,iret)
           do kpro = 1,NPRO
            if (iprocm(kpro,1).eq.ip1) then
              iprocc = iprocm(kpro,2)
              lapro = kpro
              dpro = .false.
            endif
            if (iprocm(kpro,1).eq.0.and.dpro) then
              iprocm(kpro,1) = ip1
              iprocc = 0
              lapro = kpro
              dpro = .false.
            endif 
           enddo
        endif
        IF (IERR.NE.0) GOTO 900
        lect(iplac) = ip1
        idon = idon + 1
        goto 10
      endif
 
 
 
      if (iprocc.eq.0) then
        if (lect(7).le.0) then
          moterr(1:8) ='PROCEDUR'
          call erreur(37)
          GOTO 900
        endif
        mlree1 = lect(1)
        mlree2 = lect(2)
        mlree6 = lect(6)
        mlree3 = lect(3)
        if (mlree1.le.0 .or. mlree2.le.0 .or. mlree6.le.0 .or.
     &      mlree3.le.0) then
          call erreur(201)
          GOTO 900
        endif
        segact,mlree1,mlree2,mlree6
        n1 = mlree1.prog(/1)
        n2 = mlree2.prog(/1)
        n6 = mlree6.prog(/1)
        if (n1.LE.0 .or. n1.NE.n2 .or. n1.NE.n6 .or. n2.NE.n6) then
          segdes,mlree1,mlree2,mlree6
          call ERREUR(217)
          GOTO 900
        endif
        segact,mlree3
        ma = mlree3.prog(/1)
        segdes,mlree3
        if (ma.LE.0) then
          segdes,mlree1,mlree2,mlree6
          call ERREUR(21)
          GOTO 900
        endif
        segini,mlre11=mlree3
        lect(3) = mlre11
        jg = ma
        segini mlre12
        lect(11) = mlre12
        jg = ma*ma
        segini mlree8
        lect(8) = mlree8
        jg = ma
        segini mlree9
        lect(9) = mlree9
        jg = 2
        segini mlreel
        prog(1) = -1.
        prog(2) = -1.
        lect(10) = mlreel
c
* dimensionner critere
        ip2 = mlree2
        call ecrobj('LISTREEL',ip2)
        call opobje(14)
        call maximu(1)
        call lirree(umordo,1,iretou)
        if (umordo.GT.XSPETI.and.umordo.lt.1.d0) then 
            xmordo(lapro) = umordo * umordo
        else
            xmordo(lapro) = 1.D0
        endif 
c
      else
c*    else if (iprocc.gt.0) then
        do idon = 1, NTYP
          if (lect(idon).LE.0) then
*            write(ioimp,*) '!! IPROCC =',iprocc
*            iprocc = 0
            iprocm(lapro,1) = 0
            iprocm(lapro,2) = 0
            if (idon.eq.7) then
              moterr(1:8) ='PROCEDUR'
              call erreur(37)
            else
              call erreur(201)
            endif
            GOTO 900
          endif
        enddo
 
        call levma1(mlenti)
      endif
 
      mlreel = lect(10)
      chi2 = prog(1)
      alamda = prog(2)
*
*      write(ioimp,*)'levm',xmordo(lapro),chi2,(abs(chi2)/xmordo(lapro))
      if(((abs(chi2)/xmordo(lapro)).lt.1.e-6.and.chi2.gt.0)
     &.or.iprocc.gt.100.or.alamda.gt.1.e11) goto 80
*
      iprocc = iprocc + 1
      iprocm(lapro,2) = iprocc
      ip1 = lect(11)
      call nomobj('LISTREEL','WLEVPARO',ip1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('PARO')
      ip1 = lect(9)
      call nomobj('LISTREEL','WLEVBETA',ip1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('BETA')
      ip1 = lect(8)
      call nomobj('LISTREEL','WLEVALPH',ip1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('ALPH')
      ip1 = lect(7)
      call ecrobj('PROCEDUR',ip1)
      call ECRCHA('PROC')
      ip1 = lect(10)
      call nomobj('LISTREEL','WLEVCHI2',ip1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('CHI2')
      ip1 = lect(6)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('SIGM')
      ip1 = lect(3)
      call nomobj('LISTREEL','WLEVPARA',ip1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('PARA')
      ip1 = lect(2)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('ORDO')
      ip1 = lect(1)
      call ecrobj('LISTREEL',ip1)
      call ECRCHA('ABSC')
      call ECRCHA('LEVM')
      ip1 = lect(3)
      call ecrobj('LISTREEL',ip1)
      ip1 = lect(1)
      call ecrobj('LISTREEL',ip1)
      ip1 = lect(7)
      call ecrobj('PROCEDUR',ip1)
      GOTO 900
 
 80   continue
      if (alamda.gt.1.e11) write(ioimp,*) 'je ne sais pas faire mieux'
* oubli des objets intermediares
      ip1 = lect(8)
      CALL NOMOBJ('ANNULE  ','WLEVALPH',ip1)
      ip1 = lect(9)
      CALL NOMOBJ('ANNULE  ','WLEVBETA',ip1)
      ip1 = lect(11)
      CALL NOMOBJ('ANNULE  ','WLEVPARO',ip1)
 
      mlree1 = lect(1)
      mlree2 = lect(2)
      mlree3 = lect(6)
      segdes mlree1,mlree2,mlree3
      mlree1 = lect(3)
      segini,mlree2=mlree1
      ip1 = mlree1
      CALL NOMOBJ('ANNULE  ','WLEVPARA',ip1)
      ipout = mlree2
      mlreel = lect(10)
      chi2 = prog(1)
      ip1 = mlreel
      CALL NOMOBJ('ANNULE  ','WLEVPARA',ip1)
      call ECRREE(chi2)
      call ECROBJ('LISTREEL', ipout)
      iprocm(lapro,1) = 0
      iprocm(lapro,2) = 0
 
 900  CONTINUE
      segsup mlenti
 
      RETURN
      END