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