equil.f 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750
  1. PROGRAM EQUIL
  2. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  3. PARAMETER (LDA=1000)
  4. DIMENSION A(LDA,LDA)
  5. DIMENSION B(LDA)
  6. DIMENSION IPVT(LDA)
  7. DIMENSION Z(LDA), S(LDA), PRF(LDA), PRFTOT(LDA), TPRF(LDA)
  8. DIMENSION X(LDA), IFIRM(LDA)
  9. C CHARACTER ( LEN = 9 ) STR
  10. DIMENSION XP(LDA), VP(LDA), XHAT(LDA), XM(LDA), XM1(LDA)
  11. DIMENSION P(LDA), DX(LDA), VIP(LDA), VIM(LDA)
  12. DIMENSION PP(LDA), PM(LDA), ZP(LDA), ZM(LDA)
  13. DIMENSION HM(LDA), HP(LDA), GRD(LDA), PRF0(LDA)
  14. DIMENSION TCM(LDA),TCP(LDA),TC(LDA), TCTOT(LDA)
  15. CHARACTER (LEN=5) STRNM
  16. CHARACTER (LEN=5) STRNMP
  17. CHARACTER (LEN=5) STRN
  18. COMMON IDEBUG
  19. IDEBUG = 0
  20. OPEN (UNIT = 1, FILE = "cases/conf33-s.dat" )
  21. READ (1, *) N, M, C, D
  22. NM = N * M
  23. PRINT*, N, M, C, D, NM
  24. WRITE(STRNM,'(I5)') NM
  25. WRITE(STRNMP,'(I5)') NM+1
  26. WRITE(STRN,'(I5)') N
  27. C DO 2, I=1, NM
  28. C X(I) = D * DBLE(I-1)/DBLE(NM)
  29. C2 CONTINUE
  30. 1 CONTINUE
  31. READ (1, *, END=200) IXX, Y, (IFIRM(J), J=1,NM), IX
  32. PRINT "(A7,"//STRNM//"I9)", "IFIRM", (IFIRM(J), J = 1, NM )
  33. CALL OPTLOCATION(NM, IFIRM, M, D, C, XM, XP, TGRD, IT)
  34. CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
  35. CALL XPOS(NM, XP, X)
  36. CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P)
  37. CALL PPM (NM, P, D, PP, PM)
  38. CALL ZPM (NM, P, C, D, PP, PM, XP, XM, X, ZP, ZM)
  39. CALL SHARE ( NM, C, PP, PM, XP, XM, S )
  40. CALL PROFIT (NM, P, S, PRF)
  41. CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
  42. CALL XPOS(NM, XP, X)
  43. CALL TRSPCOST(NM,IFIRM,C,ZM,ZP,TCM,TCP,TC,TCTOT,TCAGG,TCMAX)
  44. AGGPRF = 0.0D0
  45. DO 10, I = 1, N
  46. 10 AGGPRF = AGGPRF + PRFTOT(I)
  47. PRINT "(A7,"//STRNM//"F9.5)", "XP", (XP(J), J=1,NM)
  48. PRINT "(A7,"//STRNM//"F9.5)", "X" , (X(J), J=1,NM)
  49. PRINT "(A7,"//STRNM//"F9.5)", "ZM", (ZM(J), J=1,NM)
  50. PRINT "(A7,"//STRNM//"F9.5)", "ZP", (ZP(J), J=1,NM)
  51. PRINT "(A7,"//STRNM//"F9.5)", "S ", (S(J), J=1,NM)
  52. PRINT "(A7,"//STRNM//"F9.5)", "P ", (P(J), J=1,NM)
  53. PRINT "(A7,"//STRNM//"F9.5)", "PRF",(PRF(J),J=1,NM)
  54. PRINT "(A7,"//STRNM//"F9.5)", "PRFTOT",(PRFTOT(J),J=1,N)
  55. PRINT "(A7,"//STRNM//"F9.3)", "TCM",(TCM(J),J=1,NM)
  56. PRINT "(A7,"//STRNM//"F9.3)", "TCP",(TCP(J),J=1,NM)
  57. PRINT "(A7,"//STRNM//"F9.3)", "TC" ,(TC (J),J=1,NM)
  58. PRINT "(A7,"//STRNM//"F9.3)", "TCTOT" ,(TCTOT(J),J=1,N)
  59. PRINT '(I12,I12,E10.1,3F9.3)', IX, IT, TGRD, TCAGG, TCMAX, AGGPRF
  60. GOTO 1
  61. 200 CONTINUE
  62. CLOSE(1)
  63. END
  64. SUBROUTINE TRSPCOST(NM,IFIRM,C,ZM,ZP,TCM,TCP,TC,TCTOT,TCAGG,TCMAX)
  65. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  66. PARAMETER (LDA=1000)
  67. DIMENSION IFIRM(LDA), ZM(LDA), ZP(LDA)
  68. DIMENSION TCM(LDA),TCP(LDA),TC(LDA), TCTOT(LDA)
  69. TCMAX = 0.0D0
  70. TC = 0.0D0
  71. TCTOT = 0.0D0
  72. TCAGG = 0.0D0
  73. DO 1, I = 1, NM
  74. TCM ( I ) = C * ZM ( I ) ** 3 / 3.0D0
  75. TCP ( I ) = C * ZP ( I ) ** 3 / 3.0D0
  76. TC ( I ) = TC ( I ) + TCM(I) + TCP(I)
  77. TCTOT(IFIRM(I)) = TCTOT(IFIRM(I)) + TC(I)
  78. TCAGG = TCAGG + TC(I)
  79. IF (TCM(I).GT.TCMAX) TCMAX = TCM(I)
  80. IF (TCP(I).GT.TCMAX) TCMAX = TCP(I)
  81. 1 CONTINUE
  82. RETURN
  83. END
  84. FUNCTION Z(K, NM, C, D)
  85. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  86. DIMENSION X(NM), P(NM)
  87. TWO = 2.D0
  88. IF (K.LT.NM) THEN
  89. Z = (X(K)+X(K+1))/TWO + (P(K+1)-P(K))/(TWO*C*(X(K+1)-X(K)))
  90. ELSE
  91. Z = (X(K)+ D + X(1))/TWO + (P(1)-P(K))/(TWO*C*(D+X(1)-X(K)))
  92. ENDIF
  93. C Z(K) = X(K) + XP(K)/TWO + PP(K) / ( TWO * C * XP(K) )
  94. RETURN
  95. END
  96. SUBROUTINE EQPRICE (NM , XP, XM, VIP, VIM, C, D, P)
  97. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  98. PARAMETER (LDA=1000)
  99. DIMENSION A(LDA,LDA), B(LDA)
  100. DIMENSION IPVT(LDA), Z(LDA)
  101. DIMENSION XP(*), XM(*), VIP(*), VIM(*)
  102. DIMENSION P(NM)
  103. ZERO = 0.D0
  104. ONE = 1.D0
  105. TWO = 2.D0
  106. A = ZERO
  107. DO 10, I = 1, NM
  108. A(I,I) = TWO * (XP(I) + XM(I)) / ( XP(I) * XM(I) )
  109. B(I) = C * (XP(I) + XM(I))
  110. 10 CONTINUE
  111. DO 11, I=1,NM-1
  112. A(I ,I+1) = - ( VIP(I) + ONE ) / XP(I)
  113. A(I+1,I ) = - ( VIP(I) + ONE ) / XP(I)
  114. 11 CONTINUE
  115. A(NM,1 ) = - ( VIP(NM) + ONE ) / XP(NM)
  116. A(1 ,NM) = - ( VIP(NM) + ONE ) / XP(NM)
  117. C DO 100, I=1, NM
  118. C100 PRINT '(13F7.2)', (A(I,J), J=1,9), B(I)
  119. JOB = 0
  120. CALL DGECO (A, LDA, NM, IPVT, RCOND, Z)
  121. CALL DGESL (A, LDA, NM, IPVT, B, JOB)
  122. P = B
  123. RETURN
  124. END
  125. SUBROUTINE XPM ( N, X, D , XP, XM )
  126. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  127. DIMENSION XP(*), XM(*), X(*)
  128. DO 1, I = 1, N - 1
  129. IF(X(I).GT.X(I+1)) THEN
  130. XP(I) = D - X(I) + X(I+1)
  131. ELSE
  132. XP(I) = X(I+1) - X(I)
  133. ENDIF
  134. 1 CONTINUE
  135. IF(X(N).GT.X(1)) THEN
  136. XP(N) = X(1) + D - X(N)
  137. ELSE
  138. XP(N) = X(1) - X(N)
  139. ENDIF
  140. XM ( 1 ) = XP ( N )
  141. DO 2, I = 2, N
  142. 2 XM ( I ) = XP ( I - 1 )
  143. RETURN
  144. END
  145. SUBROUTINE PPM ( N, P, D , PP, PM )
  146. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  147. DIMENSION PP(*), PM(*), P(*)
  148. DO 10, I = 1, N - 1
  149. 10 PP ( I ) = P ( I+1 ) - P ( I )
  150. PP ( N ) = P ( 1 ) - P ( N )
  151. PM ( 1 ) = PP ( N )
  152. DO 20, I = 2, N
  153. 20 PM ( I ) = PP ( I - 1 )
  154. RETURN
  155. END
  156. SUBROUTINE ZPM ( N, P, C, D , PP, PM, XP, XM, X, ZP, ZM )
  157. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  158. DIMENSION ZP(*), ZM(*), XP(*), XM(*), PP(*), PM(*), P(*), X(*)
  159. ZERO = 0.D0
  160. HALF = 0.5D0
  161. TWO = 2.D0
  162. DO 10, I = 1, N
  163. ZP ( I ) = HALF*XP(I)+PP(I)/(TWO*C*XP(I))
  164. C IF (ZP(I).GT.D) ZP(I) = ZP(I) - D
  165. ZM ( I ) = HALF*XM(I)-PM(I)/(TWO*C*XM(I))
  166. C IF (ZM(I).LT.ZERO) ZM(I) = ZM(I) + D
  167. 10 CONTINUE
  168. RETURN
  169. END
  170. SUBROUTINE SHARE ( N, C , PP, PM, XP, XM, S )
  171. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  172. DIMENSION XP(*), XM(*), PP(*), PM(*), S(*)
  173. TWO = 2.D0
  174. DO 1, I = 1, N
  175. 1 S(I) = (XP(I)+XM(I))/TWO + (PP(I)/XP(I)-PM(I)/XM(I))/(TWO*C)
  176. RETURN
  177. END
  178. SUBROUTINE PROFIT ( N, P, S, PRF )
  179. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  180. DIMENSION PRF(*), P(*), S(*)
  181. DO 10, I = 1, N
  182. 10 PRF(I) = P(I) * S(I)
  183. RETURN
  184. END
  185. SUBROUTINE NEIGHBOUR ( N, IFIRM, VIP, VIM )
  186. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  187. DIMENSION IFIRM(*), VIP(*), VIM(*)
  188. DO 1, I = 1, N
  189. 1 VIP(I) = 0.D0
  190. DO 2, I = 1, N-1
  191. 2 IF(IFIRM(I).EQ.IFIRM(I+1)) VIP(I) = 1.D0
  192. IF(IFIRM(N).EQ.IFIRM(1 )) VIP(N) = 1.D0
  193. VIM(1) = VIP(N)
  194. DO 3, I = 2, N
  195. 3 VIM(I) = VIP(I-1)
  196. RETURN
  197. END
  198. SUBROUTINE TOTPROFIT (NM, IFIRM, PRF, PRFTOT)
  199. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  200. DIMENSION PRF(*), PRFTOT(*), IFIRM (*)
  201. ZERO = 0.D0
  202. DO 1, I = 1, NM
  203. 1 PRFTOT(I) = ZERO
  204. DO 2, I = 1, NM
  205. 2 PRFTOT ( IFIRM(I) ) = PRFTOT ( IFIRM (I) ) + PRF(I)
  206. RETURN
  207. END
  208. SUBROUTINE VECH (N, I, EPS, HM, HP)
  209. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  210. DIMENSION HM(*), HP(*)
  211. DO 1, J = 1, N
  212. HM(J) = 0.D0
  213. HP(J) = 0.D0
  214. 1 CONTINUE
  215. IF(I.EQ.1) THEN
  216. HM(1) = -EPS
  217. HP(1) = EPS
  218. HM(2) = EPS
  219. HP(N) = -EPS
  220. ELSEIF(I.EQ.N) THEN
  221. HM(N) = -EPS
  222. HP(N) = EPS
  223. HM(1) = EPS
  224. HP(N-1) = -EPS
  225. ELSE
  226. HM(I) = -EPS
  227. HP(I) = EPS
  228. HM(I+1) = EPS
  229. HP(I-1) = -EPS
  230. ENDIF
  231. RETURN
  232. END
  233. SUBROUTINE XPROFIT (NM, IFIRM, XM, XP, VIM, VIP, C, D, PRF, TPRF)
  234. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  235. PARAMETER (LDA=1000)
  236. DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
  237. DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), TPRF(LDA), S(LDA)
  238. CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P)
  239. CALL PPM (NM, P, D, PP, PM)
  240. CALL SHARE ( NM, C, PP, PM, XP, XM, S )
  241. CALL PROFIT (NM, P, S, PRF)
  242. CALL TOTPROFIT (NM, IFIRM, PRF, TPRF )
  243. RETURN
  244. END
  245. FUNCTION QUADVAR (N, VECT1, VECT2)
  246. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  247. PARAMETER ( LDA = 1000 )
  248. DIMENSION VECT1 (LDA), VECT2(LDA)
  249. X = 0.0D0
  250. DO 10, I = 1, N
  251. 10 X = X + ( VECT1(I) - VECT2(I) ) ** 2
  252. QUADVAR = SQRT(X)
  253. RETURN
  254. END
  255. FUNCTION VECNORM (N, VEC)
  256. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  257. PARAMETER ( LDA = 1000 )
  258. DIMENSION VEC (LDA)
  259. X = 0.0D0
  260. DO 10, I = 1, N
  261. 10 X = X + VEC(I) ** 2
  262. VECNORM = SQRT(X)
  263. RETURN
  264. END
  265. FUNCTION ISUCC (N, I)
  266. IF (I.LT.N) ISUCC = I + 1
  267. IF (I.EQ.N) ISUCC = 1
  268. RETURN
  269. END
  270. FUNCTION IPRED (N, I)
  271. IF (I.EQ.1) IPRED = N
  272. IF (I.GT.1) IPRED = I - 1
  273. RETURN
  274. END
  275. FUNCTION ICOUNT (N, K, IFIRM, I)
  276. C Returns the last adjacent outlet belonging to
  277. C firm K; if outlet I does not belong to firm K
  278. C it returns zero.
  279. PARAMETER (LDA=1000)
  280. DIMENSION IFIRM(LDA)
  281. II = 0
  282. IF( IFIRM(I).NE.K.OR.
  283. & (IFIRM(I).EQ.K.AND.IFIRM(IPRED(N,I)).EQ.K) ) THEN
  284. ICOUNT = II
  285. ELSE
  286. II = I
  287. 1 CONTINUE
  288. IF( IFIRM(ISUCC(N,II)).EQ.K ) THEN
  289. II = ISUCC(N,II)
  290. GOTO 1
  291. ENDIF
  292. ICOUNT = II
  293. ENDIF
  294. RETURN
  295. END
  296. SUBROUTINE MAXPRF1 ( NM, K, L, IFIRM, XM, XP, C, D)
  297. C Returns the optimum positions XM() and XP()
  298. C for firm L
  299. C INPUT:
  300. C NM -> total number of outlets
  301. C K -> number of outlets for each firm
  302. C L -> the firm for which profit is maximized
  303. C IFIRM (NM) -> IFIRM (I) is the index of the firm
  304. C for which outlet I belongs
  305. C INPUT/OUTPUT:
  306. C XM(NM) -> X(I) = X(I) - X(I-1)
  307. C XP(NM) -> X(I) = X(I+1) - X(I)
  308. C
  309. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  310. PARAMETER (LDA=1000)
  311. DIMENSION IFIRM (LDA), XM (LDA), XP (LDA)
  312. DIMENSION XVAR(LDA)
  313. DIMENSION YM(LDA), YP(LDA)
  314. DIMENSION A(LDA,LDA), B(LDA), IVAR(LDA)
  315. DIMENSION AA(LDA, LDA)
  316. DIMENSION IFFIRM(LDA)
  317. DIMENSION W(1000000)
  318. COMMON /VARS1/ NTOT, IX, IL, IVAR, IFFIRM
  319. COMMON /VARS2/ YM, YP
  320. COMMON /PARAM1/ CST, DIST
  321. CST = C
  322. DIST = D
  323. IX = NM
  324. IL = L
  325. ZERO = 0.0D0
  326. ONE = 1.0D0
  327. A = ZERO
  328. B = ZER0
  329. XMIN = 0.01D0
  330. C DO 10, I = 1, NM
  331. C10 PRINT* , I, IFIRM(I), ICOUNT(NM, L, IFIRM, I)
  332. I = 1
  333. IC = 0
  334. ICONS = 0
  335. IOUTLET = 0
  336. IDX = 0
  337. IVAR = 0
  338. 1 CONTINUE
  339. IC = ICOUNT( NM, L, IFIRM, I )
  340. IF (IC.EQ.0) THEN
  341. I = ISUCC(NM,I)
  342. ELSE IF (IC.GE.I) THEN
  343. IOUTLET = IOUTLET + ( IC - I + 1 )
  344. ICONS = ICONS + 1
  345. IDX = IDX + 1
  346. IVAR ( IDX ) = IPRED(NM, I )
  347. A ( IPRED(NM,I), ICONS ) = - ONE
  348. B ( ICONS ) = ZERO - XMIN
  349. DO 5, IT = I, IC
  350. ICONS = ICONS + 1
  351. IDX = IDX + 1
  352. IVAR ( IDX ) = IT
  353. A ( IT , ICONS ) = - ONE
  354. B ( ICONS ) = ZERO - XMIN
  355. 5 CONTINUE
  356. XSUM = ZERO
  357. ICONS = ICONS + 1
  358. A (IPRED(NM,I), ICONS) = ONE
  359. XSUM = XP (IPRED(NM,I) )
  360. DO 6, IT = I, IC
  361. XSUM = XSUM + XP( IT )
  362. A (IT, ICONS) = ONE
  363. 6 CONTINUE
  364. B(ICONS) = XSUM
  365. XSUM = ZERO
  366. ICONS = ICONS + 1
  367. A (IPRED(NM,I), ICONS) = - ONE
  368. XSUM = XP (IPRED(NM,I) )
  369. DO 7, IT = I, IC
  370. XSUM = XSUM + XP( IT )
  371. A (IT, ICONS) = - ONE
  372. 7 CONTINUE
  373. B(ICONS) = - XSUM
  374. I = ISUCC(NM, IC)
  375. ELSE
  376. IOUTLET = IOUTLET + ( NM - I + 1 ) + IC
  377. ICONS = ICONS + 1
  378. IDX = IDX + 1
  379. IVAR ( IDX ) = IPRED( NM, I )
  380. A ( IPRED(NM,I), ICONS ) = - ONE
  381. DO 11, IT = I, NM
  382. ICONS = ICONS + 1
  383. IDX = IDX + 1
  384. IVAR ( IDX ) = IT
  385. A ( IT , ICONS ) = - ONE
  386. B( ICONS ) = ZERO - XMIN
  387. 11 CONTINUE
  388. DO 12, IT = 1, IC
  389. ICONS = ICONS + 1
  390. IDX = IDX + 1
  391. IVAR ( IDX ) = IT
  392. A ( IT , ICONS ) = - ONE
  393. B( ICONS ) = ZERO - XMIN
  394. 12 CONTINUE
  395. XSUM = ZERO
  396. ICONS = ICONS + 1
  397. A (IPRED(NM,I), ICONS) = ONE
  398. XSUM = XP (IPRED(NM,I) )
  399. DO 13, IT = I, NM
  400. XSUM = XSUM + XP ( IT )
  401. 13 A (IT, ICONS) = ONE
  402. DO 14, IT = 1, IC
  403. XSUM = XSUM + XP ( IT )
  404. 14 A (IT, ICONS) = ONE
  405. B(ICONS) = XSUM
  406. XSUM = ZERO
  407. ICONS = ICONS + 1
  408. A (IPRED(NM,I), ICONS) = -ONE
  409. XSUM = XP (IPRED(NM,I) )
  410. DO 15, IT = I, NM
  411. XSUM = XSUM + XP ( IT )
  412. 15 A (IT, ICONS) = -ONE
  413. DO 16, IT = 1, IC
  414. XSUM = XSUM + XP ( IT )
  415. 16 A (IT, ICONS) = -ONE
  416. B(ICONS) = -XSUM
  417. I = ISUCC(NM, IC)
  418. ENDIF
  419. IF (IOUTLET.LT.K) GOTO 1
  420. DO 90, IT = 1, NM
  421. YM ( IT ) = XM ( IT )
  422. YP ( IT ) = XP ( IT )
  423. IFFIRM( IT ) = IFIRM ( IT )
  424. 90 CONTINUE
  425. C PRINT*, '----'
  426. DO 60, I1 = 1, IDX
  427. DO 60, I2 = 1, ICONS
  428. 60 AA(I1,I2) = A(IVAR(I1),I2)
  429. C PRINT *
  430. C PRINT *, ICONS
  431. C DO 50, IT = 1, NM
  432. C50 PRINT '(12F9.3,F9.3)', (A(IT,J) , J=1,ICONS)
  433. C PRINT *, '===='
  434. C PRINT '(12F6.1,F6.1)', (XP( J) , J=1,IDX)
  435. C DO 70, IT = 1, IDX
  436. C 70 PRINT '(12F6.1,F6.1)', (AA(IT, J) , J=1,ICONS)
  437. C PRINT*, '----'
  438. C PRINT '(12F6.1,F6.1)', (B( J) , J=1,ICONS)
  439. C PRINT *, '===='
  440. C PRINT*, '----'
  441. C PRINT*, IDX
  442. C PRINT *, (IVAR(J), J = 1, IDX)
  443. DO 120, IT = 1, IDX
  444. 120 XVAR(IT) = XP(IVAR(IT))
  445. NPT = IDX + 3
  446. RBEG = 1.0D0
  447. REND = 1.0D-6
  448. IPRINT = 0
  449. MAXFUN = 10000
  450. CALL LINCOA(IDX,NPT,ICONS,AA,LDA,B,XVAR,RBEG,REND,IPRINT,MAXFUN,W)
  451. C PRINT*, IDX, ICONS
  452. C PRINT*, ( IFIRM(J), J = 1, NM )
  453. C PRINT*, ( XVAR(J) , J = 1, IDX )
  454. DO 130, IT = 1, IDX
  455. 130 XP ( IVAR ( IT ) ) = XVAR ( IT )
  456. XM ( 1 ) = XP ( NM )
  457. DO 140, IT = 2, NM
  458. 140 XM ( IT ) = XP ( IT - 1 )
  459. C CALL CALFUN(IDX, XVAR, F)
  460. C PRINT*, F
  461. RETURN
  462. END
  463. SUBROUTINE CALFUN(N, XVAR, F)
  464. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  465. PARAMETER (LDA=1000)
  466. DIMENSION YM(LDA), YP(LDA)
  467. DIMENSION IVAR(LDA)
  468. DIMENSION XVAR(LDA)
  469. DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
  470. DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), TPRF(LDA), S(LDA)
  471. DIMENSION IFFIRM(LDA)
  472. COMMON /VARS1/ NTOT, IX, IL, IVAR, IFFIRM
  473. COMMON /VARS2/ YM, YP
  474. COMMON /PARAM1/ CST, DIST
  475. C = CST
  476. D = DIST
  477. NM = IX
  478. DO 1, I = 1, N
  479. 1 YP( IVAR(I) ) = XVAR(I)
  480. YM(1) = YP(NM)
  481. DO 2, I = 2, NM
  482. 2 YM(I) = YP(I-1)
  483. DO 3, I = 1, NM
  484. XM(I) = YM(I)
  485. XP(I) = YP(I)
  486. IFIRM(I) = IFFIRM(I)
  487. 3 CONTINUE
  488. CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
  489. CALL EQPRICE(NM, XP, XM, VIP, VIM, C, D, P)
  490. CALL PPM (NM, P, D, PP, PM)
  491. CALL SHARE (NM, C, PP, PM, XP, XM, S )
  492. CALL PROFIT (NM, P, S, PRF)
  493. CALL TOTPROFIT (NM, IFIRM, PRF, TPRF )
  494. F = -TPRF(IL)
  495. RETURN
  496. END
  497. SUBROUTINE NGRAD(NM, IFIRM, XM, XP, VIM, VIP, C, D, EPS, GRD)
  498. IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  499. PARAMETER (LDA=1000)
  500. DIMENSION IFIRM(LDA), XM(LDA), XP(LDA), VIM(LDA), VIP(LDA)
  501. DIMENSION GRD(LDA)
  502. DIMENSION P(LDA), PM(LDA), PP(LDA), PRF(LDA), PRFTOT(LDA), S(LDA)
  503. DIMENSION HM(LDA), HP(LDA)
  504. TWO = 2.D0
  505. DO 1, I = 1, NM
  506. CALL VECH(NM, I, EPS, HM,HP)
  507. CALL EQPRICE(NM, XP+HP, XM+HM, VIP, VIM, C, D, P)
  508. CALL PPM (NM, P, D, PP, PM)
  509. CALL SHARE ( NM, C, PP, PM, XP+HP, XM+HM, S )
  510. CALL PROFIT (NM, P, S, PRF)
  511. CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
  512. PRF1 = PRFTOT(IFIRM(I))
  513. CALL EQPRICE(NM, XP-HP, XM-HM, VIP, VIM, C, D, P)
  514. CALL PPM (NM, P, D, PP, PM)
  515. CALL SHARE ( NM, C, PP, PM, XP-HP, XM-HM, S )
  516. CALL PROFIT (NM, P, S, PRF)
  517. CALL TOTPROFIT (NM, IFIRM, PRF, PRFTOT )
  518. PRF2 = PRFTOT(IFIRM(I))
  519. GRD(I) = (PRF2-PRF1)/(TWO*EPS)
  520. 1 CONTINUE
  521. RETURN
  522. END
  523. SUBROUTINE OPTLOCATION (NM, IFIRM, M, D, C, XM, XP, TGRD, IT)
  524. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  525. C Finds locations XP that maximizes the profit of firm L
  526. PARAMETER (LDA=1000)
  527. DIMENSION XM(LDA), XP(LDA)
  528. DIMENSION P(LDA), X(LDA), XP0(LDA), GRD(LDA)
  529. DIMENSION IFIRM(LDA), VIP(LDA), VIM(LDA)
  530. CHARACTER (LEN=5) STR
  531. COMMON IDEBUG
  532. EPS = 1.0D-4
  533. ITMAX = 100
  534. WRITE(STR,'(I5)') NM
  535. C Initial solution
  536. DO 1, I=1, NM
  537. X(I) = D * DBLE(I-1)/DBLE(NM)
  538. 1 CONTINUE
  539. C Preliminary computation
  540. CALL NEIGHBOUR (NM, IFIRM, VIP, VIM)
  541. CALL XPM (NM, X, D, XP, XM)
  542. C Main loop
  543. IT = 0
  544. 10 CONTINUE
  545. IT = IT + 1
  546. XP0 = XP
  547. DO 15, L = 1, NM/M
  548. IF (IDEBUG.EQ.1) THEN
  549. PRINT*, "================================================"
  550. PRINT *, "Firm : ", L
  551. PRINT *, "Initial values"
  552. PRINT "("//STR//"F9.5)", XP(1:NM)
  553. ENDIF
  554. CALL MAXPRF1 ( NM, M, L, IFIRM, XM, XP , C, D)
  555. IF (IDEBUG.EQ.1) THEN
  556. PRINT *, "New values"
  557. PRINT "("//STR//"F9.5)", XP(1:NM)
  558. PRINT*, "================================================"
  559. ENDIF
  560. 15 CONTINUE
  561. XDIFF = QUADVAR(NM, XP0, XP)
  562. IF (XDIFF.GT.EPS.AND.IT.LT.ITMAX) GOTO 10
  563. CALL NGRAD(NM, IFIRM, XM, XP, VIM, VIP, C, D, EPS, GRD)
  564. TGRD = VECNORM ( NM, GRD )
  565. RETURN
  566. END
  567. SUBROUTINE XPOS(N, XP, X)
  568. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  569. PARAMETER(LDA=1000)
  570. DIMENSION XP(LDA), X(LDA)
  571. X(1) = 0.0D0
  572. DO 1, I = 2, N
  573. 1 X(I) = X( I - 1 ) + XP( I - 1 )
  574. RETURN
  575. END