C ASYMP     SOURCE    FD218221  26/10/07    21:15:05     12662          
      SUBROUTINE ASYMP(ITER,N,M,XVAL,XOLD1,XOLD2,XMIN,XMAX,
     &                 RAA0EPS,RAAEPS,DF0DX,DFDX,
     &                 ASYINIT,ASYDECR,ASYINCR,ASYMIN,ASYMAX,
     &                 LOW,UPP,RAA0,RAA)

C ----------------------------------------------------------------------
C Ce code est une adaptation en FORTRAN/ESOPE de l'algorithme de la GC-MMA
C (Methode des Asymptotes Mobiles Globalement Convergente) propose initialement
C en langage Matlab par K. Svarberg (sources : https://www.smoptit.se/)
C ----------------------------------------------------------------------
C    Cette subroutine met a jour les asymptotes et les parametres de courbure
C    des fonctions objectif et contraintes definissant les sous problemes convexes
C    pour la GC-MMA
C
C  Entrees :
C  ---------
C  ITER      = Numero de l'iteration courante
C  N         = Nombre de variables x_j
C  M         = Nombre de contraintes
C  XVAL      = Vecteur (N) des variables x_j courantes
C  XOLD1     = XVAL, a 1 iteration  precedente  (a condition que iter>1)
C  XOLD2     = XVAL, a 2 iterations precedentes (a condition que iter>2)
C  XMIN/XMAX = Vecteurs (N) des bornes inferieures/superieures pour les variables x_j
C  RAA0EPS   = Scalaire minimum pour RAA0
C  RAAEPS    = Vecteur (M) minimum pour RAA
C  DF0DX     = Vecteur (N) des valeurs du gradient de la fonction objectif f_0,
C              par rapport aux variables x_j, calcules en XVAL
C  DFDX      = Matrice (M x N) des valeurs du gradient des fonctions contraintes f_i,
C              par rapport aux variables x_j, calcules en XVAL
C              DFDX(i,j) = gradient de f_i par rapport a x_j
C  ASYINIT   = Scalaire pour calculer la distance initiale entre les asymptotes
C  ASYDECR   = Scalaire pour decroitre de la distance entre les asymptotes
C  ASYINCR   = Scalaire pour incrementer la distance entre les asymptotes
C  ASYMIN    = Scalaire pour calculer la distance minimale entre les asymptotes
C  ASYMAX    = Scalaire pour calculer la distance maximale entre les asymptotes
C  LOW/UPP   = Vecteurs (N) des asymptotes inferieures/superieures du sous-probleme
C
C  Sorties :
C  ---------
C  LOW/UPP   = Vecteurs (N) mis a jour
C  RAA0      = Scalaire valeur approchee de la fonction objectif par le sous probleme
C  RAA       = Vecteur (M) des valeurs approchees des fonctions contraintes par le sous probleme
C ----------------------------------------------------------------------

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

C     Arguments
      INTEGER ITER,N,M
      REAL*8 XVAL(N),XOLD1(N),XOLD2(N),XMIN(N),XMAX(N),RAA0EPS,RAAEPS(M)
      REAL*8 DF0DX(N),DFDX(M,N),ASYINIT,ASYDECR,ASYINCR,ASYMIN,ASYMAX
      REAL*8 LOW(N),UPP(N)
      REAL*8 RAA0,RAA(M)

C     Constante locale
      PARAMETER (XMAMIEPS=1.0D-5)

C     Calcul de RAA0 et RAA
      RAA0 = 0.0D0
      DO 10 I = 1,N
        XMAMI = MAX(XMAX(I)-XMIN(I),XMAMIEPS)
        RAA0 = RAA0 + ABS(DF0DX(I)) * XMAMI
 10   CONTINUE
      RAA0 = MAX(RAA0EPS,(0.1D0 / DBLE(N)) * RAA0)

      DO 20 J = 1,M
        RAAJ = 0.0D0
        DO 21 I = 1,N
          XMAMI = MAX(XMAX(I)-XMIN(I),XMAMIEPS)
          RAAJ = RAAJ + ABS(DFDX(J,I)) * XMAMI
 21     CONTINUE
        RAA(J) = MAX(RAAEPS(I),(0.1D0 / DBLE(N)) * RAAJ)
 20   CONTINUE       


C     Mise a jour de LOW et UPP
      IF (ITER.LE.2) THEN
        DO 30 I = 1,N
          XMAMI = MAX(XMAX(I)-XMIN(I),XMAMIEPS)
          LOW(I) = XVAL(I) - ASYINIT * XMAMI
          UPP(I) = XVAL(I) + ASYINIT * XMAMI
 30     CONTINUE

      ELSE
        DO 40 I = 1, N
          XMAMI = MAX(XMAX(I)-XMIN(I),XMAMIEPS)
          DXX = (XVAL(I) - XOLD1(I)) * (XOLD1(I) - XOLD2(I))
          IF (DXX.GT.0.0D0) THEN
            FACTOR = ASYINCR
          ELSEIF (DXX.LT.0.0D0) THEN
            FACTOR = ASYDECR
          ELSE
            FACTOR = 1.0D0
          ENDIF
          LOW(I) = XVAL(I) - FACTOR * (XOLD1(I) - LOW(I))
          LOW(I) = MAX(LOW(I), XVAL(I) - ASYMAX * XMAMI)
          LOW(I) = MIN(LOW(I), XVAL(I) - ASYMIN * XMAMI)
          UPP(I) = XVAL(I) + FACTOR * (UPP(I) - XOLD1(I))
          UPP(I) = MIN(UPP(I), XVAL(I) + ASYMAX * XMAMI)
          UPP(I) = MAX(UPP(I), XVAL(I) + ASYMIN * XMAMI)
 40     CONTINUE
      ENDIF

      RETURN
      END