Description of maxima_unicode_display in reference manual,
[maxima.git] / share / lbfgs / lbfgs.lisp
blob9044aaea962d44a5d1809116f26682b825c12426
1 ;;; Compiled by f2cl version:
2 ;;; ("f2cl1.l,v 1.215 2009/04/07 22:05:21 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.200 2009/01/19 02:38:17 rtoy Exp $"
7 ;;; "f2cl6.l,v 1.48 2008/08/24 00:56:27 rtoy Exp $"
8 ;;; "macros.l,v 1.112 2009/01/08 12:57:19 rtoy Exp $")
10 ;;; Using Lisp CMU Common Lisp 19f (19F)
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))
17 (in-package :common-lisp-user)
20 (let ((finish nil)
21 (iter 0)
22 (nfun 0)
23 (point 0)
24 (ispt 0)
25 (iypt 0)
26 (maxfev 0)
27 (info 0)
28 (bound 0)
29 (npt 0)
30 (cp 0)
31 (i 0)
32 (nfev 0)
33 (inmc 0)
34 (iycn 0)
35 (iscn 0)
36 (one 1.0)
37 (zero 0.0)
38 (gnorm 0.0)
39 (stp1 0.0)
40 (ftol 0.0)
41 (stp 0.0)
42 (ys 0.0)
43 (yy 0.0)
44 (sq 0.0)
45 (yr 0.0)
46 (beta 0.0)
47 (xnorm 0.0))
48 (declare (type f2cl-lib:logical finish)
49 (type (f2cl-lib:integer4) iter nfun point ispt iypt maxfev info
50 bound npt cp i nfev inmc iycn iscn)
51 (type (double-float) one zero gnorm stp1 ftol stp ys yy sq yr beta
52 xnorm))
53 (defun lbfgs (n m x f g diagco diag iprint eps xtol w iflag scache)
54 (declare (type (array f2cl-lib:integer4 (*)) iprint)
55 (type f2cl-lib:logical diagco)
56 (type (double-float) xtol eps f)
57 (type (array double-float (*)) scache w diag g x)
58 (type (f2cl-lib:integer4) iflag m n))
59 (let ()
60 (symbol-macrolet ((gtol (lb3-gtol *lb3-common-block*))
61 (lp (lb3-lp *lb3-common-block*)))
62 (f2cl-lib:with-multi-array-data
63 ((x double-float x-%data% x-%offset%)
64 (g double-float g-%data% g-%offset%)
65 (diag double-float diag-%data% diag-%offset%)
66 (w double-float w-%data% w-%offset%)
67 (scache double-float scache-%data% scache-%offset%)
68 (iprint f2cl-lib:integer4 iprint-%data% iprint-%offset%))
69 (prog ()
70 (declare)
71 (if (= iflag 0) (go label10))
72 (f2cl-lib:computed-goto (label172 label100) iflag)
73 label10
74 (setf iter 0)
75 (if (or (<= n 0) (<= m 0)) (go label196))
76 (cond
77 ((<= gtol 1.0e-4)
78 (if (> lp 0)
79 (f2cl-lib:fformat lp
80 ("~%"
81 " GTOL IS LESS THAN OR EQUAL TO 1.D-04"
82 "~%" " IT HAS BEEN RESET TO 9.D-01"
83 "~%")))
84 (setf gtol 0.9)))
85 (setf nfun 1)
86 (setf point 0)
87 (setf finish f2cl-lib:%false%)
88 (cond
89 (diagco
90 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
91 ((> i n) nil)
92 (tagbody
93 label30
94 (if
95 (<= (f2cl-lib:fref diag-%data% (i) ((1 n)) diag-%offset%)
96 zero)
97 (go label195)))))
99 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
100 ((> i n) nil)
101 (tagbody
102 label40
103 (setf (f2cl-lib:fref diag-%data% (i) ((1 n)) diag-%offset%)
104 1.0)))))
105 (setf ispt (f2cl-lib:int-add n (f2cl-lib:int-mul 2 m)))
106 (setf iypt (f2cl-lib:int-add ispt (f2cl-lib:int-mul n m)))
107 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
108 ((> i n) nil)
109 (tagbody
110 label50
111 (setf (f2cl-lib:fref w-%data%
112 ((f2cl-lib:int-add ispt i))
114 (f2cl-lib:int-add
115 (f2cl-lib:int-mul n
116 (f2cl-lib:int-add
117 (f2cl-lib:int-mul 2
120 (f2cl-lib:int-mul 2 m))))
121 w-%offset%)
122 (* (- (f2cl-lib:fref g-%data% (i) ((1 n)) g-%offset%))
123 (f2cl-lib:fref diag-%data%
125 ((1 n))
126 diag-%offset%)))))
127 (setf gnorm
128 (f2cl-lib:dsqrt
129 (multiple-value-bind
130 (ret-val var-0 var-1 var-2 var-3 var-4)
131 (ddot n g 1 g 1)
132 (declare (ignore var-1 var-2 var-3 var-4))
133 (when var-0
134 (setf n var-0))
135 ret-val)))
136 (setf stp1 (/ one gnorm))
137 (setf ftol 1.0e-4)
138 (setf maxfev 20)
140 (>= (f2cl-lib:fref iprint-%data% (1) ((1 2)) iprint-%offset%) 0)
141 (lb1 iprint iter nfun gnorm n m x f g stp finish))
142 label80
143 (setf iter (f2cl-lib:int-add iter 1))
144 (setf info 0)
145 (setf bound (f2cl-lib:int-sub iter 1))
146 (if (= iter 1) (go label165))
147 (if (> iter m) (setf bound m))
148 (setf ys
149 (multiple-value-bind
150 (ret-val var-0 var-1 var-2 var-3 var-4)
151 (ddot n
152 (f2cl-lib:array-slice w
153 double-float
154 ((+ iypt npt 1))
156 (f2cl-lib:int-add
157 (f2cl-lib:int-mul n
158 (f2cl-lib:int-add
159 (f2cl-lib:int-mul
163 (f2cl-lib:int-mul 2 m)))))
165 (f2cl-lib:array-slice w
166 double-float
167 ((+ ispt npt 1))
169 (f2cl-lib:int-add
170 (f2cl-lib:int-mul n
171 (f2cl-lib:int-add
172 (f2cl-lib:int-mul
176 (f2cl-lib:int-mul 2 m)))))
178 (declare (ignore var-1 var-2 var-3 var-4))
179 (when var-0
180 (setf n var-0))
181 ret-val))
182 (cond
183 ((not diagco)
184 (setf yy
185 (multiple-value-bind
186 (ret-val var-0 var-1 var-2 var-3 var-4)
187 (ddot n
188 (f2cl-lib:array-slice w
189 double-float
190 ((+ iypt npt 1))
192 (f2cl-lib:int-add
193 (f2cl-lib:int-mul n
194 (f2cl-lib:int-add
195 (f2cl-lib:int-mul
199 (f2cl-lib:int-mul 2 m)))))
201 (f2cl-lib:array-slice w
202 double-float
203 ((+ iypt npt 1))
205 (f2cl-lib:int-add
206 (f2cl-lib:int-mul n
207 (f2cl-lib:int-add
208 (f2cl-lib:int-mul
212 (f2cl-lib:int-mul 2 m)))))
214 (declare (ignore var-1 var-2 var-3 var-4))
215 (when var-0
216 (setf n var-0))
217 ret-val))
218 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
219 ((> i n) nil)
220 (tagbody
221 label90
222 (setf (f2cl-lib:fref diag-%data% (i) ((1 n)) diag-%offset%)
223 (/ ys yy)))))
225 (setf iflag 2)
226 (go end_label)))
227 label100
228 (cond
229 (diagco
230 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
231 ((> i n) nil)
232 (tagbody
233 label110
235 (<= (f2cl-lib:fref diag-%data% (i) ((1 n)) diag-%offset%)
236 zero)
237 (go label195))))))
238 (setf cp point)
239 (if (= point 0) (setf cp m))
240 (setf (f2cl-lib:fref w-%data%
241 ((f2cl-lib:int-add n cp))
243 (f2cl-lib:int-add
244 (f2cl-lib:int-mul n
245 (f2cl-lib:int-add
246 (f2cl-lib:int-mul 2 m)
248 (f2cl-lib:int-mul 2 m))))
249 w-%offset%)
250 (/ one ys))
251 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
252 ((> i n) nil)
253 (tagbody
254 label112
255 (setf (f2cl-lib:fref w-%data%
258 (f2cl-lib:int-add
259 (f2cl-lib:int-mul n
260 (f2cl-lib:int-add
261 (f2cl-lib:int-mul 2
264 (f2cl-lib:int-mul 2 m))))
265 w-%offset%)
266 (- (f2cl-lib:fref g-%data% (i) ((1 n)) g-%offset%)))))
267 (setf cp point)
268 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
269 ((> i bound) nil)
270 (tagbody
271 (setf cp (f2cl-lib:int-sub cp 1))
272 (if (= cp -1) (setf cp (f2cl-lib:int-sub m 1)))
273 (setf sq
274 (multiple-value-bind
275 (ret-val var-0 var-1 var-2 var-3 var-4)
276 (ddot n
277 (f2cl-lib:array-slice w
278 double-float
279 ((+ ispt
280 (f2cl-lib:int-mul cp n)
283 (f2cl-lib:int-add
284 (f2cl-lib:int-mul n
285 (f2cl-lib:int-add
286 (f2cl-lib:int-mul
290 (f2cl-lib:int-mul 2
291 m)))))
292 1 w 1)
293 (declare (ignore var-1 var-2 var-3 var-4))
294 (when var-0
295 (setf n var-0))
296 ret-val))
297 (setf inmc (f2cl-lib:int-add n m cp 1))
298 (setf iycn (f2cl-lib:int-add iypt (f2cl-lib:int-mul cp n)))
299 (setf (f2cl-lib:fref w-%data%
300 (inmc)
302 (f2cl-lib:int-add
303 (f2cl-lib:int-mul n
304 (f2cl-lib:int-add
305 (f2cl-lib:int-mul 2
308 (f2cl-lib:int-mul 2 m))))
309 w-%offset%)
311 (f2cl-lib:fref w-%data%
312 ((f2cl-lib:int-add n cp 1))
314 (f2cl-lib:int-add
315 (f2cl-lib:int-mul n
316 (f2cl-lib:int-add
317 (f2cl-lib:int-mul
321 (f2cl-lib:int-mul 2 m))))
322 w-%offset%)
323 sq))
324 (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
325 (daxpy n
327 (f2cl-lib:fref w-%data%
328 (inmc)
330 (f2cl-lib:int-add
331 (f2cl-lib:int-mul n
332 (f2cl-lib:int-add
333 (f2cl-lib:int-mul 2
336 (f2cl-lib:int-mul 2 m))))
337 w-%offset%))
338 (f2cl-lib:array-slice w
339 double-float
340 ((+ iycn 1))
342 (f2cl-lib:int-add
343 (f2cl-lib:int-mul n
344 (f2cl-lib:int-add
345 (f2cl-lib:int-mul
349 (f2cl-lib:int-mul 2 m)))))
350 1 w 1)
351 (declare (ignore var-1 var-2 var-3 var-4 var-5))
352 (when var-0
353 (setf n var-0)))
354 label125))
355 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
356 ((> i n) nil)
357 (tagbody
358 label130
359 (setf (f2cl-lib:fref w-%data%
362 (f2cl-lib:int-add
363 (f2cl-lib:int-mul n
364 (f2cl-lib:int-add
365 (f2cl-lib:int-mul 2
368 (f2cl-lib:int-mul 2 m))))
369 w-%offset%)
371 (f2cl-lib:fref diag-%data% (i) ((1 n)) diag-%offset%)
372 (f2cl-lib:fref w-%data%
375 (f2cl-lib:int-add
376 (f2cl-lib:int-mul n
377 (f2cl-lib:int-add
378 (f2cl-lib:int-mul
382 (f2cl-lib:int-mul 2 m))))
383 w-%offset%)))))
384 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
385 ((> i bound) nil)
386 (tagbody
387 (setf yr
388 (multiple-value-bind
389 (ret-val var-0 var-1 var-2 var-3 var-4)
390 (ddot n
391 (f2cl-lib:array-slice w
392 double-float
393 ((+ iypt
394 (f2cl-lib:int-mul cp n)
397 (f2cl-lib:int-add
398 (f2cl-lib:int-mul n
399 (f2cl-lib:int-add
400 (f2cl-lib:int-mul
404 (f2cl-lib:int-mul 2
405 m)))))
406 1 w 1)
407 (declare (ignore var-1 var-2 var-3 var-4))
408 (when var-0
409 (setf n var-0))
410 ret-val))
411 (setf beta
413 (f2cl-lib:fref w-%data%
414 ((f2cl-lib:int-add n cp 1))
416 (f2cl-lib:int-add
417 (f2cl-lib:int-mul n
418 (f2cl-lib:int-add
419 (f2cl-lib:int-mul
423 (f2cl-lib:int-mul 2 m))))
424 w-%offset%)
425 yr))
426 (setf inmc (f2cl-lib:int-add n m cp 1))
427 (setf beta
429 (f2cl-lib:fref w-%data%
430 (inmc)
432 (f2cl-lib:int-add
433 (f2cl-lib:int-mul n
434 (f2cl-lib:int-add
435 (f2cl-lib:int-mul
439 (f2cl-lib:int-mul 2 m))))
440 w-%offset%)
441 beta))
442 (setf iscn (f2cl-lib:int-add ispt (f2cl-lib:int-mul cp n)))
443 (multiple-value-bind (var-0 var-1 var-2 var-3 var-4 var-5)
444 (daxpy n beta
445 (f2cl-lib:array-slice w
446 double-float
447 ((+ iscn 1))
449 (f2cl-lib:int-add
450 (f2cl-lib:int-mul n
451 (f2cl-lib:int-add
452 (f2cl-lib:int-mul
456 (f2cl-lib:int-mul 2 m)))))
457 1 w 1)
458 (declare (ignore var-2 var-3 var-4 var-5))
459 (when var-0
460 (setf n var-0))
461 (when var-1
462 (setf beta var-1)))
463 (setf cp (f2cl-lib:int-add cp 1))
464 (if (= cp m) (setf cp 0))
465 label145))
466 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
467 ((> i n) nil)
468 (tagbody
469 label160
470 (setf (f2cl-lib:fref w-%data%
471 ((f2cl-lib:int-add ispt
472 (f2cl-lib:int-mul point
476 (f2cl-lib:int-add
477 (f2cl-lib:int-mul n
478 (f2cl-lib:int-add
479 (f2cl-lib:int-mul 2
482 (f2cl-lib:int-mul 2 m))))
483 w-%offset%)
484 (f2cl-lib:fref w-%data%
487 (f2cl-lib:int-add
488 (f2cl-lib:int-mul n
489 (f2cl-lib:int-add
490 (f2cl-lib:int-mul
494 (f2cl-lib:int-mul 2 m))))
495 w-%offset%))))
496 label165
497 (setf nfev 0)
498 (setf stp one)
499 (if (= iter 1) (setf stp stp1))
500 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
501 ((> i n) nil)
502 (tagbody
503 label170
504 (setf (f2cl-lib:fref w-%data%
507 (f2cl-lib:int-add
508 (f2cl-lib:int-mul n
509 (f2cl-lib:int-add
510 (f2cl-lib:int-mul 2
513 (f2cl-lib:int-mul 2 m))))
514 w-%offset%)
515 (f2cl-lib:fref g-%data% (i) ((1 n)) g-%offset%))))
516 label172
517 (multiple-value-bind
518 (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9
519 var-10 var-11)
520 (mcsrch n x f g
521 (f2cl-lib:array-slice w
522 double-float
523 ((+ ispt (f2cl-lib:int-mul point n) 1))
525 (f2cl-lib:int-add
526 (f2cl-lib:int-mul n
527 (f2cl-lib:int-add
528 (f2cl-lib:int-mul
532 (f2cl-lib:int-mul 2 m)))))
533 stp ftol xtol maxfev info nfev diag)
534 (declare (ignore var-0 var-1 var-2 var-3 var-4 var-6 var-7 var-8
535 var-11))
536 (setf stp var-5)
537 (setf info var-9)
538 (setf nfev var-10))
539 (cond
540 ((= info (f2cl-lib:int-sub 1))
541 (setf iflag 1)
542 (go end_label)))
543 (if (/= info 1) (go label190))
544 (setf nfun (f2cl-lib:int-add nfun nfev))
545 (setf npt (f2cl-lib:int-mul point n))
546 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
547 ((> i n) nil)
548 (tagbody
549 (setf (f2cl-lib:fref w-%data%
550 ((f2cl-lib:int-add ispt npt i))
552 (f2cl-lib:int-add
553 (f2cl-lib:int-mul n
554 (f2cl-lib:int-add
555 (f2cl-lib:int-mul 2
558 (f2cl-lib:int-mul 2 m))))
559 w-%offset%)
560 (* stp
561 (f2cl-lib:fref w-%data%
562 ((f2cl-lib:int-add ispt npt i))
564 (f2cl-lib:int-add
565 (f2cl-lib:int-mul n
566 (f2cl-lib:int-add
567 (f2cl-lib:int-mul
571 (f2cl-lib:int-mul 2 m))))
572 w-%offset%)))
573 label175
574 (setf (f2cl-lib:fref w-%data%
575 ((f2cl-lib:int-add iypt npt i))
577 (f2cl-lib:int-add
578 (f2cl-lib:int-mul n
579 (f2cl-lib:int-add
580 (f2cl-lib:int-mul 2
583 (f2cl-lib:int-mul 2 m))))
584 w-%offset%)
585 (- (f2cl-lib:fref g-%data% (i) ((1 n)) g-%offset%)
586 (f2cl-lib:fref w-%data%
589 (f2cl-lib:int-add
590 (f2cl-lib:int-mul n
591 (f2cl-lib:int-add
592 (f2cl-lib:int-mul
596 (f2cl-lib:int-mul 2 m))))
597 w-%offset%)))))
598 (setf point (f2cl-lib:int-add point 1))
599 (if (= point m) (setf point 0))
600 (setf gnorm
601 (f2cl-lib:dsqrt
602 (multiple-value-bind
603 (ret-val var-0 var-1 var-2 var-3 var-4)
604 (ddot n g 1 g 1)
605 (declare (ignore var-1 var-2 var-3 var-4))
606 (when var-0
607 (setf n var-0))
608 ret-val)))
609 (setf xnorm
610 (f2cl-lib:dsqrt
611 (multiple-value-bind
612 (ret-val var-0 var-1 var-2 var-3 var-4)
613 (ddot n x 1 x 1)
614 (declare (ignore var-1 var-2 var-3 var-4))
615 (when var-0
616 (setf n var-0))
617 ret-val)))
618 (setf xnorm (f2cl-lib:dmax1 1.0 xnorm))
619 (if (<= (/ gnorm xnorm) eps) (setf finish f2cl-lib:%true%))
621 (>= (f2cl-lib:fref iprint-%data% (1) ((1 2)) iprint-%offset%) 0)
622 (lb1 iprint iter nfun gnorm n m x f g stp finish))
623 (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
624 ((> i n) nil)
625 (tagbody
626 (setf (f2cl-lib:fref scache-%data% (i) ((1 n)) scache-%offset%)
627 (f2cl-lib:fref x-%data% (i) ((1 n)) x-%offset%))
628 label177))
629 (cond
630 (finish
631 (setf iflag 0)
632 (go end_label)))
633 (go label80)
634 label190
635 (setf iflag -1)
636 (if (> lp 0)
637 (f2cl-lib:fformat lp
638 ("~%" " IFLAG= -1 " "~%"
639 " LINE SEARCH FAILED. SEE"
640 " DOCUMENTATION OF ROUTINE MCSRCH" "~%"
641 " ERROR RETURN" " OF LINE SEARCH: INFO= " 1
642 (("~2D")) "~%"
643 " POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT"
644 "~%" " OR INCORRECT TOLERANCES" "~%")
645 info))
646 (go end_label)
647 label195
648 (setf iflag -2)
649 (if (> lp 0)
650 (f2cl-lib:fformat lp
651 ("~%" " IFLAG= -2" "~%" " THE" 1 (("~5D"))
652 "-TH DIAGONAL ELEMENT OF THE" "~%"
653 " INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE"
654 "~%")
656 (go end_label)
657 label196
658 (setf iflag -3)
659 (if (> lp 0)
660 (f2cl-lib:fformat lp
661 ("~%" " IFLAG= -3" "~%"
662 " IMPROPER INPUT PARAMETERS (N OR M"
663 " ARE NOT POSITIVE)" "~%")))
664 (go end_label)
665 end_label
666 (return
667 (values n nil nil nil nil nil nil nil nil nil nil iflag nil))))))))
669 (in-package #:cl-user)
670 #+#.(cl:if (cl:find-package '#:f2cl) '(and) '(or))
671 (eval-when (:load-toplevel :compile-toplevel :execute)
672 (setf (gethash 'fortran-to-lisp::lbfgs fortran-to-lisp::*f2cl-function-info*)
673 (fortran-to-lisp::make-f2cl-finfo
674 :arg-types '((fortran-to-lisp::integer4) (fortran-to-lisp::integer4)
675 (array double-float (*)) (double-float)
676 (array double-float (*)) fortran-to-lisp::logical
677 (array double-float (*))
678 (array fortran-to-lisp::integer4 (2)) (double-float)
679 (double-float) (array double-float (*))
680 (fortran-to-lisp::integer4) (array double-float (*)))
681 :return-values '(fortran-to-lisp::n nil nil nil nil nil nil nil nil
682 nil nil fortran-to-lisp::iflag nil)
683 :calls '(fortran-to-lisp::mcsrch fortran-to-lisp::lb1))))