1;;; Compiled by f2cl version: 2;;; ("f2cl1.l,v 1.222 2010/10/08 03:05:30 rtoy Exp $" 3;;; "f2cl2.l,v 1.37 2008/02/22 22:19:33 rtoy Exp $" 4;;; "f2cl3.l,v 1.6 2008/02/22 22:19:33 rtoy Exp $" 5;;; "f2cl4.l,v 1.7 2008/02/22 22:19:34 rtoy Exp $" 6;;; "f2cl5.l,v 1.204 2010/02/23 05:21:30 rtoy Exp $" 7;;; "f2cl6.l,v 1.48 2008/08/24 00:56:27 rtoy Exp $" 8;;; "macros.l,v 1.114 2010/05/17 01:42:14 rtoy Exp $") 9 10;;; Using Lisp CMU Common Lisp CVS Head 2010-09-27 22:45:24 (20B Unicode) 11;;; 12;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t) 13;;; (:coerce-assigns :as-needed) (:array-type ':array) 14;;; (:array-slicing t) (:declare-common nil) 15;;; (:float-format double-float)) 16 17(in-package :cobyla) 18 19 20(defun cobylb 21 (n m mpp x rhobeg rhoend iprint maxfun con sim simi datmat a vsig veta 22 sigbar dx w iact ierr) 23 (declare (type (array f2cl-lib:integer4 (*)) iact) 24 (type double-float rhoend rhobeg) 25 (type (array double-float (*)) w dx sigbar veta vsig a datmat simi 26 sim con x) 27 (type (f2cl-lib:integer4) ierr maxfun iprint mpp m n)) 28 (f2cl-lib:with-multi-array-data 29 ((x double-float x-%data% x-%offset%) 30 (con double-float con-%data% con-%offset%) 31 (sim double-float sim-%data% sim-%offset%) 32 (simi double-float simi-%data% simi-%offset%) 33 (datmat double-float datmat-%data% datmat-%offset%) 34 (a double-float a-%data% a-%offset%) 35 (vsig double-float vsig-%data% vsig-%offset%) 36 (veta double-float veta-%data% veta-%offset%) 37 (sigbar double-float sigbar-%data% sigbar-%offset%) 38 (dx double-float dx-%data% dx-%offset%) 39 (w double-float w-%data% w-%offset%) 40 (iact f2cl-lib:integer4 iact-%data% iact-%offset%)) 41 (prog ((cmax 0.0) (cmin 0.0) (denom 0.0) (l 0) (edgmax 0.0) (ratio 0.0) 42 (trured 0.0) (vmnew 0.0) (vmold 0.0) (prerem 0.0) (phi 0.0) 43 (prerec 0.0) (barmu 0.0) (resnew 0.0) (ifull 0) (ivmd 0) (idxnew 0) 44 (isdirn 0) (ivmc 0) (izdota 0) (iz 0) (dxsign 0.0) (sum 0.0) 45 (cvmaxm 0.0) (cvmaxp 0.0) (weta 0.0) (wsig 0.0) (pareta 0.0) 46 (parsig 0.0) (iflag 0) (error$ 0.0) (tempa 0.0) (nbest 0) 47 (phimin 0.0) (k 0) (resmax 0.0) (f 0.0) (ibrnch 0) (jdrop 0) (j 0) 48 (i 0) (temp 0.0) (nfvals 0) (parmu 0.0) (rho 0.0) (delta 0.0) 49 (gamma 0.0) (beta 0.0) (alpha 0.0) (mp 0) (np 0) (iptemp 0) 50 (iptem 0)) 51 (declare (type (f2cl-lib:integer4) iptem iptemp np mp nfvals i j jdrop 52 ibrnch k nbest iflag iz izdota ivmc 53 isdirn idxnew ivmd ifull l) 54 (type double-float alpha beta gamma delta rho parmu temp f 55 resmax phimin tempa error$ parsig pareta wsig 56 weta cvmaxp cvmaxm sum dxsign resnew barmu 57 prerec phi prerem vmold vmnew trured ratio 58 edgmax denom cmin cmax)) 59 '"" 60 '" set the initial values of some parameters. the last column of sim" 61 '" the optimal vertex of the current simplex, and the preceding n col" 62 '" hold the displacements from the optimal vertex to the other vertic" 63 '" further, simi holds the inverse of the matrix that is contained in" 64 '" first n columns of sim." 65 '"" 66 (setf ierr 0) 67 (setf iptem (f2cl-lib:min0 n 5)) 68 (setf iptemp (f2cl-lib:int-add iptem 1)) 69 (setf np (f2cl-lib:int-add n 1)) 70 (setf mp (f2cl-lib:int-add m 1)) 71 (setf alpha 0.25) 72 (setf beta 2.1) 73 (setf gamma 0.5) 74 (setf delta 1.1) 75 (setf rho rhobeg) 76 (setf parmu 0.0) 77 (if (>= iprint 2) 78 (f2cl-lib:fformat t 79 ("~%" "~3@T" "The initial value of RHO is" 1 80 (("~13,6,2,1,'*,,'EE")) "~2@T" 81 "and PARMU is set to zero." "~%") 82 rho)) 83 (setf nfvals 0) 84 (setf temp (/ 1.0 rho)) 85 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 86 ((> i n) nil) 87 (tagbody 88 (setf (f2cl-lib:fref sim-%data% (i np) ((1 n) (1 *)) sim-%offset%) 89 (f2cl-lib:fref x-%data% (i) ((1 *)) x-%offset%)) 90 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 91 ((> j n) nil) 92 (tagbody 93 (setf (f2cl-lib:fref sim-%data% (i j) ((1 n) (1 *)) sim-%offset%) 94 0.0) 95 label20 96 (setf (f2cl-lib:fref simi-%data% 97 (i j) 98 ((1 n) (1 *)) 99 simi-%offset%) 100 0.0))) 101 (setf (f2cl-lib:fref sim-%data% (i i) ((1 n) (1 *)) sim-%offset%) 102 rho) 103 label30 104 (setf (f2cl-lib:fref simi-%data% (i i) ((1 n) (1 *)) simi-%offset%) 105 temp))) 106 (setf jdrop np) 107 (setf ibrnch 0) 108 '"" 109 '" make the next call of the user-supplied subroutine calcfc. these" 110 '" instructions are also used for calling calcfc during the iteration" 111 '" the algorithm." 112 '"" 113 label40 114 (cond 115 ((and (>= nfvals maxfun) (> nfvals 0)) 116 (if (>= iprint 1) 117 (f2cl-lib:fformat t 118 ("~%" "~3@T" 119 "Return from subroutine COBYLA because the " 120 "MAXFUN limit has been reached." "~%"))) 121 (setf ierr 1) 122 (go label600))) 123 (setf nfvals (f2cl-lib:int-add nfvals 1)) 124 (multiple-value-bind (var-0 var-1 var-2 var-3 var-4) 125 (calcfc n m x f con) 126 (declare (ignore var-0 var-1 var-2 var-4)) 127 (setf f var-3)) 128 (setf resmax 0.0) 129 (cond 130 ((> m 0) 131 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 132 ((> k m) nil) 133 (tagbody 134 label60 135 (setf resmax 136 (f2cl-lib:dmax1 resmax 137 (- 138 (f2cl-lib:fref con-%data% 139 (k) 140 ((1 *)) 141 con-%offset%)))))))) 142 (cond 143 ((or (= nfvals (f2cl-lib:int-add iprint (f2cl-lib:int-sub 1))) 144 (= iprint 3)) 145 (f2cl-lib:fformat t 146 ("~%" "~3@T" "NFVALS =" 1 (("~5D")) "~3@T" "F =" 1 147 (("~13,6,2,1,'*,,'EE")) "~4@T" "MAXCV =" 1 148 (("~13,6,2,1,'*,,'EE")) "~%" "~3@T" "X =" 1 149 (("~13,6,2,1,'*,,'EE")) 4 (("~15,6,2,1,'*,,'EE")) 150 "~%") 151 nfvals 152 f 153 resmax 154 (do ((i 1 (f2cl-lib:int-add i 1)) 155 (%ret nil)) 156 ((> i iptem) (nreverse %ret)) 157 (declare (type f2cl-lib:integer4 i)) 158 (push 159 (f2cl-lib:fref x-%data% (i) ((1 *)) x-%offset%) 160 %ret))) 161 (if (< iptem n) 162 (f2cl-lib:fformat t 163 (1 (("~19,6,2,1,'*,,'EE")) 4 164 (("~15,6,2,1,'*,,'EE")) "~%") 165 (do ((i iptemp (f2cl-lib:int-add i 1)) 166 (%ret nil)) 167 ((> i n) (nreverse %ret)) 168 (declare (type f2cl-lib:integer4 i)) 169 (push 170 (f2cl-lib:fref x-%data% 171 (i) 172 ((1 *)) 173 x-%offset%) 174 %ret)))))) 175 (setf (f2cl-lib:fref con-%data% (mp) ((1 *)) con-%offset%) f) 176 (setf (f2cl-lib:fref con-%data% (mpp) ((1 *)) con-%offset%) resmax) 177 (if (= ibrnch 1) (go label440)) 178 '"" 179 '" set the recently calculated function values in a column of datmat." 180 '" array has a column for each vertex of the current simplex, the ent" 181 '" each column being the values of the constraint functions (if any)" 182 '" followed by the objective function and the greatest constraint vio" 183 '" at the vertex." 184 '"" 185 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 186 ((> k mpp) nil) 187 (tagbody 188 label90 189 (setf (f2cl-lib:fref datmat-%data% 190 (k jdrop) 191 ((1 mpp) (1 *)) 192 datmat-%offset%) 193 (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%)))) 194 (if (> nfvals np) (go label130)) 195 '"" 196 '" exchange the new vertex of the initial simplex with the optimal ve" 197 '" necessary. then, if the initial simplex is not complete, pick its" 198 '" vertex and calculate the function values there." 199 '"" 200 (cond 201 ((<= jdrop n) 202 (cond 203 ((<= (f2cl-lib:fref datmat (mp np) ((1 mpp) (1 *))) f) 204 (setf (f2cl-lib:fref x-%data% (jdrop) ((1 *)) x-%offset%) 205 (f2cl-lib:fref sim-%data% 206 (jdrop np) 207 ((1 n) (1 *)) 208 sim-%offset%))) 209 (t 210 (setf (f2cl-lib:fref sim-%data% 211 (jdrop np) 212 ((1 n) (1 *)) 213 sim-%offset%) 214 (f2cl-lib:fref x-%data% (jdrop) ((1 *)) x-%offset%)) 215 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 216 ((> k mpp) nil) 217 (tagbody 218 (setf (f2cl-lib:fref datmat-%data% 219 (k jdrop) 220 ((1 mpp) (1 *)) 221 datmat-%offset%) 222 (f2cl-lib:fref datmat-%data% 223 (k np) 224 ((1 mpp) (1 *)) 225 datmat-%offset%)) 226 label100 227 (setf (f2cl-lib:fref datmat-%data% 228 (k np) 229 ((1 mpp) (1 *)) 230 datmat-%offset%) 231 (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%)))) 232 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 233 ((> k jdrop) nil) 234 (tagbody 235 (setf (f2cl-lib:fref sim-%data% 236 (jdrop k) 237 ((1 n) (1 *)) 238 sim-%offset%) 239 (- rho)) 240 (setf temp (coerce 0.0f0 'double-float)) 241 (f2cl-lib:fdo (i k (f2cl-lib:int-add i 1)) 242 ((> i jdrop) nil) 243 (tagbody 244 label110 245 (setf temp 246 (- temp 247 (f2cl-lib:fref simi-%data% 248 (i k) 249 ((1 n) (1 *)) 250 simi-%offset%))))) 251 label120 252 (setf (f2cl-lib:fref simi-%data% 253 (jdrop k) 254 ((1 n) (1 *)) 255 simi-%offset%) 256 temp))))))) 257 (cond 258 ((<= nfvals n) 259 (setf jdrop nfvals) 260 (setf (f2cl-lib:fref x-%data% (jdrop) ((1 *)) x-%offset%) 261 (+ (f2cl-lib:fref x-%data% (jdrop) ((1 *)) x-%offset%) rho)) 262 (go label40))) 263 label130 264 (setf ibrnch 1) 265 '"" 266 '" identify the optimal vertex of the current simplex." 267 '"" 268 label140 269 (setf phimin 270 (+ 271 (f2cl-lib:fref datmat-%data% 272 (mp np) 273 ((1 mpp) (1 *)) 274 datmat-%offset%) 275 (* parmu 276 (f2cl-lib:fref datmat-%data% 277 (mpp np) 278 ((1 mpp) (1 *)) 279 datmat-%offset%)))) 280 (setf nbest np) 281 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 282 ((> j n) nil) 283 (tagbody 284 (setf temp 285 (+ 286 (f2cl-lib:fref datmat-%data% 287 (mp j) 288 ((1 mpp) (1 *)) 289 datmat-%offset%) 290 (* parmu 291 (f2cl-lib:fref datmat-%data% 292 (mpp j) 293 ((1 mpp) (1 *)) 294 datmat-%offset%)))) 295 (cond 296 ((< temp phimin) 297 (setf nbest j) 298 (setf phimin temp)) 299 ((and (= temp phimin) (= parmu 0.0)) 300 (if 301 (< 302 (f2cl-lib:fref datmat-%data% 303 (mpp j) 304 ((1 mpp) (1 *)) 305 datmat-%offset%) 306 (f2cl-lib:fref datmat-%data% 307 (mpp nbest) 308 ((1 mpp) (1 *)) 309 datmat-%offset%)) 310 (setf nbest j)))) 311 label150)) 312 '"" 313 '" switch the best vertex into pole position if it is not there alrea" 314 '" and also update sim, simi and datmat." 315 '"" 316 (cond 317 ((<= nbest n) 318 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 319 ((> i mpp) nil) 320 (tagbody 321 (setf temp 322 (f2cl-lib:fref datmat-%data% 323 (i np) 324 ((1 mpp) (1 *)) 325 datmat-%offset%)) 326 (setf (f2cl-lib:fref datmat-%data% 327 (i np) 328 ((1 mpp) (1 *)) 329 datmat-%offset%) 330 (f2cl-lib:fref datmat-%data% 331 (i nbest) 332 ((1 mpp) (1 *)) 333 datmat-%offset%)) 334 label160 335 (setf (f2cl-lib:fref datmat-%data% 336 (i nbest) 337 ((1 mpp) (1 *)) 338 datmat-%offset%) 339 temp))) 340 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 341 ((> i n) nil) 342 (tagbody 343 (setf temp 344 (f2cl-lib:fref sim-%data% 345 (i nbest) 346 ((1 n) (1 *)) 347 sim-%offset%)) 348 (setf (f2cl-lib:fref sim-%data% 349 (i nbest) 350 ((1 n) (1 *)) 351 sim-%offset%) 352 0.0) 353 (setf (f2cl-lib:fref sim-%data% (i np) ((1 n) (1 *)) sim-%offset%) 354 (+ 355 (f2cl-lib:fref sim-%data% 356 (i np) 357 ((1 n) (1 *)) 358 sim-%offset%) 359 temp)) 360 (setf tempa 0.0) 361 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 362 ((> k n) nil) 363 (tagbody 364 (setf (f2cl-lib:fref sim-%data% 365 (i k) 366 ((1 n) (1 *)) 367 sim-%offset%) 368 (- 369 (f2cl-lib:fref sim-%data% 370 (i k) 371 ((1 n) (1 *)) 372 sim-%offset%) 373 temp)) 374 label170 375 (setf tempa 376 (- tempa 377 (f2cl-lib:fref simi-%data% 378 (k i) 379 ((1 n) (1 *)) 380 simi-%offset%))))) 381 label180 382 (setf (f2cl-lib:fref simi-%data% 383 (nbest i) 384 ((1 n) (1 *)) 385 simi-%offset%) 386 tempa))))) 387 '"" 388 '" make an error return if sigi is a poor approximation to the invers" 389 '" the leading n by n submatrix of sig." 390 '"" 391 (setf error$ 0.0) 392 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 393 ((> i n) nil) 394 (tagbody 395 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 396 ((> j n) nil) 397 (tagbody 398 (setf temp 0.0) 399 (if (= i j) (setf temp (- temp 1.0))) 400 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 401 ((> k n) nil) 402 (tagbody 403 label190 404 (setf temp 405 (+ temp 406 (* 407 (f2cl-lib:fref simi-%data% 408 (i k) 409 ((1 n) (1 *)) 410 simi-%offset%) 411 (f2cl-lib:fref sim-%data% 412 (k j) 413 ((1 n) (1 *)) 414 sim-%offset%)))))) 415 (setf error$ (f2cl-lib:dmax1 error$ (abs temp))))))) 416 label200 417 (cond 418 ((> error$ 0.1) 419 (if (>= iprint 1) 420 (f2cl-lib:fformat t 421 ("~%" "~3@T" 422 "Return from subroutine COBYLA because " 423 "rounding errors are becoming damaging." "~%"))) 424 (setf ierr 2) 425 (go label600))) 426 '"" 427 '" calculate the coefficients of the linear approximations to the obj" 428 '" and constraint functions, placing minus the objective function gra" 429 '" after the constraint gradients in the array a. the vector w is use" 430 '" working space." 431 '"" 432 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 433 ((> k mp) nil) 434 (tagbody 435 (setf (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%) 436 (- 437 (f2cl-lib:fref datmat-%data% 438 (k np) 439 ((1 mpp) (1 *)) 440 datmat-%offset%))) 441 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 442 ((> j n) nil) 443 (tagbody 444 label220 445 (setf (f2cl-lib:fref w-%data% (j) ((1 *)) w-%offset%) 446 (+ 447 (f2cl-lib:fref datmat-%data% 448 (k j) 449 ((1 mpp) (1 *)) 450 datmat-%offset%) 451 (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%))))) 452 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 453 ((> i n) nil) 454 (tagbody 455 (setf temp 0.0) 456 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 457 ((> j n) nil) 458 (tagbody 459 label230 460 (setf temp 461 (+ temp 462 (* (f2cl-lib:fref w-%data% (j) ((1 *)) w-%offset%) 463 (f2cl-lib:fref simi-%data% 464 (j i) 465 ((1 n) (1 *)) 466 simi-%offset%)))))) 467 (if (= k mp) (setf temp (- temp))) 468 (setf (f2cl-lib:fref a-%data% (i k) ((1 n) (1 *)) a-%offset%) 469 temp))))) 470 label240 471 '"" 472 '" calculate the values of sigma and eta, and set iflag=0 if the curr" 473 '" simplex is not acceptable." 474 '"" 475 (setf iflag 1) 476 (setf parsig (* alpha rho)) 477 (setf pareta (* beta rho)) 478 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 479 ((> j n) nil) 480 (tagbody 481 (setf wsig 0.0) 482 (setf weta 0.0) 483 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 484 ((> i n) nil) 485 (tagbody 486 (setf wsig 487 (+ wsig 488 (expt 489 (f2cl-lib:fref simi-%data% 490 (j i) 491 ((1 n) (1 *)) 492 simi-%offset%) 493 2))) 494 label250 495 (setf weta 496 (+ weta 497 (expt 498 (f2cl-lib:fref sim-%data% 499 (i j) 500 ((1 n) (1 *)) 501 sim-%offset%) 502 2))))) 503 (setf (f2cl-lib:fref vsig-%data% (j) ((1 *)) vsig-%offset%) 504 (/ 1.0 (f2cl-lib:fsqrt wsig))) 505 (setf (f2cl-lib:fref veta-%data% (j) ((1 *)) veta-%offset%) 506 (f2cl-lib:fsqrt weta)) 507 (if 508 (or (< (f2cl-lib:fref vsig-%data% (j) ((1 *)) vsig-%offset%) parsig) 509 (> (f2cl-lib:fref veta-%data% (j) ((1 *)) veta-%offset%) 510 pareta)) 511 (setf iflag 0)) 512 label260)) 513 '"" 514 '" if a new vertex is needed to improve acceptability, then decide wh" 515 '" vertex to drop from the simplex." 516 '"" 517 (if (or (= ibrnch 1) (= iflag 1)) (go label370)) 518 (setf jdrop 0) 519 (setf temp pareta) 520 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 521 ((> j n) nil) 522 (tagbody 523 (cond 524 ((> (f2cl-lib:fref veta (j) ((1 *))) temp) 525 (setf jdrop j) 526 (setf temp 527 (f2cl-lib:fref veta-%data% (j) ((1 *)) veta-%offset%)))) 528 label270)) 529 (cond 530 ((= jdrop 0) 531 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 532 ((> j n) nil) 533 (tagbody 534 (cond 535 ((< (f2cl-lib:fref vsig (j) ((1 *))) temp) 536 (setf jdrop j) 537 (setf temp 538 (f2cl-lib:fref vsig-%data% 539 (j) 540 ((1 *)) 541 vsig-%offset%)))) 542 label280)))) 543 '"" 544 '" calculate the step to the new vertex and its sign." 545 '"" 546 (setf temp 547 (* gamma 548 rho 549 (f2cl-lib:fref vsig-%data% (jdrop) ((1 *)) vsig-%offset%))) 550 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 551 ((> i n) nil) 552 (tagbody 553 label290 554 (setf (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%) 555 (* temp 556 (f2cl-lib:fref simi-%data% 557 (jdrop i) 558 ((1 n) (1 *)) 559 simi-%offset%))))) 560 (setf cvmaxp 0.0) 561 (setf cvmaxm 0.0) 562 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 563 ((> k mp) nil) 564 (tagbody 565 (setf sum 0.0) 566 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 567 ((> i n) nil) 568 (tagbody 569 label300 570 (setf sum 571 (+ sum 572 (* 573 (f2cl-lib:fref a-%data% 574 (i k) 575 ((1 n) (1 *)) 576 a-%offset%) 577 (f2cl-lib:fref dx-%data% 578 (i) 579 ((1 *)) 580 dx-%offset%)))))) 581 (cond 582 ((< k mp) 583 (setf temp 584 (f2cl-lib:fref datmat-%data% 585 (k np) 586 ((1 mpp) (1 *)) 587 datmat-%offset%)) 588 (setf cvmaxp (f2cl-lib:dmax1 cvmaxp (- (- temp) sum))) 589 (setf cvmaxm (f2cl-lib:dmax1 cvmaxm (- sum temp))))) 590 label310)) 591 (setf dxsign 1.0) 592 (if (> (* parmu (- cvmaxp cvmaxm)) (+ sum sum)) (setf dxsign -1.0)) 593 '"" 594 '" update the elements of sim and simi, and set the next x." 595 '"" 596 (setf temp 0.0) 597 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 598 ((> i n) nil) 599 (tagbody 600 (setf (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%) 601 (* dxsign (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%))) 602 (setf (f2cl-lib:fref sim-%data% (i jdrop) ((1 n) (1 *)) sim-%offset%) 603 (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)) 604 label320 605 (setf temp 606 (+ temp 607 (* 608 (f2cl-lib:fref simi-%data% 609 (jdrop i) 610 ((1 n) (1 *)) 611 simi-%offset%) 612 (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)))))) 613 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 614 ((> i n) nil) 615 (tagbody 616 label330 617 (setf (f2cl-lib:fref simi-%data% 618 (jdrop i) 619 ((1 n) (1 *)) 620 simi-%offset%) 621 (/ 622 (f2cl-lib:fref simi-%data% 623 (jdrop i) 624 ((1 n) (1 *)) 625 simi-%offset%) 626 temp)))) 627 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 628 ((> j n) nil) 629 (tagbody 630 (cond 631 ((/= j jdrop) 632 (setf temp 0.0) 633 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 634 ((> i n) nil) 635 (tagbody 636 label340 637 (setf temp 638 (+ temp 639 (* 640 (f2cl-lib:fref simi-%data% 641 (j i) 642 ((1 n) (1 *)) 643 simi-%offset%) 644 (f2cl-lib:fref dx-%data% 645 (i) 646 ((1 *)) 647 dx-%offset%)))))) 648 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 649 ((> i n) nil) 650 (tagbody 651 label350 652 (setf (f2cl-lib:fref simi-%data% 653 (j i) 654 ((1 n) (1 *)) 655 simi-%offset%) 656 (- 657 (f2cl-lib:fref simi-%data% 658 (j i) 659 ((1 n) (1 *)) 660 simi-%offset%) 661 (* temp 662 (f2cl-lib:fref simi-%data% 663 (jdrop i) 664 ((1 n) (1 *)) 665 simi-%offset%)))))))) 666 label360 667 (setf (f2cl-lib:fref x-%data% (j) ((1 *)) x-%offset%) 668 (+ 669 (f2cl-lib:fref sim-%data% (j np) ((1 n) (1 *)) sim-%offset%) 670 (f2cl-lib:fref dx-%data% (j) ((1 *)) dx-%offset%))))) 671 (go label40) 672 '"" 673 '" calculate dx=x(*)-x(0). branch if the length of dx is less than 0." 674 '"" 675 label370 676 (setf iz 1) 677 (setf izdota (f2cl-lib:int-add iz (f2cl-lib:int-mul n n))) 678 (setf ivmc (f2cl-lib:int-add izdota n)) 679 (setf isdirn (f2cl-lib:int-add ivmc mp)) 680 (setf idxnew (f2cl-lib:int-add isdirn n)) 681 (setf ivmd (f2cl-lib:int-add idxnew n)) 682 (multiple-value-bind 683 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10 684 var-11 var-12 var-13) 685 (trstlp n m a con rho dx ifull iact 686 (f2cl-lib:array-slice w double-float (iz) ((1 *))) 687 (f2cl-lib:array-slice w double-float (izdota) ((1 *))) 688 (f2cl-lib:array-slice w double-float (ivmc) ((1 *))) 689 (f2cl-lib:array-slice w double-float (isdirn) ((1 *))) 690 (f2cl-lib:array-slice w double-float (idxnew) ((1 *))) 691 (f2cl-lib:array-slice w double-float (ivmd) ((1 *)))) 692 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-7 var-8 var-9 693 var-10 var-11 var-12 var-13)) 694 (setf ifull var-6)) 695 (cond 696 ((= ifull 0) 697 (setf temp 0.0) 698 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 699 ((> i n) nil) 700 (tagbody 701 label380 702 (setf temp 703 (+ temp 704 (expt (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%) 705 2))))) 706 (cond 707 ((< temp (* 0.25 rho rho)) 708 (setf ibrnch 1) 709 (go label550))))) 710 '"" 711 '" predict the change to f and the new maximum constraint violation i" 712 '" variables are altered from x(0) to x(0)+dx." 713 '"" 714 (setf resnew 0.0) 715 (setf (f2cl-lib:fref con-%data% (mp) ((1 *)) con-%offset%) 0.0) 716 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 717 ((> k mp) nil) 718 (tagbody 719 (setf sum (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%)) 720 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 721 ((> i n) nil) 722 (tagbody 723 label390 724 (setf sum 725 (- sum 726 (* 727 (f2cl-lib:fref a-%data% 728 (i k) 729 ((1 n) (1 *)) 730 a-%offset%) 731 (f2cl-lib:fref dx-%data% 732 (i) 733 ((1 *)) 734 dx-%offset%)))))) 735 (if (< k mp) (setf resnew (f2cl-lib:dmax1 resnew sum))) 736 label400)) 737 '"" 738 '" increase parmu if necessary and branch back if this change alters" 739 '" optimal vertex. otherwise prerem and prerec will be set to the pre" 740 '" reductions in the merit function and the maximum constraint violat" 741 '" respectively." 742 '"" 743 (setf barmu 0.0) 744 (setf prerec 745 (- 746 (f2cl-lib:fref datmat-%data% 747 (mpp np) 748 ((1 mpp) (1 *)) 749 datmat-%offset%) 750 resnew)) 751 (if (> prerec 0.0) (setf barmu (/ sum prerec))) 752 (cond 753 ((< parmu (* 1.5 barmu)) 754 (setf parmu (* 2.0 barmu)) 755 (if (>= iprint 2) 756 (f2cl-lib:fformat t 757 ("~%" "~3@T" "Increase in PARMU to" 1 758 (("~13,6,2,1,'*,,'EE")) "~%") 759 parmu)) 760 (setf phi 761 (+ 762 (f2cl-lib:fref datmat-%data% 763 (mp np) 764 ((1 mpp) (1 *)) 765 datmat-%offset%) 766 (* parmu 767 (f2cl-lib:fref datmat-%data% 768 (mpp np) 769 ((1 mpp) (1 *)) 770 datmat-%offset%)))) 771 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 772 ((> j n) nil) 773 (tagbody 774 (setf temp 775 (+ 776 (f2cl-lib:fref datmat-%data% 777 (mp j) 778 ((1 mpp) (1 *)) 779 datmat-%offset%) 780 (* parmu 781 (f2cl-lib:fref datmat-%data% 782 (mpp j) 783 ((1 mpp) (1 *)) 784 datmat-%offset%)))) 785 (if (< temp phi) (go label140)) 786 (cond 787 ((and (= temp phi) (= parmu 0.0f0)) 788 (if 789 (< 790 (f2cl-lib:fref datmat-%data% 791 (mpp j) 792 ((1 mpp) (1 *)) 793 datmat-%offset%) 794 (f2cl-lib:fref datmat-%data% 795 (mpp np) 796 ((1 mpp) (1 *)) 797 datmat-%offset%)) 798 (go label140)))) 799 label420)))) 800 (setf prerem (- (* parmu prerec) sum)) 801 '"" 802 '" calculate the constraint and objective functions at x(*). then fin" 803 '" actual reduction in the merit function." 804 '"" 805 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 806 ((> i n) nil) 807 (tagbody 808 label430 809 (setf (f2cl-lib:fref x-%data% (i) ((1 *)) x-%offset%) 810 (+ 811 (f2cl-lib:fref sim-%data% (i np) ((1 n) (1 *)) sim-%offset%) 812 (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%))))) 813 (setf ibrnch 1) 814 (go label40) 815 label440 816 (setf vmold 817 (+ 818 (f2cl-lib:fref datmat-%data% 819 (mp np) 820 ((1 mpp) (1 *)) 821 datmat-%offset%) 822 (* parmu 823 (f2cl-lib:fref datmat-%data% 824 (mpp np) 825 ((1 mpp) (1 *)) 826 datmat-%offset%)))) 827 (setf vmnew (+ f (* parmu resmax))) 828 (setf trured (- vmold vmnew)) 829 (cond 830 ((and (= parmu 0.0) 831 (= f (f2cl-lib:fref datmat (mp np) ((1 mpp) (1 *))))) 832 (setf prerem prerec) 833 (setf trured 834 (- 835 (f2cl-lib:fref datmat-%data% 836 (mpp np) 837 ((1 mpp) (1 *)) 838 datmat-%offset%) 839 resmax)))) 840 '"" 841 '" begin the operations that decide whether x(*) should replace one o" 842 '" vertices of the current simplex, the change being mandatory if tru" 843 '" positive. firstly, jdrop is set to the index of the vertex that is" 844 '" replaced." 845 '"" 846 (setf ratio 0.0) 847 (if (<= trured 0.0f0) (setf ratio (coerce 1.0f0 'double-float))) 848 (setf jdrop 0) 849 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 850 ((> j n) nil) 851 (tagbody 852 (setf temp 0.0) 853 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 854 ((> i n) nil) 855 (tagbody 856 label450 857 (setf temp 858 (+ temp 859 (* 860 (f2cl-lib:fref simi-%data% 861 (j i) 862 ((1 n) (1 *)) 863 simi-%offset%) 864 (f2cl-lib:fref dx-%data% 865 (i) 866 ((1 *)) 867 dx-%offset%)))))) 868 (setf temp (abs temp)) 869 (cond 870 ((> temp ratio) 871 (setf jdrop j) 872 (setf ratio temp))) 873 label460 874 (setf (f2cl-lib:fref sigbar-%data% (j) ((1 *)) sigbar-%offset%) 875 (* temp 876 (f2cl-lib:fref vsig-%data% (j) ((1 *)) vsig-%offset%))))) 877 '"" 878 '" calculate the value of ell." 879 '"" 880 (setf edgmax (* delta rho)) 881 (setf l 0) 882 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 883 ((> j n) nil) 884 (tagbody 885 (cond 886 ((or (>= (f2cl-lib:fref sigbar (j) ((1 *))) parsig) 887 (>= (f2cl-lib:fref sigbar (j) ((1 *))) 888 (f2cl-lib:fref vsig (j) ((1 *))))) 889 (setf temp (f2cl-lib:fref veta-%data% (j) ((1 *)) veta-%offset%)) 890 (cond 891 ((> trured 0.0) 892 (setf temp 0.0) 893 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 894 ((> i n) nil) 895 (tagbody 896 label470 897 (setf temp 898 (+ temp 899 (expt 900 (- 901 (f2cl-lib:fref dx-%data% 902 (i) 903 ((1 *)) 904 dx-%offset%) 905 (f2cl-lib:fref sim-%data% 906 (i j) 907 ((1 n) (1 *)) 908 sim-%offset%)) 909 2))))) 910 (setf temp (f2cl-lib:fsqrt temp)))) 911 (cond 912 ((> temp edgmax) 913 (setf l j) 914 (setf edgmax temp))))) 915 label480)) 916 (if (> l 0) (setf jdrop l)) 917 (if (= jdrop 0) (go label550)) 918 '"" 919 '" revise the simplex by updating the elements of sim, simi and datma" 920 '"" 921 (setf temp 0.0) 922 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 923 ((> i n) nil) 924 (tagbody 925 (setf (f2cl-lib:fref sim-%data% (i jdrop) ((1 n) (1 *)) sim-%offset%) 926 (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)) 927 label490 928 (setf temp 929 (+ temp 930 (* 931 (f2cl-lib:fref simi-%data% 932 (jdrop i) 933 ((1 n) (1 *)) 934 simi-%offset%) 935 (f2cl-lib:fref dx-%data% (i) ((1 *)) dx-%offset%)))))) 936 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 937 ((> i n) nil) 938 (tagbody 939 label500 940 (setf (f2cl-lib:fref simi-%data% 941 (jdrop i) 942 ((1 n) (1 *)) 943 simi-%offset%) 944 (/ 945 (f2cl-lib:fref simi-%data% 946 (jdrop i) 947 ((1 n) (1 *)) 948 simi-%offset%) 949 temp)))) 950 (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1)) 951 ((> j n) nil) 952 (tagbody 953 (cond 954 ((/= j jdrop) 955 (setf temp 0.0) 956 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 957 ((> i n) nil) 958 (tagbody 959 label510 960 (setf temp 961 (+ temp 962 (* 963 (f2cl-lib:fref simi-%data% 964 (j i) 965 ((1 n) (1 *)) 966 simi-%offset%) 967 (f2cl-lib:fref dx-%data% 968 (i) 969 ((1 *)) 970 dx-%offset%)))))) 971 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 972 ((> i n) nil) 973 (tagbody 974 label520 975 (setf (f2cl-lib:fref simi-%data% 976 (j i) 977 ((1 n) (1 *)) 978 simi-%offset%) 979 (- 980 (f2cl-lib:fref simi-%data% 981 (j i) 982 ((1 n) (1 *)) 983 simi-%offset%) 984 (* temp 985 (f2cl-lib:fref simi-%data% 986 (jdrop i) 987 ((1 n) (1 *)) 988 simi-%offset%)))))))) 989 label530)) 990 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 991 ((> k mpp) nil) 992 (tagbody 993 label540 994 (setf (f2cl-lib:fref datmat-%data% 995 (k jdrop) 996 ((1 mpp) (1 *)) 997 datmat-%offset%) 998 (f2cl-lib:fref con-%data% (k) ((1 *)) con-%offset%)))) 999 '"" 1000 '" branch back for further iterations with the current rho." 1001 '"" 1002 (if (and (> trured 0.0) (>= trured (* 0.1 prerem))) (go label140)) 1003 label550 1004 (cond 1005 ((= iflag 0) 1006 (setf ibrnch 0) 1007 (go label140))) 1008 '"" 1009 '" otherwise reduce rho if it is not at its least value and reset par" 1010 '"" 1011 (cond 1012 ((> rho rhoend) 1013 (setf rho (* 0.5 rho)) 1014 (if (<= rho (* 1.5 rhoend)) (setf rho rhoend)) 1015 (cond 1016 ((> parmu 0.0) 1017 (setf denom 0.0) 1018 (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1)) 1019 ((> k mp) nil) 1020 (tagbody 1021 (setf cmin 1022 (f2cl-lib:fref datmat-%data% 1023 (k np) 1024 ((1 mpp) (1 *)) 1025 datmat-%offset%)) 1026 (setf cmax cmin) 1027 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 1028 ((> i n) nil) 1029 (tagbody 1030 (setf cmin 1031 (f2cl-lib:dmin1 cmin 1032 (f2cl-lib:fref datmat-%data% 1033 (k i) 1034 ((1 mpp) (1 *)) 1035 datmat-%offset%))) 1036 label560 1037 (setf cmax 1038 (f2cl-lib:dmax1 cmax 1039 (f2cl-lib:fref datmat-%data% 1040 (k i) 1041 ((1 mpp) (1 *)) 1042 datmat-%offset%))))) 1043 (cond 1044 ((and (<= k m) (< cmin (* 0.5 cmax))) 1045 (setf temp (- (f2cl-lib:dmax1 cmax 0.0) cmin)) 1046 (cond 1047 ((<= denom 0.0) 1048 (setf denom temp)) 1049 (t 1050 (setf denom (f2cl-lib:dmin1 denom temp)))))) 1051 label570)) 1052 (cond 1053 ((= denom 0.0) 1054 (setf parmu 0.0)) 1055 ((< (+ cmax (- cmin)) (* parmu denom)) 1056 (setf parmu (/ (- cmax cmin) denom)))))) 1057 (if (>= iprint 2) 1058 (f2cl-lib:fformat t 1059 ("~%" "~3@T" "Reduction in RHO to" 1 1060 (("~13,6,2,1,'*,,'EE")) " and PARMU =" 1 1061 (("~13,6,2,1,'*,,'EE")) "~%") 1062 rho 1063 parmu)) 1064 (cond 1065 ((= iprint 2) 1066 (f2cl-lib:fformat t 1067 ("~%" "~3@T" "NFVALS =" 1 (("~5D")) "~3@T" "F =" 1068 1 (("~13,6,2,1,'*,,'EE")) "~4@T" "MAXCV =" 1 1069 (("~13,6,2,1,'*,,'EE")) "~%" "~3@T" "X =" 1 1070 (("~13,6,2,1,'*,,'EE")) 4 1071 (("~15,6,2,1,'*,,'EE")) "~%") 1072 nfvals 1073 (f2cl-lib:fref datmat-%data% 1074 (mp np) 1075 ((1 mpp) (1 *)) 1076 datmat-%offset%) 1077 (f2cl-lib:fref datmat-%data% 1078 (mpp np) 1079 ((1 mpp) (1 *)) 1080 datmat-%offset%) 1081 (do ((i 1 (f2cl-lib:int-add i 1)) 1082 (%ret nil)) 1083 ((> i iptem) (nreverse %ret)) 1084 (declare (type f2cl-lib:integer4 i)) 1085 (push 1086 (f2cl-lib:fref sim-%data% 1087 (i np) 1088 ((1 n) (1 *)) 1089 sim-%offset%) 1090 %ret))) 1091 (if (< iptem n) 1092 (f2cl-lib:fformat t 1093 (1 (("~19,6,2,1,'*,,'EE")) 4 1094 (("~15,6,2,1,'*,,'EE")) "~%") 1095 (do ((i iptemp (f2cl-lib:int-add i 1)) 1096 (%ret nil)) 1097 ((> i n) (nreverse %ret)) 1098 (declare (type f2cl-lib:integer4 i)) 1099 (push 1100 (f2cl-lib:fref x-%data% 1101 (i) 1102 ((1 *)) 1103 x-%offset%) 1104 %ret)))))) 1105 (go label140))) 1106 '"" 1107 '" return the best calculated values of the variables." 1108 '"" 1109 (if (>= iprint 1) 1110 (f2cl-lib:fformat t 1111 ("~%" "~3@T" "Normal return from subroutine COBYLA" 1112 "~%"))) 1113 (if (= ifull 1) (go label620)) 1114 label600 1115 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1)) 1116 ((> i n) nil) 1117 (tagbody 1118 label610 1119 (setf (f2cl-lib:fref x-%data% (i) ((1 *)) x-%offset%) 1120 (f2cl-lib:fref sim-%data% 1121 (i np) 1122 ((1 n) (1 *)) 1123 sim-%offset%)))) 1124 (setf f 1125 (f2cl-lib:fref datmat-%data% 1126 (mp np) 1127 ((1 mpp) (1 *)) 1128 datmat-%offset%)) 1129 (setf resmax 1130 (f2cl-lib:fref datmat-%data% 1131 (mpp np) 1132 ((1 mpp) (1 *)) 1133 datmat-%offset%)) 1134 label620 1135 (cond 1136 ((>= iprint 1) 1137 (f2cl-lib:fformat t 1138 ("~%" "~3@T" "NFVALS =" 1 (("~5D")) "~3@T" "F =" 1 1139 (("~13,6,2,1,'*,,'EE")) "~4@T" "MAXCV =" 1 1140 (("~13,6,2,1,'*,,'EE")) "~%" "~3@T" "X =" 1 1141 (("~13,6,2,1,'*,,'EE")) 4 (("~15,6,2,1,'*,,'EE")) 1142 "~%") 1143 nfvals 1144 f 1145 resmax 1146 (do ((i 1 (f2cl-lib:int-add i 1)) 1147 (%ret nil)) 1148 ((> i iptem) (nreverse %ret)) 1149 (declare (type f2cl-lib:integer4 i)) 1150 (push 1151 (f2cl-lib:fref x-%data% (i) ((1 *)) x-%offset%) 1152 %ret))) 1153 (if (< iptem n) 1154 (f2cl-lib:fformat t 1155 (1 (("~19,6,2,1,'*,,'EE")) 4 1156 (("~15,6,2,1,'*,,'EE")) "~%") 1157 (do ((i iptemp (f2cl-lib:int-add i 1)) 1158 (%ret nil)) 1159 ((> i n) (nreverse %ret)) 1160 (declare (type f2cl-lib:integer4 i)) 1161 (push 1162 (f2cl-lib:fref x-%data% 1163 (i) 1164 ((1 *)) 1165 x-%offset%) 1166 %ret)))))) 1167 (setf maxfun nfvals) 1168 (go end_label) 1169 end_label 1170 (return 1171 (values nil 1172 nil 1173 nil 1174 nil 1175 nil 1176 nil 1177 nil 1178 maxfun 1179 nil 1180 nil 1181 nil 1182 nil 1183 nil 1184 nil 1185 nil 1186 nil 1187 nil 1188 nil 1189 nil 1190 ierr))))) 1191 1192(in-package #-gcl #:cl-user #+gcl "CL-USER") 1193#+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or)) 1194(eval-when (:load-toplevel :compile-toplevel :execute) 1195 (setf (gethash 'fortran-to-lisp::cobylb 1196 fortran-to-lisp::*f2cl-function-info*) 1197 (fortran-to-lisp::make-f2cl-finfo 1198 :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4) 1199 (fortran-to-lisp::integer4) (array double-float (*)) 1200 double-float double-float (fortran-to-lisp::integer4) 1201 (fortran-to-lisp::integer4) (array double-float (*)) 1202 (array double-float (*)) (array double-float (*)) 1203 (array double-float (*)) (array double-float (*)) 1204 (array double-float (*)) (array double-float (*)) 1205 (array double-float (*)) (array double-float (*)) 1206 (array double-float (*)) 1207 (array fortran-to-lisp::integer4 (*)) 1208 (fortran-to-lisp::integer4)) 1209 :return-values '(nil nil nil nil nil nil nil fortran-to-lisp::maxfun 1210 nil nil nil nil nil nil nil nil nil nil nil 1211 fortran-to-lisp::ierr) 1212 :calls '(fortran-to-lisp::trstlp fortran-to-lisp::calcfc)))) 1213 1214