PROGRAM EQUIL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      PARAMETER (LDA=1000)
      DIMENSION A(LDA,LDA)
      DIMENSION B(LDA)
      DIMENSION IPVT(LDA)
      DIMENSION Z(LDA), S(LDA), PRF(LDA), PRFTOT(LDA), TPRF(LDA)
      DIMENSION X(LDA), IFIRM(LDA)
C     CHARACTER ( LEN = 9 ) STR
      DIMENSION XP(LDA), VP(LDA), XHAT(LDA), XM(LDA), XM1(LDA)
      DIMENSION P(LDA), DX(LDA), VIP(LDA), VIM(LDA)
      DIMENSION PP(LDA), PM(LDA), ZP(LDA), ZM(LDA)
      DIMENSION HM(LDA), HP(LDA), GRD(LDA), PRF0(LDA)
      DIMENSION TCM(LDA),TCP(LDA),TC(LDA), TCTOT(LDA)
      CHARACTER (LEN=5) STRNM
      CHARACTER (LEN=5) STRNMP
      CHARACTER (LEN=5) STRN
      COMMON IDEBUG

      IDEBUG = 0

      OPEN (UNIT = 1, FILE = "cases/conf33-s.dat" )

      READ (1, *) N, M, C, D
      NM = N * M
      PRINT*, N, M, C, D, NM

      WRITE(STRNM,'(I5)') NM
      WRITE(STRNMP,'(I5)') NM+1
      WRITE(STRN,'(I5)') N 

C     DO 2, I=1, NM
C     X(I) = D * DBLE(I-1)/DBLE(NM) 
C2    CONTINUE 
      
 1    CONTINUE

      READ (1, *, END=200) IXX, Y, (IFIRM(J), J=1,NM), IX

      PRINT "(A7,"//STRNM//"I9)", "IFIRM", (IFIRM(J), J = 1, NM )
      CALL OPTLOCATION(NM, IFIRM, M, D, C, XM, XP, TGRD, IT)
      CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
      CALL XPOS(NM, XP, X)
      CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P) 
      CALL PPM (NM, P, D, PP, PM) 
      CALL ZPM (NM, P, C, D, PP, PM, XP, XM, X, ZP, ZM)
      CALL SHARE ( NM, C, PP, PM, XP, XM, S )
      CALL PROFIT (NM, P, S, PRF)
      CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
      CALL XPOS(NM, XP, X)
      CALL TRSPCOST(NM,IFIRM,C,ZM,ZP,TCM,TCP,TC,TCTOT,TCAGG,TCMAX)
      AGGPRF = 0.0D0
      DO 10, I = 1, N
 10   AGGPRF = AGGPRF + PRFTOT(I)
      PRINT "(A7,"//STRNM//"F9.5)", "XP", (XP(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "X" , (X(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "ZM", (ZM(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "ZP", (ZP(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "S ", (S(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "P ", (P(J), J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "PRF",(PRF(J),J=1,NM)
      PRINT "(A7,"//STRNM//"F9.5)", "PRFTOT",(PRFTOT(J),J=1,N)
      PRINT "(A7,"//STRNM//"F9.3)", "TCM",(TCM(J),J=1,NM)
      PRINT "(A7,"//STRNM//"F9.3)", "TCP",(TCP(J),J=1,NM)
      PRINT "(A7,"//STRNM//"F9.3)", "TC" ,(TC (J),J=1,NM)
      PRINT "(A7,"//STRNM//"F9.3)", "TCTOT" ,(TCTOT(J),J=1,N)
      PRINT '(I12,I12,E10.1,3F9.3)', IX, IT, TGRD, TCAGG, TCMAX, AGGPRF

      GOTO 1

200   CONTINUE 

      CLOSE(1) 


      END

      SUBROUTINE TRSPCOST(NM,IFIRM,C,ZM,ZP,TCM,TCP,TC,TCTOT,TCAGG,TCMAX)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (LDA=1000)
      DIMENSION IFIRM(LDA), ZM(LDA), ZP(LDA)
      DIMENSION TCM(LDA),TCP(LDA),TC(LDA), TCTOT(LDA)

      TCMAX = 0.0D0
      TC    = 0.0D0
      TCTOT = 0.0D0
      TCAGG = 0.0D0

      DO 1, I = 1, NM
      TCM ( I ) = C * ZM ( I ) ** 3 / 3.0D0
      TCP ( I ) = C * ZP ( I ) ** 3 / 3.0D0
      TC  ( I ) = TC ( I ) + TCM(I) + TCP(I)
      TCTOT(IFIRM(I)) = TCTOT(IFIRM(I)) + TC(I)
      TCAGG = TCAGG + TC(I)
      IF (TCM(I).GT.TCMAX) TCMAX = TCM(I)
      IF (TCP(I).GT.TCMAX) TCMAX = TCP(I)
 1    CONTINUE

      RETURN
      END

      FUNCTION Z(K, NM, C, D)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(NM), P(NM)

      TWO = 2.D0

      IF (K.LT.NM) THEN
        Z = (X(K)+X(K+1))/TWO + (P(K+1)-P(K))/(TWO*C*(X(K+1)-X(K)))
      ELSE
        Z = (X(K)+ D + X(1))/TWO + (P(1)-P(K))/(TWO*C*(D+X(1)-X(K)))
      ENDIF
C     Z(K) = X(K) + XP(K)/TWO + PP(K) / ( TWO * C * XP(K) )

      RETURN
      END

      SUBROUTINE EQPRICE (NM , XP, XM, VIP, VIM, C, D, P)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      PARAMETER (LDA=1000)
      DIMENSION A(LDA,LDA), B(LDA)
      DIMENSION IPVT(LDA), Z(LDA)
      
      DIMENSION XP(*), XM(*), VIP(*), VIM(*)
      DIMENSION P(NM)

      ZERO = 0.D0
      ONE  = 1.D0
      TWO  = 2.D0

      A = ZERO

      DO 10, I = 1, NM
      A(I,I) = TWO * (XP(I) + XM(I)) / ( XP(I) * XM(I) )
      B(I)   =   C * (XP(I) + XM(I))
  10  CONTINUE
       
      DO 11, I=1,NM-1
      A(I  ,I+1) = - ( VIP(I) + ONE ) / XP(I)
      A(I+1,I  ) = - ( VIP(I) + ONE ) / XP(I)
  11  CONTINUE

      A(NM,1 ) = - ( VIP(NM) + ONE ) / XP(NM)
      A(1 ,NM) = - ( VIP(NM) + ONE ) / XP(NM) 

C     DO 100, I=1, NM
C100  PRINT '(13F7.2)', (A(I,J), J=1,9), B(I)

      JOB = 0

      CALL DGECO (A, LDA, NM, IPVT, RCOND, Z) 
      CALL DGESL (A, LDA, NM, IPVT, B, JOB)

      P = B

      RETURN
      END

      SUBROUTINE XPM ( N, X, D , XP, XM )
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION XP(*), XM(*), X(*)

      DO 1, I = 1, N - 1
      IF(X(I).GT.X(I+1)) THEN
        XP(I) = D - X(I) + X(I+1)
      ELSE
        XP(I) = X(I+1) - X(I)
      ENDIF
 1    CONTINUE
      IF(X(N).GT.X(1)) THEN
        XP(N) = X(1) + D - X(N)
      ELSE
        XP(N) = X(1) - X(N)
      ENDIF

      XM ( 1 ) = XP ( N )
      DO 2, I = 2, N 
 2    XM ( I ) = XP ( I - 1 )

      RETURN
      END

      SUBROUTINE PPM ( N, P, D , PP, PM )
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION PP(*), PM(*), P(*)

      DO 10, I = 1, N - 1
10    PP ( I ) = P ( I+1 ) - P ( I )
      PP ( N ) = P ( 1   ) - P ( N )

      PM ( 1 ) = PP ( N )
      DO 20, I = 2, N 
20    PM ( I ) = PP ( I - 1 )

      RETURN
      END

      SUBROUTINE ZPM ( N, P, C, D , PP, PM, XP, XM, X, ZP, ZM )
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION ZP(*), ZM(*), XP(*), XM(*), PP(*), PM(*), P(*), X(*)

      ZERO = 0.D0
      HALF = 0.5D0
      TWO = 2.D0

      DO 10, I = 1, N 
      ZP ( I ) = HALF*XP(I)+PP(I)/(TWO*C*XP(I)) 
C     IF (ZP(I).GT.D) ZP(I) = ZP(I) - D
      ZM ( I ) = HALF*XM(I)-PM(I)/(TWO*C*XM(I)) 
C     IF (ZM(I).LT.ZERO) ZM(I) = ZM(I) + D
10    CONTINUE

      RETURN
      END

      SUBROUTINE SHARE ( N, C , PP, PM, XP, XM, S )
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION XP(*), XM(*), PP(*), PM(*), S(*)

      TWO = 2.D0

      DO 1, I = 1, N
1     S(I) = (XP(I)+XM(I))/TWO + (PP(I)/XP(I)-PM(I)/XM(I))/(TWO*C) 

      RETURN
      END

      SUBROUTINE PROFIT ( N, P, S, PRF )
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION PRF(*), P(*), S(*)

      DO 10, I = 1, N
10    PRF(I) = P(I) * S(I)

      RETURN
      END

      SUBROUTINE NEIGHBOUR ( N, IFIRM, VIP, VIM )
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IFIRM(*), VIP(*), VIM(*)

      DO 1, I = 1, N
 1    VIP(I) = 0.D0 

      DO 2, I = 1, N-1
 2    IF(IFIRM(I).EQ.IFIRM(I+1)) VIP(I) = 1.D0
      IF(IFIRM(N).EQ.IFIRM(1  )) VIP(N) = 1.D0

      VIM(1) = VIP(N)
      DO 3, I = 2, N
 3    VIM(I) = VIP(I-1)
      
      RETURN
      END

      SUBROUTINE TOTPROFIT (NM, IFIRM, PRF, PRFTOT)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 
      DIMENSION PRF(*), PRFTOT(*), IFIRM (*)

      ZERO = 0.D0

      DO 1, I = 1, NM
 1    PRFTOT(I) = ZERO

      DO 2, I = 1, NM
 2    PRFTOT ( IFIRM(I) ) = PRFTOT ( IFIRM (I) ) + PRF(I) 

      RETURN
      END

      SUBROUTINE VECH (N, I, EPS, HM, HP)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

      DIMENSION HM(*), HP(*)

      DO 1, J = 1, N
      HM(J) = 0.D0
      HP(J) = 0.D0
 1    CONTINUE
      
      IF(I.EQ.1) THEN
        HM(1) = -EPS
        HP(1) =  EPS
        HM(2) =  EPS
        HP(N) = -EPS
      ELSEIF(I.EQ.N) THEN
        HM(N)   = -EPS
        HP(N)   =  EPS
        HM(1)   =  EPS
        HP(N-1) = -EPS
      ELSE
        HM(I)   = -EPS
        HP(I)   =  EPS
        HM(I+1) =  EPS
        HP(I-1) = -EPS
      ENDIF

      RETURN
      END
        
      SUBROUTINE XPROFIT (NM, IFIRM, XM, XP, VIM, VIP, C, D, PRF, TPRF)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

      PARAMETER (LDA=1000)
      DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
      DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), TPRF(LDA), S(LDA)

      CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P) 
      CALL PPM (NM, P, D, PP, PM)
      CALL SHARE ( NM, C, PP, PM, XP, XM, S )
      CALL PROFIT (NM, P, S, PRF)
      CALL TOTPROFIT (NM, IFIRM, PRF, TPRF ) 
 
      RETURN
      END

      FUNCTION QUADVAR (N, VECT1, VECT2)
      IMPLICIT DOUBLE PRECISION  (A-H, O-Z)

      PARAMETER ( LDA = 1000 )
      DIMENSION VECT1 (LDA), VECT2(LDA)

      X = 0.0D0

      DO 10, I = 1, N
 10   X = X + ( VECT1(I) - VECT2(I) ) ** 2 

      QUADVAR = SQRT(X)

      RETURN
      END

      FUNCTION VECNORM (N, VEC)
      IMPLICIT DOUBLE PRECISION  (A-H, O-Z)

      PARAMETER ( LDA = 1000 )
      DIMENSION VEC (LDA)

      X = 0.0D0

      DO 10, I = 1, N
 10   X = X + VEC(I) ** 2 

      VECNORM = SQRT(X)

      RETURN
      END

      FUNCTION ISUCC (N, I)

      IF (I.LT.N) ISUCC = I + 1
      IF (I.EQ.N) ISUCC = 1 

      RETURN
      END

      FUNCTION IPRED (N, I) 

      IF (I.EQ.1) IPRED = N 
      IF (I.GT.1) IPRED = I - 1

      RETURN
      END

      FUNCTION ICOUNT (N, K, IFIRM, I)
C     Returns the last adjacent outlet belonging to
C     firm K; if outlet I does not belong to firm K
C     it returns zero.

      PARAMETER (LDA=1000)
      DIMENSION IFIRM(LDA)
      
      II = 0
      IF(   IFIRM(I).NE.K.OR.
     &     (IFIRM(I).EQ.K.AND.IFIRM(IPRED(N,I)).EQ.K) ) THEN

        ICOUNT = II 

      ELSE 

        II = I 
 1      CONTINUE
        IF( IFIRM(ISUCC(N,II)).EQ.K ) THEN
          II = ISUCC(N,II)
          GOTO 1
        ENDIF
        ICOUNT = II

      ENDIF

      RETURN
      END

      SUBROUTINE MAXPRF1 ( NM, K, L, IFIRM, XM, XP, C, D)
C   Returns the optimum positions XM() and XP()
C   for firm L 

C   INPUT:
C      NM -> total number of outlets
C      K  -> number of outlets for each firm
C      L  -> the firm for which profit is maximized
C      IFIRM (NM) -> IFIRM (I) is the index of the firm
C                    for which outlet I belongs

C   INPUT/OUTPUT:
C      XM(NM) -> X(I) = X(I)   - X(I-1)
C      XP(NM) -> X(I) = X(I+1) - X(I)
C      
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      
      PARAMETER (LDA=1000)
      DIMENSION IFIRM (LDA), XM (LDA), XP (LDA)
      DIMENSION XVAR(LDA)
      DIMENSION YM(LDA), YP(LDA)
      DIMENSION A(LDA,LDA), B(LDA), IVAR(LDA)
      DIMENSION AA(LDA, LDA)
      DIMENSION IFFIRM(LDA)
      DIMENSION W(1000000)

      COMMON /VARS1/ NTOT, IX, IL, IVAR, IFFIRM
      COMMON /VARS2/ YM, YP 
      COMMON /PARAM1/ CST, DIST

      CST = C
      DIST = D
      IX = NM 
      IL = L

      ZERO = 0.0D0
      ONE  = 1.0D0

      A = ZERO
      B = ZER0
      XMIN = 0.01D0

C     DO 10, I = 1, NM
C10   PRINT* , I, IFIRM(I), ICOUNT(NM, L, IFIRM, I)

      I = 1
      IC = 0
      ICONS = 0
      IOUTLET = 0
      IDX = 0 
      IVAR = 0

 1    CONTINUE

      IC = ICOUNT( NM, L, IFIRM, I )

      IF (IC.EQ.0) THEN

        I = ISUCC(NM,I)

      ELSE IF (IC.GE.I) THEN

        IOUTLET = IOUTLET + ( IC - I + 1 )

        ICONS = ICONS + 1 
        IDX = IDX + 1
        IVAR ( IDX ) = IPRED(NM, I )
        A ( IPRED(NM,I), ICONS ) = - ONE
        B ( ICONS ) = ZERO - XMIN

        DO 5, IT = I, IC 
        ICONS = ICONS + 1
        IDX = IDX + 1
        IVAR ( IDX ) = IT
        A ( IT , ICONS ) = - ONE
        B ( ICONS ) = ZERO - XMIN
  5     CONTINUE

        XSUM = ZERO
        ICONS = ICONS + 1
        A (IPRED(NM,I), ICONS) = ONE 
        XSUM = XP (IPRED(NM,I) ) 
        DO 6, IT = I, IC
        XSUM = XSUM + XP( IT )
        A (IT, ICONS) = ONE
  6     CONTINUE
        B(ICONS) = XSUM

        XSUM = ZERO
        ICONS = ICONS + 1
        A (IPRED(NM,I), ICONS) = - ONE 
        XSUM = XP (IPRED(NM,I) ) 
        DO 7, IT = I, IC
        XSUM = XSUM + XP( IT )
        A (IT, ICONS) = - ONE
  7     CONTINUE
        B(ICONS) = - XSUM 

        I = ISUCC(NM, IC)

      ELSE

        IOUTLET = IOUTLET + ( NM - I + 1 ) + IC

        ICONS = ICONS + 1 
        IDX = IDX + 1
        IVAR ( IDX ) = IPRED( NM, I )
        A ( IPRED(NM,I), ICONS ) = - ONE

        DO 11, IT = I, NM
        ICONS = ICONS + 1
        IDX = IDX + 1
        IVAR ( IDX ) =  IT 
        A ( IT , ICONS ) = - ONE
        B( ICONS ) = ZERO - XMIN
  11    CONTINUE

        DO 12, IT = 1, IC
        ICONS = ICONS + 1
        IDX = IDX + 1
        IVAR ( IDX ) =  IT 
        A ( IT , ICONS ) = - ONE
        B( ICONS ) = ZERO - XMIN
  12    CONTINUE 
        
        XSUM = ZERO
        ICONS = ICONS + 1
        A (IPRED(NM,I), ICONS) = ONE 
        XSUM = XP (IPRED(NM,I) ) 
        DO 13, IT = I, NM
        XSUM = XSUM + XP ( IT ) 
  13    A (IT, ICONS) = ONE
        DO 14, IT = 1, IC
        XSUM = XSUM + XP ( IT ) 
  14    A (IT, ICONS) = ONE
        B(ICONS) = XSUM

        XSUM = ZERO
        ICONS = ICONS + 1
        A (IPRED(NM,I), ICONS) = -ONE 
        XSUM = XP (IPRED(NM,I) ) 
        DO 15, IT = I, NM
        XSUM = XSUM + XP ( IT ) 
  15    A (IT, ICONS) = -ONE
        DO 16, IT = 1, IC
        XSUM = XSUM + XP ( IT ) 
  16    A (IT, ICONS) = -ONE
        B(ICONS) = -XSUM

        I = ISUCC(NM, IC)

      ENDIF


      IF (IOUTLET.LT.K) GOTO 1 


      DO 90, IT = 1, NM
      YM ( IT ) = XM ( IT )
      YP ( IT ) = XP ( IT )
      IFFIRM( IT ) = IFIRM ( IT )
 90   CONTINUE


C     PRINT*, '----'
      DO 60, I1 = 1, IDX
      DO 60, I2 = 1, ICONS
  60  AA(I1,I2) = A(IVAR(I1),I2)

C     PRINT *
C     PRINT *, ICONS
C     DO 50, IT = 1, NM
C50   PRINT '(12F9.3,F9.3)', (A(IT,J) , J=1,ICONS) 

C     PRINT *, '===='
C     PRINT '(12F6.1,F6.1)', (XP(   J)   , J=1,IDX) 
C     DO 70, IT = 1, IDX
C 70  PRINT '(12F6.1,F6.1)', (AA(IT,  J) , J=1,ICONS) 
C     PRINT*, '----'
C     PRINT '(12F6.1,F6.1)', (B(   J)    , J=1,ICONS) 
C     PRINT *, '===='

C     PRINT*, '----'
C     PRINT*, IDX
C     PRINT *, (IVAR(J), J = 1, IDX) 

      DO 120, IT = 1, IDX
 120  XVAR(IT) = XP(IVAR(IT))

      NPT = IDX + 3
      RBEG = 1.0D0
      REND = 1.0D-6 
      IPRINT = 0
      MAXFUN = 10000
      CALL LINCOA(IDX,NPT,ICONS,AA,LDA,B,XVAR,RBEG,REND,IPRINT,MAXFUN,W)
       
C     PRINT*, IDX, ICONS

C     PRINT*, ( IFIRM(J), J = 1, NM )
C     PRINT*, ( XVAR(J) , J = 1, IDX )
      DO 130, IT = 1, IDX
 130  XP ( IVAR ( IT ) ) = XVAR ( IT )
      XM ( 1 ) = XP ( NM )
      DO 140, IT = 2, NM
 140  XM ( IT ) = XP ( IT - 1 )
C     CALL CALFUN(IDX, XVAR, F)
C     PRINT*, F 

      RETURN
      END

      SUBROUTINE CALFUN(N, XVAR, F)
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)

      PARAMETER (LDA=1000)

      DIMENSION YM(LDA), YP(LDA)
      DIMENSION IVAR(LDA)
      DIMENSION XVAR(LDA)
      DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
      DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), TPRF(LDA), S(LDA)
      DIMENSION IFFIRM(LDA)

      COMMON /VARS1/ NTOT, IX, IL, IVAR, IFFIRM
      COMMON /VARS2/ YM, YP
      COMMON /PARAM1/ CST, DIST

      C = CST
      D = DIST
      NM = IX 

      DO 1, I = 1, N
 1    YP( IVAR(I) ) = XVAR(I)

      YM(1) = YP(NM)
      DO 2, I = 2, NM
 2    YM(I) = YP(I-1)

      DO 3, I = 1, NM
      XM(I) = YM(I)
      XP(I) = YP(I)
      IFIRM(I) = IFFIRM(I)
 3    CONTINUE 

      CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
      CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P) 
      CALL PPM (NM, P, D, PP, PM)
      CALL SHARE (NM, C, PP, PM, XP, XM, S )
      CALL PROFIT (NM, P, S, PRF)
      CALL TOTPROFIT (NM, IFIRM, PRF, TPRF ) 

      F = -TPRF(IL) 

      RETURN
      END

      SUBROUTINE NGRAD(NM, IFIRM, XM, XP, VIM, VIP, C, D, EPS, GRD)    
      IMPLICIT DOUBLE PRECISION (A-H, O-Z) 

      PARAMETER (LDA=1000)
      DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
      DIMENSION GRD(LDA)
      DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), PRFTOT(LDA), S(LDA)
      DIMENSION HM(LDA), HP(LDA)

      TWO = 2.D0

      DO 1, I = 1, NM

      CALL VECH(NM, I, EPS, HM,HP)
      CALL EQPRICE(NM, XP+HP, XM+HM, VIP, VIM, C, D, P) 
      CALL PPM (NM, P, D, PP, PM)
      CALL SHARE ( NM, C, PP, PM, XP+HP, XM+HM, S )
      CALL PROFIT (NM, P, S, PRF)
      CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
      PRF1 = PRFTOT(IFIRM(I))

      CALL EQPRICE(NM, XP-HP, XM-HM, VIP, VIM, C, D, P) 
      CALL PPM (NM, P, D, PP, PM)
      CALL SHARE ( NM, C, PP, PM, XP-HP, XM-HM, S )
      CALL PROFIT (NM, P, S, PRF)
      CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
      PRF2 = PRFTOT(IFIRM(I))

      GRD(I) = (PRF2-PRF1)/(TWO*EPS)

1     CONTINUE
 
      RETURN
      END

      SUBROUTINE OPTLOCATION (NM, IFIRM, M, D, C, XM, XP, TGRD, IT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C Finds locations XP that maximizes the profit of firm L
      PARAMETER (LDA=1000)
      DIMENSION XM(LDA), XP(LDA)
      DIMENSION P(LDA), X(LDA), XP0(LDA), GRD(LDA)
      DIMENSION IFIRM(LDA), VIP(LDA), VIM(LDA)
      CHARACTER (LEN=5) STR
      COMMON IDEBUG

      EPS = 1.0D-4
      ITMAX = 100

      WRITE(STR,'(I5)') NM
C Initial solution
      DO 1, I=1, NM
      X(I) = D * DBLE(I-1)/DBLE(NM) 
 1    CONTINUE 

C Preliminary computation
      CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
      CALL XPM (NM, X, D, XP, XM)

C Main loop
      IT = 0
 10   CONTINUE
      IT = IT + 1
      XP0 = XP 
      DO 15, L = 1, NM/M
      IF (IDEBUG.EQ.1) THEN
        PRINT*, "================================================"
        PRINT *, "Firm : ", L
        PRINT *, "Initial values"
        PRINT "("//STR//"F9.5)", XP(1:NM)
      ENDIF
      CALL MAXPRF1 ( NM, M, L, IFIRM, XM, XP , C, D)
      IF (IDEBUG.EQ.1) THEN
        PRINT *, "New values"
        PRINT "("//STR//"F9.5)", XP(1:NM)
        PRINT*, "================================================"
      ENDIF
 15   CONTINUE
      XDIFF = QUADVAR(NM, XP0, XP)
      IF (XDIFF.GT.EPS.AND.IT.LT.ITMAX) GOTO 10

      CALL NGRAD(NM, IFIRM, XM, XP, VIM, VIP, C, D, EPS, GRD)
      TGRD = VECNORM ( NM, GRD ) 

      RETURN
      END

      SUBROUTINE XPOS(N, XP, X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 

      PARAMETER(LDA=1000)
      DIMENSION XP(LDA), X(LDA)

      X(1) = 0.0D0
      DO 1, I = 2, N
 1    X(I) = X( I - 1 ) + XP( I - 1 )

      RETURN
      END