cobylb.f 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444
  1. SUBROUTINE COBYLB (N,M,MPP,X,RHOBEG,RHOEND,IPRINT,MAXFUN,
  2. 1 CON,SIM,SIMI,DATMAT,A,VSIG,VETA,SIGBAR,DX,W,IACT,ierr)
  3. IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  4. DIMENSION X(*),CON(*),SIM(N,*),SIMI(N,*),DATMAT(MPP,*),
  5. 1 A(N,*),VSIG(*),VETA(*),SIGBAR(*),DX(*),W(*),IACT(*)
  6. C
  7. C Set the initial values of some parameters. The last column of SIM holds
  8. C the optimal vertex of the current simplex, and the preceding N columns
  9. C hold the displacements from the optimal vertex to the other vertices.
  10. C Further, SIMI holds the inverse of the matrix that is contained in the
  11. C first N columns of SIM.
  12. C
  13. ierr = 0
  14. IPTEM=MIN0(N,5)
  15. IPTEMP=IPTEM+1
  16. NP=N+1
  17. MP=M+1
  18. ALPHA=0.25d0
  19. BETA=2.1d0
  20. GAMMA=0.5d0
  21. DELTA=1.1d0
  22. RHO=RHOBEG
  23. PARMU=0.0d0
  24. IF (IPRINT .GE. 2) PRINT 10, RHO
  25. 10 FORMAT (/3X,'The initial value of RHO is',1PE13.6,2X,
  26. 1 'and PARMU is set to zero.')
  27. NFVALS=0
  28. TEMP=1.0d0/RHO
  29. DO 30 I=1,N
  30. SIM(I,NP)=X(I)
  31. DO 20 J=1,N
  32. SIM(I,J)=0.0d0
  33. 20 SIMI(I,J)=0.0d0
  34. SIM(I,I)=RHO
  35. 30 SIMI(I,I)=TEMP
  36. JDROP=NP
  37. IBRNCH=0
  38. C
  39. C Make the next call of the user-supplied subroutine CALCFC. These
  40. C instructions are also used for calling CALCFC during the iterations of
  41. C the algorithm.
  42. C
  43. 40 IF (NFVALS .GE. MAXFUN .AND. NFVALS .GT. 0) THEN
  44. IF (IPRINT .GE. 1) PRINT 50
  45. 50 FORMAT (/3X,'Return from subroutine COBYLA because the ',
  46. 1 'MAXFUN limit has been reached.')
  47. ierr = 1
  48. GOTO 600
  49. END IF
  50. NFVALS=NFVALS+1
  51. CALL CALCFC (N,M,X,F,CON)
  52. RESMAX=0.0d0
  53. IF (M .GT. 0) THEN
  54. DO 60 K=1,M
  55. 60 RESMAX=DMAX1(RESMAX,-CON(K))
  56. END IF
  57. IF (NFVALS .EQ. IPRINT-1 .OR. IPRINT .EQ. 3) THEN
  58. PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
  59. 70 FORMAT (/3X,'NFVALS =',I5,3X,'F =',1PE13.6,4X,'MAXCV =',
  60. 1 1PE13.6/3X,'X =',1PE13.6,1P4E15.6)
  61. IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
  62. 80 FORMAT (1PE19.6,1P4E15.6)
  63. END IF
  64. CON(MP)=F
  65. CON(MPP)=RESMAX
  66. IF (IBRNCH .EQ. 1) GOTO 440
  67. C
  68. C Set the recently calculated function values in a column of DATMAT. This
  69. C array has a column for each vertex of the current simplex, the entries of
  70. C each column being the values of the constraint functions (if any)
  71. C followed by the objective function and the greatest constraint violation
  72. C at the vertex.
  73. C
  74. DO 90 K=1,MPP
  75. 90 DATMAT(K,JDROP)=CON(K)
  76. IF (NFVALS .GT. NP) GOTO 130
  77. C
  78. C Exchange the new vertex of the initial simplex with the optimal vertex if
  79. C necessary. Then, if the initial simplex is not complete, pick its next
  80. C vertex and calculate the function values there.
  81. C
  82. IF (JDROP .LE. N) THEN
  83. IF (DATMAT(MP,NP) .LE. F) THEN
  84. X(JDROP)=SIM(JDROP,NP)
  85. ELSE
  86. SIM(JDROP,NP)=X(JDROP)
  87. DO 100 K=1,MPP
  88. DATMAT(K,JDROP)=DATMAT(K,NP)
  89. 100 DATMAT(K,NP)=CON(K)
  90. DO 120 K=1,JDROP
  91. SIM(JDROP,K)=-RHO
  92. TEMP=0.0
  93. DO 110 I=K,JDROP
  94. 110 TEMP=TEMP-SIMI(I,K)
  95. 120 SIMI(JDROP,K)=TEMP
  96. END IF
  97. END IF
  98. IF (NFVALS .LE. N) THEN
  99. JDROP=NFVALS
  100. X(JDROP)=X(JDROP)+RHO
  101. GOTO 40
  102. END IF
  103. 130 IBRNCH=1
  104. C
  105. C Identify the optimal vertex of the current simplex.
  106. C
  107. 140 PHIMIN=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
  108. NBEST=NP
  109. DO 150 J=1,N
  110. TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
  111. IF (TEMP .LT. PHIMIN) THEN
  112. NBEST=J
  113. PHIMIN=TEMP
  114. ELSE IF (TEMP .EQ. PHIMIN .AND. PARMU .EQ. 0.0d0) THEN
  115. IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NBEST)) NBEST=J
  116. END IF
  117. 150 CONTINUE
  118. C
  119. C Switch the best vertex into pole position if it is not there already,
  120. C and also update SIM, SIMI and DATMAT.
  121. C
  122. IF (NBEST .LE. N) THEN
  123. DO 160 I=1,MPP
  124. TEMP=DATMAT(I,NP)
  125. DATMAT(I,NP)=DATMAT(I,NBEST)
  126. 160 DATMAT(I,NBEST)=TEMP
  127. DO 180 I=1,N
  128. TEMP=SIM(I,NBEST)
  129. SIM(I,NBEST)=0.0d0
  130. SIM(I,NP)=SIM(I,NP)+TEMP
  131. TEMPA=0.0d0
  132. DO 170 K=1,N
  133. SIM(I,K)=SIM(I,K)-TEMP
  134. 170 TEMPA=TEMPA-SIMI(K,I)
  135. 180 SIMI(NBEST,I)=TEMPA
  136. END IF
  137. C
  138. C Make an error return if SIGI is a poor approximation to the inverse of
  139. C the leading N by N submatrix of SIG.
  140. C
  141. ERROR=0.0d0
  142. DO 200 I=1,N
  143. DO 200 J=1,N
  144. TEMP=0.0d0
  145. IF (I .EQ. J) TEMP=TEMP-1.0d0
  146. DO 190 K=1,N
  147. 190 TEMP=TEMP+SIMI(I,K)*SIM(K,J)
  148. 200 ERROR=DMAX1(ERROR,ABS(TEMP))
  149. IF (ERROR .GT. 0.1d0) THEN
  150. IF (IPRINT .GE. 1) PRINT 210
  151. 210 FORMAT (/3X,'Return from subroutine COBYLA because ',
  152. 1 'rounding errors are becoming damaging.')
  153. ierr = 2
  154. GOTO 600
  155. END IF
  156. C
  157. C Calculate the coefficients of the linear approximations to the objective
  158. C and constraint functions, placing minus the objective function gradient
  159. C after the constraint gradients in the array A. The vector W is used for
  160. C working space.
  161. C
  162. DO 240 K=1,MP
  163. CON(K)=-DATMAT(K,NP)
  164. DO 220 J=1,N
  165. 220 W(J)=DATMAT(K,J)+CON(K)
  166. DO 240 I=1,N
  167. TEMP=0.0d0
  168. DO 230 J=1,N
  169. 230 TEMP=TEMP+W(J)*SIMI(J,I)
  170. IF (K .EQ. MP) TEMP=-TEMP
  171. 240 A(I,K)=TEMP
  172. C
  173. C Calculate the values of sigma and eta, and set IFLAG=0 if the current
  174. C simplex is not acceptable.
  175. C
  176. IFLAG=1
  177. PARSIG=ALPHA*RHO
  178. PARETA=BETA*RHO
  179. DO 260 J=1,N
  180. WSIG=0.0d0
  181. WETA=0.0d0
  182. DO 250 I=1,N
  183. WSIG=WSIG+SIMI(J,I)**2
  184. 250 WETA=WETA+SIM(I,J)**2
  185. VSIG(J)=1.0d0/SQRT(WSIG)
  186. VETA(J)=SQRT(WETA)
  187. IF (VSIG(J) .LT. PARSIG .OR. VETA(J) .GT. PARETA) IFLAG=0
  188. 260 CONTINUE
  189. C
  190. C If a new vertex is needed to improve acceptability, then decide which
  191. C vertex to drop from the simplex.
  192. C
  193. IF (IBRNCH .EQ. 1 .OR. IFLAG .EQ. 1) GOTO 370
  194. JDROP=0
  195. TEMP=PARETA
  196. DO 270 J=1,N
  197. IF (VETA(J) .GT. TEMP) THEN
  198. JDROP=J
  199. TEMP=VETA(J)
  200. END IF
  201. 270 CONTINUE
  202. IF (JDROP .EQ. 0) THEN
  203. DO 280 J=1,N
  204. IF (VSIG(J) .LT. TEMP) THEN
  205. JDROP=J
  206. TEMP=VSIG(J)
  207. END IF
  208. 280 CONTINUE
  209. END IF
  210. C
  211. C Calculate the step to the new vertex and its sign.
  212. C
  213. TEMP=GAMMA*RHO*VSIG(JDROP)
  214. DO 290 I=1,N
  215. 290 DX(I)=TEMP*SIMI(JDROP,I)
  216. CVMAXP=0.0d0
  217. CVMAXM=0.0d0
  218. DO 310 K=1,MP
  219. SUM=0.0d0
  220. DO 300 I=1,N
  221. 300 SUM=SUM+A(I,K)*DX(I)
  222. IF (K .LT. MP) THEN
  223. TEMP=DATMAT(K,NP)
  224. CVMAXP=DMAX1(CVMAXP,-SUM-TEMP)
  225. CVMAXM=DMAX1(CVMAXM,SUM-TEMP)
  226. END IF
  227. 310 CONTINUE
  228. DXSIGN=1.0d0
  229. IF (PARMU*(CVMAXP-CVMAXM) .GT. SUM+SUM) DXSIGN=-1.0d0
  230. C
  231. C Update the elements of SIM and SIMI, and set the next X.
  232. C
  233. TEMP=0.0d0
  234. DO 320 I=1,N
  235. DX(I)=DXSIGN*DX(I)
  236. SIM(I,JDROP)=DX(I)
  237. 320 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
  238. DO 330 I=1,N
  239. 330 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
  240. DO 360 J=1,N
  241. IF (J .NE. JDROP) THEN
  242. TEMP=0.0d0
  243. DO 340 I=1,N
  244. 340 TEMP=TEMP+SIMI(J,I)*DX(I)
  245. DO 350 I=1,N
  246. 350 SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
  247. END IF
  248. 360 X(J)=SIM(J,NP)+DX(J)
  249. GOTO 40
  250. C
  251. C Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO.
  252. C
  253. 370 IZ=1
  254. IZDOTA=IZ+N*N
  255. IVMC=IZDOTA+N
  256. ISDIRN=IVMC+MP
  257. IDXNEW=ISDIRN+N
  258. IVMD=IDXNEW+N
  259. CALL TRSTLP (N,M,A,CON,RHO,DX,IFULL,IACT,W(IZ),W(IZDOTA),
  260. 1 W(IVMC),W(ISDIRN),W(IDXNEW),W(IVMD))
  261. IF (IFULL .EQ. 0) THEN
  262. TEMP=0.0d0
  263. DO 380 I=1,N
  264. 380 TEMP=TEMP+DX(I)**2
  265. IF (TEMP .LT. 0.25d0*RHO*RHO) THEN
  266. IBRNCH=1
  267. GOTO 550
  268. END IF
  269. END IF
  270. C
  271. C Predict the change to F and the new maximum constraint violation if the
  272. C variables are altered from x(0) to x(0)+DX.
  273. C
  274. RESNEW=0.0d0
  275. CON(MP)=0.0d0
  276. DO 400 K=1,MP
  277. SUM=CON(K)
  278. DO 390 I=1,N
  279. 390 SUM=SUM-A(I,K)*DX(I)
  280. IF (K .LT. MP) RESNEW=DMAX1(RESNEW,SUM)
  281. 400 CONTINUE
  282. C
  283. C Increase PARMU if necessary and branch back if this change alters the
  284. C optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
  285. C reductions in the merit function and the maximum constraint violation
  286. C respectively.
  287. C
  288. BARMU=0.0d0
  289. PREREC=DATMAT(MPP,NP)-RESNEW
  290. IF (PREREC .GT. 0.0d0) BARMU=SUM/PREREC
  291. IF (PARMU .LT. 1.5d0*BARMU) THEN
  292. PARMU=2.0d0*BARMU
  293. IF (IPRINT .GE. 2) PRINT 410, PARMU
  294. 410 FORMAT (/3X,'Increase in PARMU to',1PE13.6)
  295. PHI=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
  296. DO 420 J=1,N
  297. TEMP=DATMAT(MP,J)+PARMU*DATMAT(MPP,J)
  298. IF (TEMP .LT. PHI) GOTO 140
  299. IF (TEMP .EQ. PHI .AND. PARMU .EQ. 0.0) THEN
  300. IF (DATMAT(MPP,J) .LT. DATMAT(MPP,NP)) GOTO 140
  301. END IF
  302. 420 CONTINUE
  303. END IF
  304. PREREM=PARMU*PREREC-SUM
  305. C
  306. C Calculate the constraint and objective functions at x(*). Then find the
  307. C actual reduction in the merit function.
  308. C
  309. DO 430 I=1,N
  310. 430 X(I)=SIM(I,NP)+DX(I)
  311. IBRNCH=1
  312. GOTO 40
  313. 440 VMOLD=DATMAT(MP,NP)+PARMU*DATMAT(MPP,NP)
  314. VMNEW=F+PARMU*RESMAX
  315. TRURED=VMOLD-VMNEW
  316. IF (PARMU .EQ. 0.0d0 .AND. F .EQ. DATMAT(MP,NP)) THEN
  317. PREREM=PREREC
  318. TRURED=DATMAT(MPP,NP)-RESMAX
  319. END IF
  320. C
  321. C Begin the operations that decide whether x(*) should replace one of the
  322. C vertices of the current simplex, the change being mandatory if TRURED is
  323. C positive. Firstly, JDROP is set to the index of the vertex that is to be
  324. C replaced.
  325. C
  326. RATIO=0.0d0
  327. IF (TRURED .LE. 0.0) RATIO=1.0
  328. JDROP=0
  329. DO 460 J=1,N
  330. TEMP=0.0d0
  331. DO 450 I=1,N
  332. 450 TEMP=TEMP+SIMI(J,I)*DX(I)
  333. TEMP=ABS(TEMP)
  334. IF (TEMP .GT. RATIO) THEN
  335. JDROP=J
  336. RATIO=TEMP
  337. END IF
  338. 460 SIGBAR(J)=TEMP*VSIG(J)
  339. C
  340. C Calculate the value of ell.
  341. C
  342. EDGMAX=DELTA*RHO
  343. L=0
  344. DO 480 J=1,N
  345. IF (SIGBAR(J) .GE. PARSIG .OR. SIGBAR(J) .GE. VSIG(J)) THEN
  346. TEMP=VETA(J)
  347. IF (TRURED .GT. 0.0d0) THEN
  348. TEMP=0.0d0
  349. DO 470 I=1,N
  350. 470 TEMP=TEMP+(DX(I)-SIM(I,J))**2
  351. TEMP=SQRT(TEMP)
  352. END IF
  353. IF (TEMP .GT. EDGMAX) THEN
  354. L=J
  355. EDGMAX=TEMP
  356. END IF
  357. END IF
  358. 480 CONTINUE
  359. IF (L .GT. 0) JDROP=L
  360. IF (JDROP .EQ. 0) GOTO 550
  361. C
  362. C Revise the simplex by updating the elements of SIM, SIMI and DATMAT.
  363. C
  364. TEMP=0.0d0
  365. DO 490 I=1,N
  366. SIM(I,JDROP)=DX(I)
  367. 490 TEMP=TEMP+SIMI(JDROP,I)*DX(I)
  368. DO 500 I=1,N
  369. 500 SIMI(JDROP,I)=SIMI(JDROP,I)/TEMP
  370. DO 530 J=1,N
  371. IF (J .NE. JDROP) THEN
  372. TEMP=0.0d0
  373. DO 510 I=1,N
  374. 510 TEMP=TEMP+SIMI(J,I)*DX(I)
  375. DO 520 I=1,N
  376. 520 SIMI(J,I)=SIMI(J,I)-TEMP*SIMI(JDROP,I)
  377. END IF
  378. 530 CONTINUE
  379. DO 540 K=1,MPP
  380. 540 DATMAT(K,JDROP)=CON(K)
  381. C
  382. C Branch back for further iterations with the current RHO.
  383. C
  384. IF (TRURED .GT. 0.0d0 .AND. TRURED .GE. 0.1d0*PREREM) GOTO 140
  385. 550 IF (IFLAG .EQ. 0) THEN
  386. IBRNCH=0
  387. GOTO 140
  388. END IF
  389. C
  390. C Otherwise reduce RHO if it is not at its least value and reset PARMU.
  391. C
  392. IF (RHO .GT. RHOEND) THEN
  393. RHO=0.5d0*RHO
  394. IF (RHO .LE. 1.5d0*RHOEND) RHO=RHOEND
  395. IF (PARMU .GT. 0.0d0) THEN
  396. DENOM=0.0d0
  397. DO 570 K=1,MP
  398. CMIN=DATMAT(K,NP)
  399. CMAX=CMIN
  400. DO 560 I=1,N
  401. CMIN=DMIN1(CMIN,DATMAT(K,I))
  402. 560 CMAX=DMAX1(CMAX,DATMAT(K,I))
  403. IF (K .LE. M .AND. CMIN .LT. 0.5d0*CMAX) THEN
  404. TEMP=DMAX1(CMAX,0.0d0)-CMIN
  405. IF (DENOM .LE. 0.0d0) THEN
  406. DENOM=TEMP
  407. ELSE
  408. DENOM=DMIN1(DENOM,TEMP)
  409. END IF
  410. END IF
  411. 570 CONTINUE
  412. IF (DENOM .EQ. 0.0d0) THEN
  413. PARMU=0.0d0
  414. ELSE IF (CMAX-CMIN .LT. PARMU*DENOM) THEN
  415. PARMU=(CMAX-CMIN)/DENOM
  416. END IF
  417. END IF
  418. IF (IPRINT .GE. 2) PRINT 580, RHO,PARMU
  419. 580 FORMAT (/3X,'Reduction in RHO to',1PE13.6,' and PARMU =',
  420. 1 1PE13.6)
  421. IF (IPRINT .EQ. 2) THEN
  422. PRINT 70, NFVALS,DATMAT(MP,NP),DATMAT(MPP,NP),
  423. 1 (SIM(I,NP),I=1,IPTEM)
  424. IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
  425. END IF
  426. GOTO 140
  427. END IF
  428. C
  429. C Return the best calculated values of the variables.
  430. C
  431. IF (IPRINT .GE. 1) PRINT 590
  432. 590 FORMAT (/3X,'Normal return from subroutine COBYLA')
  433. IF (IFULL .EQ. 1) GOTO 620
  434. 600 DO 610 I=1,N
  435. 610 X(I)=SIM(I,NP)
  436. F=DATMAT(MP,NP)
  437. RESMAX=DATMAT(MPP,NP)
  438. 620 IF (IPRINT .GE. 1) THEN
  439. PRINT 70, NFVALS,F,RESMAX,(X(I),I=1,IPTEM)
  440. IF (IPTEM .LT. N) PRINT 80, (X(I),I=IPTEMP,N)
  441. END IF
  442. MAXFUN=NFVALS
  443. RETURN
  444. END