#include <cstdio>#include <cstring>#include <cmath>#include <cstdlib>Include dependency graph for lbfgsb.cpp:

Go to the source code of this file.
Classes | |
| struct | cilist |
| struct | olist |
Defines | |
| #define | TRUE_ (1) |
| #define | FALSE_ (0) |
| #define | min_(a, b) ((a) <= (b) ? (a) : (b)) |
| #define | max_(a, b) ((a) >= (b) ? (a) : (b)) |
| #define | abs_(x) ((x) >= 0 ? (x) : -(x)) |
Functions | |
| long int | s_cmp (char *str1, const char *const str2, long int l1, long int l2) |
| int | s_copy (char *str1, const char *const str2, long int l1, long int l2) |
| int | setulb_ (long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, double *g, double *factr, double *pgtol, double *wa, long int *iwa, char *task, long int *iprint, char *csave, long int *lsave, long int *isave, double *dsave, long int task_len, long int csave_len) |
| int | mainlb_ (long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f, double *g, double *factr, double *pgtol, double *ws, double *wy, double *sy, double *ss, double *yy, double *wt, double *wn, double *snd, double *z__, double *r__, double *d__, double *t, double *wa, double *sg, double *sgo, double *yg, double *ygo, long int *index, long int *iwhere, long int *indx2, char *task, long int *iprint, char *csave, long int *lsave, long int *isave, double *dsave, long int task_len, long int csave_len) |
| int | active_ (long int *n, double *l, double *u, long int *nbd, double *x, long int *iwhere, long int *iprint, long int *prjctd, long int *cnstnd, long int *boxed) |
| int | bmv_ (long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info) |
| int | cauchy_ (long int *n, double *x, double *l, double *u, long int *nbd, double *g, long int *iorder, long int *iwhere, double *t, double *d__, double *xcp, long int *m, double *wy, double *ws, double *sy, double *wt, double *theta, long int *col, long int *head, double *p, double *c__, double *wbp, double *v, long int *nint, double *sg, double *yg, long int *iprint, double *sbgnrm, long int *info) |
| int | cmprlb_ (long int *n, long int *m, double *x, double *g, double *ws, double *wy, double *sy, double *wt, double *z__, double *r__, double *wa, long int *index, double *theta, long int *col, long int *head, long int *nfree, long int *cnstnd, long int *info) |
| int | errclb_ (long int *n, long int *m, double *factr, double *l, double *u, long int *nbd, char *task, long int *info, long int *k, long int task_len) |
| int | formk_ (long int *n, long int *nsub, long int *ind, long int *nenter, long int *ileave, long int *indx2, long int *iupdat, long int *updatd, double *wn, double *wn1, long int *m, double *ws, double *wy, double *sy, double *theta, long int *col, long int *head, long int *info) |
| int | formt_ (long int *m, double *wt, double *sy, double *ss, long int *col, double *theta, long int *info) |
| int | freev_ (long int *n, long int *nfree, long int *index, long int *nenter, long int *ileave, long int *indx2, long int *iwhere, long int *wrk, long int *updatd, long int *cnstnd, long int *iprint, long int *iter) |
| int | hpsolb_ (long int *n, double *t, long int *iorder, long int *iheap) |
| int | lnsrlb_ (long int *n, double *l, double *u, long int *nbd, double *x, double *f, double *fold, double *gd, double *gdold, double *g, double *d__, double *r__, double *t, double *z__, double *stp, double *dnorm, double *dtd, double *xstep, double *stpmx, long int *iter, long int *ifun, long int *iback, long int *nfgv, long int *info, char *task, long int *boxed, long int *cnstnd, char *csave, long int *isave, double *dsave, long int task_len, long int csave_len) |
| int | matupd_ (long int *n, long int *m, double *ws, double *wy, double *sy, double *ss, double *d__, double *r__, long int *itail, long int *iupdat, long int *col, long int *head, double *theta, double *rr, double *dr, double *stp, double *dtd) |
| int | prn1lb_ (long int *n, long int *m, double *l, double *u, double *x, long int *iprint, long int *itfile, double *epsmch) |
| int | prn2lb_ (long int *n, double *x, double *f, double *g, long int *iprint, long int *itfile, long int *iter, long int *nfgv, long int *nact, double *sbgnrm, long int *nint, char *word, long int *iword, long int *iback, double *stp, double *xstep, long int word_len) |
| int | prn3lb_ (long int *n, double *x, double *f, char *task, long int *iprint, long int *info, long int *itfile, long int *iter, long int *nfgv, long int *nintol, long int *nskip, long int *nact, double *sbgnrm, double *time, long int *nint, char *word, long int *iback, double *stp, double *xstep, long int *k, double *cachyt, double *sbtime, double *lnscht, long int task_len, long int word_len) |
| int | projgr_ (long int *n, double *l, double *u, long int *nbd, double *x, double *g, double *sbgnrm) |
| int | subsm_ (long int *n, long int *m, long int *nsub, long int *ind, double *l, double *u, long int *nbd, double *x, double *d__, double *ws, double *wy, double *theta, long int *col, long int *head, long int *iword, double *wv, double *wn, long int *iprint, long int *info) |
| int | dcsrch_ (double *f, double *g, double *stp, double *ftol, double *gtol, double *xtol, double *stpmin, double *stpmax, char *task, long int *isave, double *dsave, long int task_len) |
| int | dcstep_ (double *stx, double *fx, double *dx, double *sty, double *fy, double *dy, double *stp, double *fp, double *dp, long int *brackt, double *stpmin, double *stpmax) |
| int | timer_ (double *ttime) |
| double | dnrm2_ (long int *n, double *x, long int *incx) |
| double | dpmeps_ () |
| int | daxpy_ (long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy) |
| int | dcopy_ (long int *n, double *dx, long int *incx, double *dy, long int *incy) |
| double | ddot_ (long int *n, double *dx, long int *incx, double *dy, long int *incy) |
| int | dpofa_ (double *a, long int *lda, long int *n, long int *info) |
| int | dscal_ (long int *n, double *da, double *dx, long int *incx) |
| int | dtrsl_ (double *t, long int *ldt, long int *n, double *b, long int *job, long int *info) |
Variables | |
| double | c_b9 = 0. |
| long int | c__1 = 1 |
| long int | c__11 = 11 |
| double | c_b275 = .001 |
| double | c_b276 = .9 |
| double | c_b277 = .1 |
| const int | SIXTY = 60 |
|
|
Definition at line 10 of file lbfgsb.cpp. Referenced by cauchy_(), dcsrch_(), dcstep_(), dnrm2_(), mainlb_(), and projgr_(). |
|
|
Definition at line 7 of file lbfgsb.cpp. |
|
|
Definition at line 9 of file lbfgsb.cpp. Referenced by dcsrch_(), dcstep_(), dnrm2_(), mainlb_(), and projgr_(). |
|
|
Definition at line 8 of file lbfgsb.cpp. Referenced by dcsrch_(), dcstep_(), formt_(), lnsrlb_(), projgr_(), s_cmp(), and s_copy(). |
|
|
Definition at line 6 of file lbfgsb.cpp. |
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 1497 of file lbfgsb.cpp. References x. Referenced by mainlb_(). 01505 {
01506 /* Format strings *//*
01507 static char fmt_1001[] = "(/,\002At X0 \002,i9,\002 variables are exactl\
01508 y at the bounds\002)";*/
01509
01510 /* System generated locals */
01511 long int i__1;
01512
01513 /* Builtin functions */
01514 //long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
01515
01516 /* Local variables */
01517 static long int nbdd, i__;
01518
01519 /* Fortran I/O blocks */
01520 /* static cilist io___81 = { 0, 6, 0, 0, 0 };
01521 static cilist io___82 = { 0, 6, 0, 0, 0 };
01522 static cilist io___83 = { 0, 6, 0, fmt_1001, 0 };*/
01523
01524
01525 /* ************ */
01526
01527 /* Subroutine active */
01528
01529 /* This subroutine initializes iwhere and projects the initial x to */
01530 /* the feasible set if necessary. */
01531
01532 /* iwhere is an long int array of dimension n. */
01533 /* On entry iwhere is unspecified. */
01534 /* On exit iwhere(i)=-1 if x(i) has no bounds */
01535 /* 3 if l(i)=u(i) */
01536 /* 0 otherwise. */
01537 /* In cauchy, iwhere is given finer gradations. */
01538
01539
01540 /* * * * */
01541
01542 /* NEOS, November 1994. (Latest revision June 1996.) */
01543 /* Optimization Technology Center. */
01544 /* Argonne National Laboratory and Northwestern University. */
01545 /* Written by */
01546 /* Ciyou Zhu */
01547 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
01548
01549
01550 /* ************ */
01551 /* Initialize nbdd, prjctd, cnstnd and boxed. */
01552 /* Parameter adjustments */
01553 --iwhere;
01554 --x;
01555 --nbd;
01556 --u;
01557 --l;
01558
01559 /* Function Body */
01560 nbdd = 0;
01561 *prjctd = FALSE_;
01562 *cnstnd = FALSE_;
01563 *boxed = TRUE_;
01564 /* Project the initial x to the easible set if necessary. */
01565 i__1 = *n;
01566 for (i__ = 1; i__ <= i__1; ++i__) {
01567 if (nbd[i__] > 0) {
01568 if (nbd[i__] <= 2 && x[i__] <= l[i__]) {
01569 if (x[i__] < l[i__]) {
01570 *prjctd = TRUE_;
01571 x[i__] = l[i__];
01572 }
01573 ++nbdd;
01574 } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) {
01575 if (x[i__] > u[i__]) {
01576 *prjctd = TRUE_;
01577 x[i__] = u[i__];
01578 }
01579 ++nbdd;
01580 }
01581 }
01582 /* L10: */
01583 }
01584 /* Initialize iwhere and assign values to cnstnd and boxed. */
01585 i__1 = *n;
01586 for (i__ = 1; i__ <= i__1; ++i__) {
01587 if (nbd[i__] != 2) {
01588 *boxed = FALSE_;
01589 }
01590 if (nbd[i__] == 0) {
01591 /* this variable is always free */
01592 iwhere[i__] = -1;
01593 /* otherwise set x(i)=mid(x(i), u(i), l(i)). */
01594 } else {
01595 *cnstnd = TRUE_;
01596 if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) {
01597 /* this variable is always fixed */
01598 iwhere[i__] = 3;
01599 } else {
01600 iwhere[i__] = 0;
01601 }
01602 }
01603 /* L20: */
01604 }
01605 if (*iprint >= 0) {
01606 /* if (*prjctd) {
01607 s_wsle(&io___81);
01608 do_lio(&c__9, &c__1, "The initial X is infeasible. Restart with\
01609 its projection.", (long int)58);
01610 e_wsle();
01611 }
01612 if (! (*cnstnd)) {
01613 s_wsle(&io___82);
01614 do_lio(&c__9, &c__1, "This problem is unconstrained.", (long int)30)
01615 ;
01616 e_wsle();
01617 }*/
01618 }
01619 if (*iprint > 0) {
01620 /* s_wsfe(&io___83);
01621 do_fio(&c__1, (char *)&nbdd, (long int)sizeof(long int));
01622 e_wsfe();*/
01623 }
01624 return 0;
01625 } /* active_ */
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 1628 of file lbfgsb.cpp. References b, c__1, c__11, dtrsl_(), sqrt(), t, and v. Referenced by cauchy_(), and cmprlb_(). 01634 {
01635 /* System generated locals */
01636 long int sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2;
01637
01638 /* Builtin functions */
01639 //double sqrt();
01640
01641 /* Local variables */
01642 static long int i__, k;
01643 extern /* Subroutine */ int dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
01644 static long int i2;
01645 static double sum;
01646
01647 /* ************ */
01648
01649 /* Subroutine bmv */
01650
01651 /* This subroutine computes the product of the 2m x 2m middle matrix */
01652 /* in the compact L-BFGS formula of B and a 2m vector v; */
01653 /* it returns the product in p. */
01654
01655 /* m is an long int variable. */
01656 /* On entry m is the maximum number of variable metric corrections */
01657 /* used to define the limited memory matrix. */
01658 /* On exit m is unchanged. */
01659
01660 /* sy is a double precision array of dimension m x m. */
01661 /* On entry sy specifies the matrix S'Y. */
01662 /* On exit sy is unchanged. */
01663
01664 /* wt is a double precision array of dimension m x m. */
01665 /* On entry wt specifies the upper triangular matrix J' which is */
01666 /* the Cholesky factor of (thetaS'S+LD^(-1)L'). */
01667 /* On exit wt is unchanged. */
01668
01669 /* col is an long int variable. */
01670 /* On entry col specifies the number of s-vectors (or y-vectors) */
01671 /* stored in the compact L-BFGS formula. */
01672 /* On exit col is unchanged. */
01673
01674 /* v is a double precision array of dimension 2col. */
01675 /* On entry v specifies vector v. */
01676 /* On exit v is unchanged. */
01677
01678 /* p is a double precision array of dimension 2col. */
01679 /* On entry p is unspecified. */
01680 /* On exit p is the product Mv. */
01681
01682 /* info is an long int variable. */
01683 /* On entry info is unspecified. */
01684 /* On exit info = 0 for normal return, */
01685 /* = nonzero for abnormal return when the system */
01686 /* to be solved by dtrsl is singular. */
01687
01688 /* Subprograms called: */
01689
01690 /* Linpack ... dtrsl. */
01691
01692
01693 /* * * * */
01694
01695 /* NEOS, November 1994. (Latest revision June 1996.) */
01696 /* Optimization Technology Center. */
01697 /* Argonne National Laboratory and Northwestern University. */
01698 /* Written by */
01699 /* Ciyou Zhu */
01700 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
01701
01702
01703 /* ************ */
01704 /* Parameter adjustments */
01705 wt_dim1 = *m;
01706 wt_offset = 1 + wt_dim1 * 1;
01707 wt -= wt_offset;
01708 sy_dim1 = *m;
01709 sy_offset = 1 + sy_dim1 * 1;
01710 sy -= sy_offset;
01711 --p;
01712 --v;
01713
01714 /* Function Body */
01715 if (*col == 0) {
01716 return 0;
01717 }
01718 /* PART I: solve [ D^(1/2) O ] [ p1 ] = [ v1 ] */
01719 /* [ -L*D^(-1/2) J ] [ p2 ] [ v2 ]. */
01720 /* solve Jp2=v2+LD^(-1)v1. */
01721 p[*col + 1] = v[*col + 1];
01722 i__1 = *col;
01723 for (i__ = 2; i__ <= i__1; ++i__) {
01724 i2 = *col + i__;
01725 sum = 0.;
01726 i__2 = i__ - 1;
01727 for (k = 1; k <= i__2; ++k) {
01728 sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1];
01729 /* L10: */
01730 }
01731 p[i2] = v[i2] + sum;
01732 /* L20: */
01733 }
01734 /* Solve the triangular system */
01735 dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info);
01736 if (*info != 0) {
01737 return 0;
01738 }
01739 /* solve D^(1/2)p1=v1. */
01740 i__1 = *col;
01741 for (i__ = 1; i__ <= i__1; ++i__) {
01742 p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
01743 /* L30: */
01744 }
01745 /* PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] */
01746 /* [ 0 J' ] [ p2 ] [ p2 ]. */
01747 /* solve J^Tp2=p2. */
01748 dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info);
01749 if (*info != 0) {
01750 return 0;
01751 }
01752 /* compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */
01753 /* =-D^(-1/2)p1+D^(-1)L'p2. */
01754 i__1 = *col;
01755 for (i__ = 1; i__ <= i__1; ++i__) {
01756 p[i__] = -p[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
01757 /* L40: */
01758 }
01759 i__1 = *col;
01760 for (i__ = 1; i__ <= i__1; ++i__) {
01761 sum = 0.;
01762 i__2 = *col;
01763 for (k = i__ + 1; k <= i__2; ++k) {
01764 sum += sy[k + i__ * sy_dim1] * p[*col + k] / sy[i__ + i__ *
01765 sy_dim1];
01766 /* L50: */
01767 }
01768 p[i__] += sum;
01769 /* L60: */
01770 }
01771 return 0;
01772 } /* bmv_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 1775 of file lbfgsb.cpp. References abs_, bmv_(), c__1, daxpy_(), dcopy_(), ddot_(), dscal_(), dt, hpsolb_(), t, theta, v, and x. Referenced by mainlb_(). 01794 {
01795 /* Format strings *//*
01796 static char fmt_3010[] = "(/,\002---------------- CAUCHY entered--------\
01797 -----------\002)";
01798 static char fmt_1010[] = "(\002Cauchy X = \002,/,(4x,1p,6(1x,d11.4)))";
01799 static char fmt_4011[] = "(/,\002Piece \002,i3,\002 --f1, f2 at start\
01800 point \002,1p,2(1x,d11.4))";
01801 static char fmt_5010[] = "(\002Distance to the next break point = \002,\
01802 1p,d11.4)";
01803 static char fmt_6010[] = "(\002Distance to the stationary point = \002,\
01804 1p,d11.4)";
01805 static char fmt_4010[] = "(\002Piece \002,i3,\002 --f1, f2 at start p\
01806 oint \002,1p,2(1x,d11.4))";
01807 static char fmt_2010[] = "(/,\002---------------- exit CAUCHY-----------\
01808 -----------\002,/)";*/
01809
01810 /* System generated locals */
01811 long int wy_dim1, wy_offset, ws_dim1, ws_offset, sy_dim1, sy_offset,
01812 wt_dim1, wt_offset, i__1, i__2;
01813 double d__1;
01814
01815 /* Builtin functions */
01816 // long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int);
01817
01818 /* Local variables */
01819 static double dibp;
01820 extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
01821 static long int iter;
01822 static double zibp, tsum, dibp2;
01823 static long int i__, j;
01824 static long int bnded;
01825 extern /* Subroutine */ int dscal_(long int *n, double *da, double *dx, long int *incx);
01826 static double neggi;
01827 static long int nfree;
01828 static double bkmin;
01829 static long int nleft;
01830 extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy),
01831 daxpy_(long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy);
01832 static double f1, f2, dt, tj, tl;
01833 static long int nbreak, ibkmin;
01834 static double tu;
01835 extern /* Subroutine */ int hpsolb_(long int *n, double *t, long int *iorder, long int *iheap);
01836 static long int pointr;
01837 static double tj0;
01838 static long int xlower, xupper;
01839 static long int ibp;
01840 static double dtm;
01841 extern /* Subroutine */ int bmv_(long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info);
01842 static double wmc, wmp, wmw;
01843 static long int col2;
01844
01845 /* Fortran I/O blocks */
01846 /* static cilist io___88 = { 0, 6, 0, 0, 0 };
01847 static cilist io___96 = { 0, 6, 0, fmt_3010, 0 };
01848 static cilist io___105 = { 0, 6, 0, fmt_1010, 0 };
01849 static cilist io___109 = { 0, 6, 0, 0, 0 };
01850 static cilist io___116 = { 0, 6, 0, fmt_4011, 0 };
01851 static cilist io___117 = { 0, 6, 0, fmt_5010, 0 };
01852 static cilist io___118 = { 0, 6, 0, fmt_6010, 0 };
01853 static cilist io___121 = { 0, 6, 0, 0, 0 };
01854 static cilist io___126 = { 0, 6, 0, 0, 0 };
01855 static cilist io___127 = { 0, 6, 0, 0, 0 };
01856 static cilist io___128 = { 0, 6, 0, fmt_4010, 0 };
01857 static cilist io___129 = { 0, 6, 0, fmt_6010, 0 };
01858 static cilist io___130 = { 0, 6, 0, fmt_1010, 0 };
01859 static cilist io___131 = { 0, 6, 0, fmt_2010, 0 };*/
01860
01861
01862 /* ************ */
01863
01864 /* Subroutine cauchy */
01865
01866 /* For given x, l, u, g (with sbgnrm > 0), and a limited memory */
01867 /* BFGS matrix B defined in terms of matrices WY, WS, WT, and */
01868 /* scalars head, col, and theta, this subroutine computes the */
01869 /* generalized Cauchy point (GCP), defined as the first local */
01870 /* minimizer of the quadratic */
01871
01872 /* Q(x + s) = g's + 1/2 s'Bs */
01873
01874 /* along the projected gradient direction P(x-tg,l,u). */
01875 /* The routine returns the GCP in xcp. */
01876
01877 /* n is an long int variable. */
01878 /* On entry n is the dimension of the problem. */
01879 /* On exit n is unchanged. */
01880
01881 /* x is a double precision array of dimension n. */
01882 /* On entry x is the starting point for the GCP computation. */
01883 /* On exit x is unchanged. */
01884
01885 /* l is a double precision array of dimension n. */
01886 /* On entry l is the lower bound of x. */
01887 /* On exit l is unchanged. */
01888
01889 /* u is a double precision array of dimension n. */
01890 /* On entry u is the upper bound of x. */
01891 /* On exit u is unchanged. */
01892
01893 /* nbd is an long int array of dimension n. */
01894 /* On entry nbd represents the type of bounds imposed on the */
01895 /* variables, and must be specified as follows: */
01896 /* nbd(i)=0 if x(i) is unbounded, */
01897 /* 1 if x(i) has only a lower bound, */
01898 /* 2 if x(i) has both lower and upper bounds, and */
01899 /* 3 if x(i) has only an upper bound. */
01900 /* On exit nbd is unchanged. */
01901
01902 /* g is a double precision array of dimension n. */
01903 /* On entry g is the gradient of f(x). g must be a nonzero vector. */
01904 /* On exit g is unchanged. */
01905
01906 /* iorder is an long int working array of dimension n. */
01907 /* iorder will be used to store the breakpoints in the piecewise */
01908 /* linear path and free variables encountered. On exit, */
01909 /* iorder(1),...,iorder(nleft) are indices of breakpoints */
01910 /* which have not been encountered; */
01911 /* iorder(nleft+1),...,iorder(nbreak) are indices of */
01912 /* encountered breakpoints; and */
01913 /* iorder(nfree),...,iorder(n) are indices of variables which */
01914 /* have no bound constraits along the search direction. */
01915
01916 /* iwhere is an long int array of dimension n. */
01917 /* On entry iwhere indicates only the permanently fixed (iwhere=3) */
01918 /* or free (iwhere= -1) components of x. */
01919 /* On exit iwhere records the status of the current x variables. */
01920 /* iwhere(i)=-3 if x(i) is free and has bounds, but is not moved */
01921 /* 0 if x(i) is free and has bounds, and is moved */
01922 /* 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) */
01923 /* 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) */
01924 /* 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) */
01925 /* -1 if x(i) is always free, i.e., it has no bounds. */
01926
01927 /* t is a double precision working array of dimension n. */
01928 /* t will be used to store the break points. */
01929
01930 /* d is a double precision array of dimension n used to store */
01931 /* the Cauchy direction P(x-tg)-x. */
01932
01933 /* xcp is a double precision array of dimension n used to return the */
01934 /* GCP on exit. */
01935
01936 /* m is an long int variable. */
01937 /* On entry m is the maximum number of variable metric corrections */
01938 /* used to define the limited memory matrix. */
01939 /* On exit m is unchanged. */
01940
01941 /* ws, wy, sy, and wt are double precision arrays. */
01942 /* On entry they store information that defines the */
01943 /* limited memory BFGS matrix: */
01944 /* ws(n,m) stores S, a set of s-vectors; */
01945 /* wy(n,m) stores Y, a set of y-vectors; */
01946 /* sy(m,m) stores S'Y; */
01947 /* wt(m,m) stores the */
01948 /* Cholesky factorization of (theta*S'S+LD^(-1)L'). */
01949 /* On exit these arrays are unchanged. */
01950
01951 /* theta is a double precision variable. */
01952 /* On entry theta is the scaling factor specifying B_0 = theta I. */
01953 /* On exit theta is unchanged. */
01954
01955 /* col is an long int variable. */
01956 /* On entry col is the actual number of variable metric */
01957 /* corrections stored so far. */
01958 /* On exit col is unchanged. */
01959
01960 /* head is an long int variable. */
01961 /* On entry head is the location of the first s-vector (or y-vector) */
01962 /* in S (or Y). */
01963 /* On exit col is unchanged. */
01964
01965 /* p is a double precision working array of dimension 2m. */
01966 /* p will be used to store the vector p = W^(T)d. */
01967
01968 /* c is a double precision working array of dimension 2m. */
01969 /* c will be used to store the vector c = W^(T)(xcp-x). */
01970
01971 /* wbp is a double precision working array of dimension 2m. */
01972 /* wbp will be used to store the row of W corresponding */
01973 /* to a breakpoint. */
01974
01975 /* v is a double precision working array of dimension 2m. */
01976
01977 /* nint is an long int variable. */
01978 /* On exit nint records the number of quadratic segments explored */
01979 /* in searching for the GCP. */
01980
01981 /* sg and yg are double precision arrays of dimension m. */
01982 /* On entry sg and yg store S'g and Y'g correspondingly. */
01983 /* On exit they are unchanged. */
01984
01985 /* iprint is an long int variable that must be set by the user. */
01986 /* It controls the frequency and type of output generated: */
01987 /* iprint<0 no output is generated; */
01988 /* iprint=0 print only one line at the last iteration; */
01989 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
01990 /* iprint=99 print details of every iteration except n-vectors; */
01991 /* iprint=100 print also the changes of active set and final x; */
01992 /* iprint>100 print details of every iteration including x and g; */
01993 /* When iprint > 0, the file iterate.dat will be created to */
01994 /* summarize the iteration. */
01995
01996 /* sbgnrm is a double precision variable. */
01997 /* On entry sbgnrm is the norm of the projected gradient at x. */
01998 /* On exit sbgnrm is unchanged. */
01999
02000 /* info is an long int variable. */
02001 /* On entry info is 0. */
02002 /* On exit info = 0 for normal return, */
02003 /* = nonzero for abnormal return when the the system */
02004 /* used in routine bmv is singular. */
02005
02006 /* Subprograms called: */
02007
02008 /* L-BFGS-B Library ... hpsolb, bmv. */
02009
02010 /* Linpack ... dscal dcopy, daxpy. */
02011
02012
02013 /* References: */
02014
02015 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
02016 /* memory algorithm for bound constrained optimization'', */
02017 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
02018
02019 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
02020 /* Subroutines for Large Scale Bound Constrained Optimization'' */
02021 /* Tech. Report, NAM-11, EECS Department, Northwestern University, */
02022 /* 1994. */
02023
02024 /* (Postscript files of these papers are available via anonymous */
02025 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
02026
02027 /* * * * */
02028
02029 /* NEOS, November 1994. (Latest revision June 1996.) */
02030 /* Optimization Technology Center. */
02031 /* Argonne National Laboratory and Northwestern University. */
02032 /* Written by */
02033 /* Ciyou Zhu */
02034 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02035
02036
02037 /* ************ */
02038 /* Check the status of the variables, reset iwhere(i) if necessary; */
02039 /* compute the Cauchy direction d and the breakpoints t; initialize */
02040 /* the derivative f1 and the vector p = W'd (for theta = 1). */
02041 /* Parameter adjustments */
02042 --xcp;
02043 --d__;
02044 --t;
02045 --iwhere;
02046 --iorder;
02047 --g;
02048 --nbd;
02049 --u;
02050 --l;
02051 --x;
02052 --yg;
02053 --sg;
02054 --v;
02055 --wbp;
02056 --c__;
02057 --p;
02058 wt_dim1 = *m;
02059 wt_offset = 1 + wt_dim1 * 1;
02060 wt -= wt_offset;
02061 sy_dim1 = *m;
02062 sy_offset = 1 + sy_dim1 * 1;
02063 sy -= sy_offset;
02064 ws_dim1 = *n;
02065 ws_offset = 1 + ws_dim1 * 1;
02066 ws -= ws_offset;
02067 wy_dim1 = *n;
02068 wy_offset = 1 + wy_dim1 * 1;
02069 wy -= wy_offset;
02070
02071 /* Function Body */
02072 if (*sbgnrm <= 0.) {
02073 if (*iprint >= 0) {
02074 /* s_wsle(&io___88);
02075 do_lio(&c__9, &c__1, "Subgnorm = 0. GCP = X.", (long int)23);
02076 e_wsle();*/
02077 }
02078 dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
02079 return 0;
02080 }
02081 bnded = TRUE_;
02082 nfree = *n + 1;
02083 nbreak = 0;
02084 ibkmin = 0;
02085 bkmin = 0.;
02086 col2 = *col << 1;
02087 f1 = 0.;
02088 if (*iprint >= 99) {
02089 /* s_wsfe(&io___96);
02090 e_wsfe();*/
02091 }
02092 /* We set p to zero and build it up as we determine d. */
02093 i__1 = col2;
02094 for (i__ = 1; i__ <= i__1; ++i__) {
02095 p[i__] = 0.;
02096 /* L20: */
02097 }
02098 /* In the following loop we determine for each variable its bound */
02099 /* status and its breakpoint, and update p accordingly. */
02100 /* Smallest breakpoint is identified. */
02101 i__1 = *n;
02102 for (i__ = 1; i__ <= i__1; ++i__) {
02103 neggi = -g[i__];
02104 if (iwhere[i__] != 3 && iwhere[i__] != -1) {
02105 /* if x(i) is not a constant and has bounds, */
02106 /* compute the difference between x(i) and its bounds. */
02107 if (nbd[i__] <= 2) {
02108 tl = x[i__] - l[i__];
02109 }
02110 if (nbd[i__] >= 2) {
02111 tu = u[i__] - x[i__];
02112 }
02113 /* If a variable is close enough to a bound */
02114 /* we treat it as at bound. */
02115 xlower = nbd[i__] <= 2 && tl <= 0.;
02116 xupper = nbd[i__] >= 2 && tu <= 0.;
02117 /* reset iwhere(i). */
02118 iwhere[i__] = 0;
02119 if (xlower) {
02120 if (neggi <= 0.) {
02121 iwhere[i__] = 1;
02122 }
02123 } else if (xupper) {
02124 if (neggi >= 0.) {
02125 iwhere[i__] = 2;
02126 }
02127 } else {
02128 if (abs_(neggi) <= 0.) {
02129 iwhere[i__] = -3;
02130 }
02131 }
02132 }
02133 pointr = *head;
02134 if (iwhere[i__] != 0 && iwhere[i__] != -1) {
02135 d__[i__] = 0.;
02136 } else {
02137 d__[i__] = neggi;
02138 f1 -= neggi * neggi;
02139 /* calculate p := p - W'e_i* (g_i). */
02140 i__2 = *col;
02141 for (j = 1; j <= i__2; ++j) {
02142 p[j] += wy[i__ + pointr * wy_dim1] * neggi;
02143 p[*col + j] += ws[i__ + pointr * ws_dim1] * neggi;
02144 pointr = pointr % *m + 1;
02145 /* L40: */
02146 }
02147 if (nbd[i__] <= 2 && nbd[i__] != 0 && neggi < 0.) {
02148 /* x(i) + d(i) is bounded; compute t(i). */
02149 ++nbreak;
02150 iorder[nbreak] = i__;
02151 t[nbreak] = tl / (-neggi);
02152 if (nbreak == 1 || t[nbreak] < bkmin) {
02153 bkmin = t[nbreak];
02154 ibkmin = nbreak;
02155 }
02156 } else if (nbd[i__] >= 2 && neggi > 0.) {
02157 /* x(i) + d(i) is bounded; compute t(i). */
02158 ++nbreak;
02159 iorder[nbreak] = i__;
02160 t[nbreak] = tu / neggi;
02161 if (nbreak == 1 || t[nbreak] < bkmin) {
02162 bkmin = t[nbreak];
02163 ibkmin = nbreak;
02164 }
02165 } else {
02166 /* x(i) + d(i) is not bounded. */
02167 --nfree;
02168 iorder[nfree] = i__;
02169 if (abs_(neggi) > 0.) {
02170 bnded = FALSE_;
02171 }
02172 }
02173 }
02174 /* L50: */
02175 }
02176 /* The indices of the nonzero components of d are now stored */
02177 /* in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */
02178 /* The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */
02179 if (*theta != 1.) {
02180 /* complete the initialization of p for theta not= one. */
02181 dscal_(col, theta, &p[*col + 1], &c__1);
02182 }
02183 /* Initialize GCP xcp = x. */
02184 dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
02185 if (nbreak == 0 && nfree == *n + 1) {
02186 /* is a zero vector, return with the initial xcp as GCP. */
02187 if (*iprint > 100) {
02188 /* s_wsfe(&io___105);
02189 i__1 = *n;
02190 for (i__ = 1; i__ <= i__1; ++i__) {
02191 do_fio(&c__1, (char *)&xcp[i__], (long int)sizeof(double));
02192 }
02193 e_wsfe();*/
02194 }
02195 return 0;
02196 }
02197 /* Initialize c = W'(xcp - x) = 0. */
02198 i__1 = col2;
02199 for (j = 1; j <= i__1; ++j) {
02200 c__[j] = 0.;
02201 /* L60: */
02202 }
02203 /* Initialize derivative f2. */
02204 f2 = -(*theta) * f1;
02205 if (*col > 0) {
02206 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info);
02207 if (*info != 0) {
02208 return 0;
02209 }
02210 f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
02211 }
02212 dtm = -f1 / f2;
02213 tsum = 0.;
02214 *nint = 1;
02215 if (*iprint >= 99) {
02216 /* s_wsle(&io___109);
02217 do_lio(&c__9, &c__1, "There are ", (long int)10);
02218 do_lio(&c__3, &c__1, (char *)&nbreak, (long int)sizeof(long int));
02219 do_lio(&c__9, &c__1, " breakpoints ", (long int)14);
02220 e_wsle();*/
02221 }
02222 /* If there are no breakpoints, locate the GCP and return. */
02223 if (nbreak == 0) {
02224 goto L888;
02225 }
02226 nleft = nbreak;
02227 iter = 1;
02228 tj = 0.;
02229 /* ------------------- the beginning of the loop ------------------------- */
02230 L777:
02231 /* Find the next smallest breakpoint; */
02232 /* compute dt = t(nleft) - t(nleft + 1). */
02233 tj0 = tj;
02234 if (iter == 1) {
02235 /* Since we already have the smallest breakpoint we need not do */
02236 /* heapsort yet. Often only one breakpoint is used and the */
02237 /* cost of heapsort is avoided. */
02238 tj = bkmin;
02239 ibp = iorder[ibkmin];
02240 } else {
02241 if (iter == 2) {
02242 /* Replace the already used smallest breakpoint with the */
02243 /* breakpoint numbered nbreak > nlast, before heapsort call. */
02244 if (ibkmin != nbreak) {
02245 t[ibkmin] = t[nbreak];
02246 iorder[ibkmin] = iorder[nbreak];
02247 }
02248 /* Update heap structure of breakpoints */
02249 /* (if iter=2, initialize heap). */
02250 }
02251 i__1 = iter - 2;
02252 hpsolb_(&nleft, &t[1], &iorder[1], &i__1);
02253 tj = t[nleft];
02254 ibp = iorder[nleft];
02255 }
02256 dt = tj - tj0;
02257 if (dt != 0. && *iprint >= 100) {
02258 /* s_wsfe(&io___116);
02259 do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
02260 do_fio(&c__1, (char *)&f1, (long int)sizeof(double));
02261 do_fio(&c__1, (char *)&f2, (long int)sizeof(double));
02262 e_wsfe();
02263 s_wsfe(&io___117);
02264 do_fio(&c__1, (char *)&dt, (long int)sizeof(double));
02265 e_wsfe();
02266 s_wsfe(&io___118);
02267 do_fio(&c__1, (char *)&dtm, (long int)sizeof(double));
02268 e_wsfe();*/
02269 }
02270 /* If a minimizer is within this interval, locate the GCP and return. */
02271 if (dtm < dt) {
02272 goto L888;
02273 }
02274 /* Otherwise fix one variable and */
02275 /* reset the corresponding component of d to zero. */
02276 tsum += dt;
02277 --nleft;
02278 ++iter;
02279 dibp = d__[ibp];
02280 d__[ibp] = 0.;
02281 if (dibp > 0.) {
02282 zibp = u[ibp] - x[ibp];
02283 xcp[ibp] = u[ibp];
02284 iwhere[ibp] = 2;
02285 } else {
02286 zibp = l[ibp] - x[ibp];
02287 xcp[ibp] = l[ibp];
02288 iwhere[ibp] = 1;
02289 }
02290 if (*iprint >= 100) {
02291 /* s_wsle(&io___121);
02292 do_lio(&c__9, &c__1, "Variable ", (long int)10);
02293 do_lio(&c__3, &c__1, (char *)&ibp, (long int)sizeof(long int));
02294 do_lio(&c__9, &c__1, " is fixed.", (long int)11);
02295 e_wsle();*/
02296 }
02297 if (nleft == 0 && nbreak == *n) {
02298 /* all n variables are fixed, */
02299 /* return with xcp as GCP. */
02300 dtm = dt;
02301 goto L999;
02302 }
02303 /* Update the derivative information. */
02304 ++(*nint);
02305 /* Computing 2nd power */
02306 d__1 = dibp;
02307 dibp2 = d__1 * d__1;
02308 /* Update f1 and f2. */
02309 /* temporarily set f1 and f2 for col=0. */
02310 f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp;
02311 f2 -= *theta * dibp2;
02312 if (*col > 0) {
02313 /* update c = c + dt*p. */
02314 daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
02315 /* choose wbp, */
02316 /* the row of W corresponding to the breakpoint encountered. */
02317 pointr = *head;
02318 i__1 = *col;
02319 for (j = 1; j <= i__1; ++j) {
02320 wbp[j] = wy[ibp + pointr * wy_dim1];
02321 wbp[*col + j] = *theta * ws[ibp + pointr * ws_dim1];
02322 pointr = pointr % *m + 1;
02323 /* L70: */
02324 }
02325 /* compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */
02326 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info);
02327 if (*info != 0) {
02328 return 0;
02329 }
02330 wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
02331 wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
02332 wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
02333 /* update p = p - dibp*wbp. */
02334 d__1 = -dibp;
02335 daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
02336 /* complete updating f1 and f2 while col > 0. */
02337 f1 += dibp * wmc;
02338 f2 = f2 + dibp * 2. * wmp - dibp2 * wmw;
02339 }
02340 if (nleft > 0) {
02341 dtm = -f1 / f2;
02342 goto L777;
02343 /* to repeat the loop for unsearched intervals. */
02344 } else if (bnded) {
02345 f1 = 0.;
02346 f2 = 0.;
02347 dtm = 0.;
02348 } else {
02349 dtm = -f1 / f2;
02350 }
02351 /* ------------------- the end of the loop ------------------------------- */
02352 L888:
02353 if (*iprint >= 99) {
02354 /* s_wsle(&io___126);
02355 e_wsle();
02356 s_wsle(&io___127);
02357 do_lio(&c__9, &c__1, "GCP found in this segment", (long int)25);
02358 e_wsle();
02359 s_wsfe(&io___128);
02360 do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
02361 do_fio(&c__1, (char *)&f1, (long int)sizeof(double));
02362 do_fio(&c__1, (char *)&f2, (long int)sizeof(double));
02363 e_wsfe();
02364 s_wsfe(&io___129);
02365 do_fio(&c__1, (char *)&dtm, (long int)sizeof(double));
02366 e_wsfe();*/
02367 }
02368 if (dtm <= 0.) {
02369 dtm = 0.;
02370 }
02371 tsum += dtm;
02372 /* Move free variables (i.e., the ones w/o breakpoints) and */
02373 /* the variables whose breakpoints haven't been reached. */
02374 daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
02375 L999:
02376 /* Update c = c + dtm*p = W'(x^c - x) */
02377 /* which will be used in computing r = Z'(B(x^c - x) + g). */
02378 if (*col > 0) {
02379 daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
02380 }
02381 if (*iprint > 100) {
02382 /* s_wsfe(&io___130);
02383 i__1 = *n;
02384 for (i__ = 1; i__ <= i__1; ++i__) {
02385 do_fio(&c__1, (char *)&xcp[i__], (long int)sizeof(double));
02386 }
02387 e_wsfe();*/
02388 }
02389 if (*iprint >= 99) {
02390 /* s_wsfe(&io___131);
02391 e_wsfe();*/
02392 }
02393 return 0;
02394 } /* cauchy_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2397 of file lbfgsb.cpp. References bmv_(), theta, v, and x. Referenced by mainlb_(). 02407 {
02408 /* System generated locals */
02409 long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
02410 wt_dim1, wt_offset, i__1, i__2;
02411
02412 /* Local variables */
02413 static long int i__, j, k;
02414 static double a1, a2;
02415 static long int pointr;
02416 extern /* Subroutine */ int bmv_(long int *m, double *sy, double *wt, long int *col, double *v, double *p, long int *info);
02417
02418 /* ************ */
02419
02420 /* Subroutine cmprlb */
02421
02422 /* This subroutine computes r=-Z'B(xcp-xk)-Z'g by using */
02423 /* wa(2m+1)=W'(xcp-x) from subroutine cauchy. */
02424
02425 /* Subprograms called: */
02426
02427 /* L-BFGS-B Library ... bmv. */
02428
02429
02430 /* * * * */
02431
02432 /* NEOS, November 1994. (Latest revision June 1996.) */
02433 /* Optimization Technology Center. */
02434 /* Argonne National Laboratory and Northwestern University. */
02435 /* Written by */
02436 /* Ciyou Zhu */
02437 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02438
02439
02440 /* ************ */
02441 /* Parameter adjustments */
02442 --index;
02443 --r__;
02444 --z__;
02445 --g;
02446 --x;
02447 --wa;
02448 wt_dim1 = *m;
02449 wt_offset = 1 + wt_dim1 * 1;
02450 wt -= wt_offset;
02451 sy_dim1 = *m;
02452 sy_offset = 1 + sy_dim1 * 1;
02453 sy -= sy_offset;
02454 wy_dim1 = *n;
02455 wy_offset = 1 + wy_dim1 * 1;
02456 wy -= wy_offset;
02457 ws_dim1 = *n;
02458 ws_offset = 1 + ws_dim1 * 1;
02459 ws -= ws_offset;
02460
02461 /* Function Body */
02462 if (! (*cnstnd) && *col > 0) {
02463 i__1 = *n;
02464 for (i__ = 1; i__ <= i__1; ++i__) {
02465 r__[i__] = -g[i__];
02466 /* L26: */
02467 }
02468 } else {
02469 i__1 = *nfree;
02470 for (i__ = 1; i__ <= i__1; ++i__) {
02471 k = index[i__];
02472 r__[i__] = -(*theta) * (z__[k] - x[k]) - g[k];
02473 /* L30: */
02474 }
02475 bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wa[(*m << 1) + 1], &wa[
02476 1], info);
02477 if (*info != 0) {
02478 *info = -8;
02479 return 0;
02480 }
02481 pointr = *head;
02482 i__1 = *col;
02483 for (j = 1; j <= i__1; ++j) {
02484 a1 = wa[j];
02485 a2 = *theta * wa[*col + j];
02486 i__2 = *nfree;
02487 for (i__ = 1; i__ <= i__2; ++i__) {
02488 k = index[i__];
02489 r__[i__] = r__[i__] + wy[k + pointr * wy_dim1] * a1 + ws[k +
02490 pointr * ws_dim1] * a2;
02491 /* L32: */
02492 }
02493 pointr = pointr % *m + 1;
02494 /* L34: */
02495 }
02496 }
02497 return 0;
02498 } /* cmprlb_ */
|
|
||||||||||||||||||||||||||||
|
Definition at line 5409 of file lbfgsb.cpp. Referenced by cauchy_(), and dtrsl_(). 05415 {
05416 /* System generated locals */
05417 long int i__1;
05418
05419 /* Local variables */
05420 static long int i__, m, ix, iy, mp1;
05421
05422
05423 /* constant times a vector plus a vector. */
05424 /* uses unrolled loops for increments equal to one. */
05425 /* jack dongarra, linpack, 3/11/78. */
05426
05427
05428 /* Parameter adjustments */
05429 --dy;
05430 --dx;
05431
05432 /* Function Body */
05433 if (*n <= 0) {
05434 return 0;
05435 }
05436 if (*da == 0.) {
05437 return 0;
05438 }
05439 if (*incx == 1 && *incy == 1) {
05440 goto L20;
05441 }
05442
05443 /* code for unequal increments or equal increments */
05444 /* not equal to 1 */
05445
05446 ix = 1;
05447 iy = 1;
05448 if (*incx < 0) {
05449 ix = (-(*n) + 1) * *incx + 1;
05450 }
05451 if (*incy < 0) {
05452 iy = (-(*n) + 1) * *incy + 1;
05453 }
05454 i__1 = *n;
05455 for (i__ = 1; i__ <= i__1; ++i__) {
05456 dy[iy] += *da * dx[ix];
05457 ix += *incx;
05458 iy += *incy;
05459 /* L10: */
05460 }
05461 return 0;
05462
05463 /* code for both increments equal to 1 */
05464
05465
05466 /* clean-up loop */
05467
05468 L20:
05469 m = *n % 4;
05470 if (m == 0) {
05471 goto L40;
05472 }
05473 i__1 = m;
05474 for (i__ = 1; i__ <= i__1; ++i__) {
05475 dy[i__] += *da * dx[i__];
05476 /* L30: */
05477 }
05478 if (*n < 4) {
05479 return 0;
05480 }
05481 L40:
05482 mp1 = m + 1;
05483 i__1 = *n;
05484 for (i__ = mp1; i__ <= i__1; i__ += 4) {
05485 dy[i__] += *da * dx[i__];
05486 dy[i__ + 1] += *da * dx[i__ + 1];
05487 dy[i__ + 2] += *da * dx[i__ + 2];
05488 dy[i__ + 3] += *da * dx[i__ + 3];
05489 /* L50: */
05490 }
05491 return 0;
05492 } /* daxpy_ */
|
|
||||||||||||||||||||||||
|
Definition at line 5495 of file lbfgsb.cpp. Referenced by cauchy_(), formk_(), lnsrlb_(), mainlb_(), and matupd_(). 05501 {
05502 /* System generated locals */
05503 long int i__1;
05504
05505 /* Local variables */
05506 static long int i__, m, ix, iy, mp1;
05507
05508
05509 /* copies a vector, x, to a vector, y. */
05510 /* uses unrolled loops for increments equal to one. */
05511 /* jack dongarra, linpack, 3/11/78. */
05512
05513
05514 /* Parameter adjustments */
05515 --dy;
05516 --dx;
05517
05518 /* Function Body */
05519 if (*n <= 0) {
05520 return 0;
05521 }
05522 if (*incx == 1 && *incy == 1) {
05523 goto L20;
05524 }
05525
05526 /* code for unequal increments or equal increments */
05527 /* not equal to 1 */
05528
05529 ix = 1;
05530 iy = 1;
05531 if (*incx < 0) {
05532 ix = (-(*n) + 1) * *incx + 1;
05533 }
05534 if (*incy < 0) {
05535 iy = (-(*n) + 1) * *incy + 1;
05536 }
05537 i__1 = *n;
05538 for (i__ = 1; i__ <= i__1; ++i__) {
05539 dy[iy] = dx[ix];
05540 ix += *incx;
05541 iy += *incy;
05542 /* L10: */
05543 }
05544 return 0;
05545
05546 /* code for both increments equal to 1 */
05547
05548
05549 /* clean-up loop */
05550
05551 L20:
05552 m = *n % 7;
05553 if (m == 0) {
05554 goto L40;
05555 }
05556 i__1 = m;
05557 for (i__ = 1; i__ <= i__1; ++i__) {
05558 dy[i__] = dx[i__];
05559 /* L30: */
05560 }
05561 if (*n < 7) {
05562 return 0;
05563 }
05564 L40:
05565 mp1 = m + 1;
05566 i__1 = *n;
05567 for (i__ = mp1; i__ <= i__1; i__ += 7) {
05568 dy[i__] = dx[i__];
05569 dy[i__ + 1] = dx[i__ + 1];
05570 dy[i__ + 2] = dx[i__ + 2];
05571 dy[i__ + 3] = dx[i__ + 3];
05572 dy[i__ + 4] = dx[i__ + 4];
05573 dy[i__ + 5] = dx[i__ + 5];
05574 dy[i__ + 6] = dx[i__ + 6];
05575 /* L50: */
05576 }
05577 return 0;
05578 } /* dcopy_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 4548 of file lbfgsb.cpp. References abs_, dcstep_(), max_, min_, s_cmp(), and s_copy(). Referenced by lnsrlb_(). 04555 {
04556 /* System generated locals */
04557 double d__1;
04558
04559 /* Builtin functions */
04560 long int s_cmp(char *, const char *const, long int, long int);
04561 /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
04562
04563 /* Local variables */
04564 static long int stage;
04565 static double finit, ginit, width, ftest, gtest, stmin, stmax, width1,
04566 fm, gm, fx, fy, gx, gy;
04567 static long int brackt;
04568 //extern /* Subroutine */ int dcstep_();
04569 extern int dcstep_(double *stx, double *fx, double *dx, double *sty, double *fy, double *dy,
04570 double *stp, double *fp, double *dp, long int *brackt, double *stpmin, double *stpmax);
04571
04572 static double fxm, fym, gxm, gym, stx, sty;
04573
04574 /* ********** */
04575
04576 /* Subroutine dcsrch */
04577
04578 /* This subroutine finds a step that satisfies a sufficient */
04579 /* decrease condition and a curvature condition. */
04580
04581 /* Each call of the subroutine updates an interval with */
04582 /* endpoints stx and sty. The interval is initially chosen */
04583 /* so that it contains a minimizer of the modified function */
04584
04585 /* psi(stp) = f(stp) - f(0) - ftol*stp*f'(0). */
04586
04587 /* If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
04588 /* interval is chosen so that it contains a minimizer of f. */
04589
04590 /* The algorithm is designed to find a step that satisfies */
04591 /* the sufficient decrease condition */
04592
04593 /* f(stp) <= f(0) + ftol*stp*f'(0), */
04594
04595 /* and the curvature condition */
04596
04597 /* abs_(f'(stp)) <= gtol*abs_(f'(0)). */
04598
04599 /* If ftol is less than gtol and if, for example, the function */
04600 /* is bounded below, then there is always a step which satisfies */
04601 /* both conditions. */
04602
04603 /* If no step can be found that satisfies both conditions, then */
04604 /* the algorithm stops with a warning. In this case stp only */
04605 /* satisfies the sufficient decrease condition. */
04606
04607 /* A typical invocation of dcsrch has the following outline: */
04608
04609 /* task = 'START' */
04610 /* 10 continue */
04611 /* call dcsrch( ... ) */
04612 /* if (task .eq. 'FG') then */
04613 /* Evaluate the function and the gradient at stp */
04614 /* goto 10 */
04615 /* end if */
04616
04617 /* NOTE: The user must no alter work arrays between calls. */
04618
04619 /* The subroutine statement is */
04620
04621 /* subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, */
04622 /* task,isave,dsave) */
04623 /* where */
04624
04625 /* f is a double precision variable. */
04626 /* On initial entry f is the value of the function at 0. */
04627 /* On subsequent entries f is the value of the */
04628 /* function at stp. */
04629 /* On exit f is the value of the function at stp. */
04630
04631 /* g is a double precision variable. */
04632 /* On initial entry g is the derivative of the function at 0. */
04633 /* On subsequent entries g is the derivative of the */
04634 /* function at stp. */
04635 /* On exit g is the derivative of the function at stp. */
04636
04637 /* stp is a double precision variable. */
04638 /* On entry stp is the current estimate of a satisfactory */
04639 /* step. On initial entry, a positive initial estimate */
04640 /* must be provided. */
04641 /* On exit stp is the current estimate of a satisfactory step */
04642 /* if task = 'FG'. If task = 'CONV' then stp satisfies */
04643 /* the sufficient decrease and curvature condition. */
04644
04645 /* ftol is a double precision variable. */
04646 /* On entry ftol specifies a nonnegative tolerance for the */
04647 /* sufficient decrease condition. */
04648 /* On exit ftol is unchanged. */
04649
04650 /* gtol is a double precision variable. */
04651 /* On entry gtol specifies a nonnegative tolerance for the */
04652 /* curvature condition. */
04653 /* On exit gtol is unchanged. */
04654
04655 /* xtol is a double precision variable. */
04656 /* On entry xtol specifies a nonnegative relative tolerance */
04657 /* for an acceptable step. The subroutine exits with a */
04658 /* warning if the relative difference between sty and stx */
04659 /* is less than xtol. */
04660 /* On exit xtol is unchanged. */
04661
04662 /* stpmin is a double precision variable. */
04663 /* On entry stpmin is a nonnegative lower bound for the step. */
04664 /* On exit stpmin is unchanged. */
04665
04666 /* stpmax is a double precision variable. */
04667 /* On entry stpmax is a nonnegative upper bound for the step. */
04668 /* On exit stpmax is unchanged. */
04669
04670 /* task is a character variable of length at least 60. */
04671 /* On initial entry task must be set to 'START'. */
04672 /* On exit task indicates the required action: */
04673
04674 /* If task(1:2) = 'FG' then evaluate the function and */
04675 /* derivative at stp and call dcsrch again. */
04676
04677 /* If task(1:4) = 'CONV' then the search is successful. */
04678
04679 /* If task(1:4) = 'WARN' then the subroutine is not able */
04680 /* to satisfy the convergence conditions. The exit value of */
04681 /* stp contains the best point found during the search. */
04682
04683 /* If task(1:5) = 'ERROR' then there is an error in the */
04684 /* input arguments. */
04685
04686 /* On exit with convergence, a warning or an error, the */
04687 /* variable task contains additional information. */
04688
04689 /* isave is an long int work array of dimension 2. */
04690
04691 /* dsave is a double precision work array of dimension 13. */
04692
04693 /* Subprograms called */
04694
04695 /* MINPACK-2 ... dcstep */
04696
04697 /* MINPACK-1 Project. June 1983. */
04698 /* Argonne National Laboratory. */
04699 /* Jorge J. More' and David J. Thuente. */
04700
04701 /* MINPACK-2 Project. October 1993. */
04702 /* Argonne National Laboratory and University of Minnesota. */
04703 /* Brett M. Averick, Richard G. Carter, and Jorge J. More'. */
04704
04705 /* ********** */
04706 /* Initialization block. */
04707 /* Parameter adjustments */
04708 --dsave;
04709 --isave;
04710
04711 /* Function Body */
04712 if (s_cmp(task, "START", (long int)5, (long int)5) == 0) {
04713 /* Check the input arguments for errors. */
04714 if (*stp < *stpmin) {
04715 s_copy(task, "ERROR: STP .LT. STPMIN", task_len, (long int)22);
04716 }
04717 if (*stp > *stpmax) {
04718 s_copy(task, "ERROR: STP .GT. STPMAX", task_len, (long int)22);
04719 }
04720 if (*g >= 0.) {
04721 s_copy(task, "ERROR: INITIAL G .GE. ZERO", task_len, (long int)26);
04722 }
04723 if (*ftol < 0.) {
04724 s_copy(task, "ERROR: FTOL .LT. ZERO", task_len, (long int)21);
04725 }
04726 if (*gtol < 0.) {
04727 s_copy(task, "ERROR: GTOL .LT. ZERO", task_len, (long int)21);
04728 }
04729 if (*xtol < 0.) {
04730 s_copy(task, "ERROR: XTOL .LT. ZERO", task_len, (long int)21);
04731 }
04732 if (*stpmin < 0.) {
04733 s_copy(task, "ERROR: STPMIN .LT. ZERO", task_len, (long int)23);
04734 }
04735 if (*stpmax < *stpmin) {
04736 s_copy(task, "ERROR: STPMAX .LT. STPMIN", task_len, (long int)25);
04737 }
04738 /* Exit if there are errors on input. */
04739 if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
04740 return 0;
04741 }
04742 /* Initialize local variables. */
04743 brackt = FALSE_;
04744 stage = 1;
04745 finit = *f;
04746 ginit = *g;
04747 gtest = *ftol * ginit;
04748 width = *stpmax - *stpmin;
04749 width1 = width / .5;
04750 /* The variables stx, fx, gx contain the values of the step, */
04751 /* function, and derivative at the best step. */
04752 /* The variables sty, fy, gy contain the value of the step, */
04753 /* function, and derivative at sty. */
04754 /* The variables stp, f, g contain the values of the step, */
04755 /* function, and derivative at stp. */
04756 stx = 0.;
04757 fx = finit;
04758 gx = ginit;
04759 sty = 0.;
04760 fy = finit;
04761 gy = ginit;
04762 stmin = 0.;
04763 stmax = *stp + *stp * 4.;
04764 s_copy(task, "FG", task_len, (long int)2);
04765 goto L1000;
04766 } else {
04767 /* Restore local variables. */
04768 if (isave[1] == 1) {
04769 brackt = TRUE_;
04770 } else {
04771 brackt = FALSE_;
04772 }
04773 stage = isave[2];
04774 ginit = dsave[1];
04775 gtest = dsave[2];
04776 gx = dsave[3];
04777 gy = dsave[4];
04778 finit = dsave[5];
04779 fx = dsave[6];
04780 fy = dsave[7];
04781 stx = dsave[8];
04782 sty = dsave[9];
04783 stmin = dsave[10];
04784 stmax = dsave[11];
04785 width = dsave[12];
04786 width1 = dsave[13];
04787 }
04788 /* If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the */
04789 /* algorithm enters the second stage. */
04790 ftest = finit + *stp * gtest;
04791 if (stage == 1 && *f <= ftest && *g >= 0.) {
04792 stage = 2;
04793 }
04794 /* Test for warnings. */
04795 if (brackt && (*stp <= stmin || *stp >= stmax)) {
04796 s_copy(task, "WARNING: ROUNDING ERRORS PREVENT PROGRESS", task_len, (
04797 long int)41);
04798 }
04799 if (brackt && stmax - stmin <= *xtol * stmax) {
04800 s_copy(task, "WARNING: XTOL TEST SATISFIED", task_len, (long int)28);
04801 }
04802 if (*stp == *stpmax && *f <= ftest && *g <= gtest) {
04803 s_copy(task, "WARNING: STP = STPMAX", task_len, (long int)21);
04804 }
04805 if (*stp == *stpmin && (*f > ftest || *g >= gtest)) {
04806 s_copy(task, "WARNING: STP = STPMIN", task_len, (long int)21);
04807 }
04808 /* Test for convergence. */
04809 if (*f <= ftest && abs_(*g) <= *gtol * (-ginit)) {
04810 s_copy(task, "CONVERGENCE", task_len, (long int)11);
04811 }
04812 /* Test for termination. */
04813 if (s_cmp(task, "WARN", (long int)4, (long int)4) == 0 || s_cmp(task, "CONV",
04814 (long int)4, (long int)4) == 0) {
04815 goto L1000;
04816 }
04817 /* A modified function is used to predict the step during the */
04818 /* first stage if a lower function value has been obtained but */
04819 /* the decrease is not sufficient. */
04820 if (stage == 1 && *f <= fx && *f > ftest) {
04821 /* Define the modified function and derivative values. */
04822 fm = *f - *stp * gtest;
04823 fxm = fx - stx * gtest;
04824 fym = fy - sty * gtest;
04825 gm = *g - gtest;
04826 gxm = gx - gtest;
04827 gym = gy - gtest;
04828 /* Call dcstep to update stx, sty, and to compute the new step. */
04829 dcstep_(&stx, &fxm, &gxm, &sty, &fym, &gym, stp, &fm, &gm, &brackt, &
04830 stmin, &stmax);
04831 /* Reset the function and derivative values for f. */
04832 fx = fxm + stx * gtest;
04833 fy = fym + sty * gtest;
04834 gx = gxm + gtest;
04835 gy = gym + gtest;
04836 } else {
04837 /* Call dcstep to update stx, sty, and to compute the new step. */
04838 dcstep_(&stx, &fx, &gx, &sty, &fy, &gy, stp, f, g, &brackt, &stmin, &
04839 stmax);
04840 }
04841 /* Decide if a bisection step is needed. */
04842 if (brackt) {
04843 if ((d__1 = sty - stx, abs_(d__1)) >= width1 * .66) {
04844 *stp = stx + (sty - stx) * .5;
04845 }
04846 width1 = width;
04847 width = (d__1 = sty - stx, abs_(d__1));
04848 }
04849 /* Set the minimum and maximum steps allowed for stp. */
04850 if (brackt) {
04851 stmin = min_(stx,sty);
04852 stmax = max_(stx,sty);
04853 } else {
04854 stmin = *stp + (*stp - stx) * 1.1;
04855 stmax = *stp + (*stp - stx) * 4.;
04856 }
04857 /* Force the step to be within the bounds stpmax and stpmin. */
04858 *stp = max_(*stp,*stpmin);
04859 *stp = min_(*stp,*stpmax);
04860 /* If further progress is not possible, let stp be the best */
04861 /* point obtained during the search. */
04862 if (brackt && (*stp <= stmin || *stp >= stmax) || brackt && stmax - stmin
04863 <= *xtol * stmax) {
04864 *stp = stx;
04865 }
04866 /* Obtain another function and derivative. */
04867 s_copy(task, "FG", task_len, (long int)2);
04868 L1000:
04869 /* Save local variables. */
04870 if (brackt) {
04871 isave[1] = 1;
04872 } else {
04873 isave[1] = 0;
04874 }
04875 isave[2] = stage;
04876 dsave[1] = ginit;
04877 dsave[2] = gtest;
04878 dsave[3] = gx;
04879 dsave[4] = gy;
04880 dsave[5] = finit;
04881 dsave[6] = fx;
04882 dsave[7] = fy;
04883 dsave[8] = stx;
04884 dsave[9] = sty;
04885 dsave[10] = stmin;
04886 dsave[11] = stmax;
04887 dsave[12] = width;
04888 dsave[13] = width1;
04889 return 0;
04890 } /* dcsrch_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 4893 of file lbfgsb.cpp. References abs_, max_, min_, q, sqrt(), and theta. Referenced by dcsrch_(). 04898 {
04899 /* System generated locals */
04900 double d__1, d__2, d__3;
04901
04902 /* Builtin functions */
04903 //double sqrt();
04904
04905 /* Local variables */
04906 static double sgnd, stpc, stpf, stpq, p, q, gamma, r__, s, theta;
04907
04908 /* ********** */
04909
04910 /* Subroutine dcstep */
04911
04912 /* This subroutine computes a safeguarded step for a search */
04913 /* procedure and updates an interval that contains a step that */
04914 /* satisfies a sufficient decrease and a curvature condition. */
04915
04916 /* The parameter stx contains the step with the least function */
04917 /* value. If brackt is set to .true. then a minimizer has */
04918 /* been bracketed in an interval with endpoints stx and sty. */
04919 /* The parameter stp contains the current step. */
04920 /* The subroutine assumes that if brackt is set to .true. then */
04921
04922 /* min_(stx,sty) < stp < max_(stx,sty), */
04923
04924 /* and that the derivative at stx is negative in the direction */
04925 /* of the step. */
04926
04927 /* The subroutine statement is */
04928
04929 /* subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, */
04930 /* stpmin,stpmax) */
04931
04932 /* where */
04933
04934 /* stx is a double precision variable. */
04935 /* On entry stx is the best step obtained so far and is an */
04936 /* endpoint of the interval that contains the minimizer. */
04937 /* On exit stx is the updated best step. */
04938
04939 /* fx is a double precision variable. */
04940 /* On entry fx is the function at stx. */
04941 /* On exit fx is the function at stx. */
04942
04943 /* dx is a double precision variable. */
04944 /* On entry dx is the derivative of the function at */
04945 /* stx. The derivative must be negative in the direction of */
04946 /* the step, that is, dx and stp - stx must have opposite */
04947 /* signs. */
04948 /* On exit dx is the derivative of the function at stx. */
04949
04950 /* sty is a double precision variable. */
04951 /* On entry sty is the second endpoint of the interval that */
04952 /* contains the minimizer. */
04953 /* On exit sty is the updated endpoint of the interval that */
04954 /* contains the minimizer. */
04955
04956 /* fy is a double precision variable. */
04957 /* On entry fy is the function at sty. */
04958 /* On exit fy is the function at sty. */
04959
04960 /* dy is a double precision variable. */
04961 /* On entry dy is the derivative of the function at sty. */
04962 /* On exit dy is the derivative of the function at the exit sty. */
04963
04964 /* stp is a double precision variable. */
04965 /* On entry stp is the current step. If brackt is set to .true. */
04966 /* then on input stp must be between stx and sty. */
04967 /* On exit stp is a new trial step. */
04968
04969 /* fp is a double precision variable. */
04970 /* On entry fp is the function at stp */
04971 /* On exit fp is unchanged. */
04972
04973 /* dp is a double precision variable. */
04974 /* On entry dp is the the derivative of the function at stp. */
04975 /* On exit dp is unchanged. */
04976
04977 /* brackt is an long int variable. */
04978 /* On entry brackt specifies if a minimizer has been bracketed. */
04979 /* Initially brackt must be set to .false. */
04980 /* On exit brackt specifies if a minimizer has been bracketed. */
04981 /* When a minimizer is bracketed brackt is set to .true. */
04982
04983 /* stpmin is a double precision variable. */
04984 /* On entry stpmin is a lower bound for the step. */
04985 /* On exit stpmin is unchanged. */
04986
04987 /* stpmax is a double precision variable. */
04988 /* On entry stpmax is an upper bound for the step. */
04989 /* On exit stpmax is unchanged. */
04990
04991 /* MINPACK-1 Project. June 1983 */
04992 /* Argonne National Laboratory. */
04993 /* Jorge J. More' and David J. Thuente. */
04994
04995 /* MINPACK-2 Project. October 1993. */
04996 /* Argonne National Laboratory and University of Minnesota. */
04997 /* Brett M. Averick and Jorge J. More'. */
04998
04999 /* ********** */
05000 sgnd = *dp * (*dx / abs_(*dx));
05001 /* First case: A higher function value. The minimum is bracketed. */
05002 /* If the cubic step is closer to stx than the quadratic step, the */
05003 /* cubic step is taken, otherwise the average of the cubic and */
05004 /* quadratic steps is taken. */
05005 if (*fp > *fx) {
05006 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05007 /* Computing MAX */
05008 d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05009 *dp);
05010 s = max_(d__1,d__2);
05011 /* Computing 2nd power */
05012 d__1 = theta / s;
05013 gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
05014 if (*stp < *stx) {
05015 gamma = -gamma;
05016 }
05017 p = gamma - *dx + theta;
05018 q = gamma - *dx + gamma + *dp;
05019 r__ = p / q;
05020 stpc = *stx + r__ * (*stp - *stx);
05021 stpq = *stx + *dx / ((*fx - *fp) / (*stp - *stx) + *dx) / 2. * (*stp
05022 - *stx);
05023 if ((d__1 = stpc - *stx, abs_(d__1)) < (d__2 = stpq - *stx, abs_(d__2)))
05024 {
05025 stpf = stpc;
05026 } else {
05027 stpf = stpc + (stpq - stpc) / 2.;
05028 }
05029 *brackt = TRUE_;
05030 /* Second case: A lower function value and derivatives of opposite */
05031 /* sign. The minimum is bracketed. If the cubic step is farther from */
05032 /* stp than the secant step, the cubic step is taken, otherwise the */
05033 /* secant step is taken. */
05034 } else if (sgnd < 0.) {
05035 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05036 /* Computing MAX */
05037 d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05038 *dp);
05039 s = max_(d__1,d__2);
05040 /* Computing 2nd power */
05041 d__1 = theta / s;
05042 gamma = s * sqrt(d__1 * d__1 - *dx / s * (*dp / s));
05043 if (*stp > *stx) {
05044 gamma = -gamma;
05045 }
05046 p = gamma - *dp + theta;
05047 q = gamma - *dp + gamma + *dx;
05048 r__ = p / q;
05049 stpc = *stp + r__ * (*stx - *stp);
05050 stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
05051 if ((d__1 = stpc - *stp, abs_(d__1)) > (d__2 = stpq - *stp, abs_(d__2)))
05052 {
05053 stpf = stpc;
05054 } else {
05055 stpf = stpq;
05056 }
05057 *brackt = TRUE_;
05058 /* Third case: A lower function value, derivatives of the same sign, */
05059 /* and the magnitude of the derivative decreases. */
05060 } else if (abs_(*dp) < abs_(*dx)) {
05061 /* The cubic step is computed only if the cubic tends to infinity */
05062 /* in the direction of the step or if the minimum of the cubic */
05063 /* is beyond stp. Otherwise the cubic step is defined to be the */
05064 /* secant step. */
05065 theta = (*fx - *fp) * 3. / (*stp - *stx) + *dx + *dp;
05066 /* Computing MAX */
05067 d__1 = abs_(theta), d__2 = abs_(*dx), d__1 = max_(d__1,d__2), d__2 = abs_(
05068 *dp);
05069 s = max_(d__1,d__2);
05070 /* The case gamma = 0 only arises if the cubic does not tend */
05071 /* to infinity in the direction of the step. */
05072 /* Computing MAX */
05073 /* Computing 2nd power */
05074 d__3 = theta / s;
05075 d__1 = 0., d__2 = d__3 * d__3 - *dx / s * (*dp / s);
05076 gamma = s * sqrt((max_(d__1,d__2)));
05077 if (*stp > *stx) {
05078 gamma = -gamma;
05079 }
05080 p = gamma - *dp + theta;
05081 q = gamma + (*dx - *dp) + gamma;
05082 r__ = p / q;
05083 if (r__ < 0. && gamma != 0.) {
05084 stpc = *stp + r__ * (*stx - *stp);
05085 } else if (*stp > *stx) {
05086 stpc = *stpmax;
05087 } else {
05088 stpc = *stpmin;
05089 }
05090 stpq = *stp + *dp / (*dp - *dx) * (*stx - *stp);
05091 if (*brackt) {
05092 /* A minimizer has been bracketed. If the cubic step is */
05093 /* closer to stp than the secant step, the cubic step is */
05094 /* taken, otherwise the secant step is taken. */
05095 if ((d__1 = stpc - *stp, abs_(d__1)) < (d__2 = stpq - *stp, abs_(
05096 d__2))) {
05097 stpf = stpc;
05098 } else {
05099 stpf = stpq;
05100 }
05101 if (*stp > *stx) {
05102 /* Computing MIN */
05103 d__1 = *stp + (*sty - *stp) * .66;
05104 stpf = min_(d__1,stpf);
05105 } else {
05106 /* Computing MAX */
05107 d__1 = *stp + (*sty - *stp) * .66;
05108 stpf = max_(d__1,stpf);
05109 }
05110 } else {
05111 /* A minimizer has not been bracketed. If the cubic step is */
05112 /* farther from stp than the secant step, the cubic step is */
05113 /* taken, otherwise the secant step is taken. */
05114 if ((d__1 = stpc - *stp, abs_(d__1)) > (d__2 = stpq - *stp, abs_(
05115 d__2))) {
05116 stpf = stpc;
05117 } else {
05118 stpf = stpq;
05119 }
05120 stpf = min_(*stpmax,stpf);
05121 stpf = max_(*stpmin,stpf);
05122 }
05123 /* Fourth case: A lower function value, derivatives of the same sign, */
05124 /* and the magnitude of the derivative does not decrease. If the */
05125 /* minimum is not bracketed, the step is either stpmin or stpmax, */
05126 /* otherwise the cubic step is taken. */
05127 } else {
05128 if (*brackt) {
05129 theta = (*fp - *fy) * 3. / (*sty - *stp) + *dy + *dp;
05130 /* Computing MAX */
05131 d__1 = abs_(theta), d__2 = abs_(*dy), d__1 = max_(d__1,d__2), d__2 =
05132 abs_(*dp);
05133 s = max_(d__1,d__2);
05134 /* Computing 2nd power */
05135 d__1 = theta / s;
05136 gamma = s * sqrt(d__1 * d__1 - *dy / s * (*dp / s));
05137 if (*stp > *sty) {
05138 gamma = -gamma;
05139 }
05140 p = gamma - *dp + theta;
05141 q = gamma - *dp + gamma + *dy;
05142 r__ = p / q;
05143 stpc = *stp + r__ * (*sty - *stp);
05144 stpf = stpc;
05145 } else if (*stp > *stx) {
05146 stpf = *stpmax;
05147 } else {
05148 stpf = *stpmin;
05149 }
05150 }
05151 /* Update the interval which contains a minimizer. */
05152 if (*fp > *fx) {
05153 *sty = *stp;
05154 *fy = *fp;
05155 *dy = *dp;
05156 } else {
05157 if (sgnd < 0.) {
05158 *sty = *stx;
05159 *fy = *fx;
05160 *dy = *dx;
05161 }
05162 *stx = *stp;
05163 *fx = *fp;
05164 *dx = *dp;
05165 }
05166 /* Compute the new step. */
05167 *stp = stpf;
05168 return 0;
05169 } /* dcstep_ */
|
|
||||||||||||||||||||||||
|
Definition at line 5581 of file lbfgsb.cpp. Referenced by cauchy_(), dpofa_(), dtrsl_(), formk_(), lnsrlb_(), mainlb_(), and matupd_(). 05587 {
05588 /* System generated locals */
05589 long int i__1;
05590 double ret_val;
05591
05592 /* Local variables */
05593 static long int i__, m;
05594 static double dtemp;
05595 static long int ix, iy, mp1;
05596
05597
05598 /* forms the dot product of two vectors. */
05599 /* uses unrolled loops for increments equal to one. */
05600 /* jack dongarra, linpack, 3/11/78. */
05601
05602
05603 /* Parameter adjustments */
05604 --dy;
05605 --dx;
05606
05607 /* Function Body */
05608 ret_val = 0.;
05609 dtemp = 0.;
05610 if (*n <= 0) {
05611 return ret_val;
05612 }
05613 if (*incx == 1 && *incy == 1) {
05614 goto L20;
05615 }
05616
05617 /* code for unequal increments or equal increments */
05618 /* not equal to 1 */
05619
05620 ix = 1;
05621 iy = 1;
05622 if (*incx < 0) {
05623 ix = (-(*n) + 1) * *incx + 1;
05624 }
05625 if (*incy < 0) {
05626 iy = (-(*n) + 1) * *incy + 1;
05627 }
05628 i__1 = *n;
05629 for (i__ = 1; i__ <= i__1; ++i__) {
05630 dtemp += dx[ix] * dy[iy];
05631 ix += *incx;
05632 iy += *incy;
05633 /* L10: */
05634 }
05635 ret_val = dtemp;
05636 return ret_val;
05637
05638 /* code for both increments equal to 1 */
05639
05640
05641 /* clean-up loop */
05642
05643 L20:
05644 m = *n % 5;
05645 if (m == 0) {
05646 goto L40;
05647 }
05648 i__1 = m;
05649 for (i__ = 1; i__ <= i__1; ++i__) {
05650 dtemp += dx[i__] * dy[i__];
05651 /* L30: */
05652 }
05653 if (*n < 5) {
05654 goto L60;
05655 }
05656 L40:
05657 mp1 = m + 1;
05658 i__1 = *n;
05659 for (i__ = mp1; i__ <= i__1; i__ += 5) {
05660 dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
05661 i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ +
05662 4] * dy[i__ + 4];
05663 /* L50: */
05664 }
05665 L60:
05666 ret_val = dtemp;
05667 return ret_val;
05668 } /* ddot_ */
|
|
||||||||||||||||
|
Definition at line 5210 of file lbfgsb.cpp. References abs_, max_, sqrt(), and x. 05214 {
05215 /* System generated locals */
05216 long int i__1, i__2;
05217 double ret_val, d__1, d__2, d__3;
05218
05219 /* Builtin functions */
05220 //double sqrt();
05221
05222 /* Local variables */
05223 static long int i__;
05224 static double scale;
05225
05226 /* ********** */
05227
05228 /* Function dnrm2 */
05229
05230 /* Given a vector x of length n, this function calculates the */
05231 /* Euclidean norm of x with stride incx. */
05232
05233 /* The function statement is */
05234
05235 /* double precision function dnrm2(n,x,incx) */
05236
05237 /* where */
05238
05239 /* n is a positive long int input variable. */
05240
05241 /* x is an input array of length n. */
05242
05243 /* incx is a positive long int variable that specifies the */
05244 /* stride of the vector. */
05245
05246 /* Subprograms called */
05247
05248 /* FORTRAN-supplied ... abs, max, sqrt */
05249
05250 /* MINPACK-2 Project. February 1991. */
05251 /* Argonne National Laboratory. */
05252 /* Brett M. Averick. */
05253
05254 /* ********** */
05255 /* Parameter adjustments */
05256 --x;
05257
05258 /* Function Body */
05259 ret_val = 0.;
05260 scale = 0.;
05261 i__1 = *n;
05262 i__2 = *incx;
05263 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
05264 /* Computing MAX */
05265 d__2 = scale, d__3 = (d__1 = x[i__], abs_(d__1));
05266 scale = max_(d__2,d__3);
05267 /* L10: */
05268 }
05269 if (scale == 0.) {
05270 return ret_val;
05271 }
05272 i__2 = *n;
05273 i__1 = *incx;
05274 for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
05275 /* Computing 2nd power */
05276 d__1 = x[i__] / scale;
05277 ret_val += d__1 * d__1;
05278 /* L20: */
05279 }
05280 ret_val = scale * sqrt(ret_val);
05281 return ret_val;
05282 } /* dnrm2_ */
|
|
|
Definition at line 5285 of file lbfgsb.cpp. References b. Referenced by mainlb_(). 05286 {
05287 /* Initialized data */
05288
05289 static double zero = 0.;
05290 static double one = 1.;
05291 static double two = 2.;
05292
05293 /* System generated locals */
05294 long int i__1;
05295 double ret_val;
05296
05297 /* Local variables */
05298 static double beta;
05299 static long int irnd;
05300 static volatile double temp, temp1, a, b;
05301 static long int i__;
05302 static double betah;
05303 static long int ibeta, negep;
05304 static double tempa;
05305 static long int itemp, it;
05306 static double betain;
05307
05308 /* ********** */
05309
05310 /* Subroutine dpeps */
05311
05312 /* This subroutine computes the machine precision parameter */
05313 /* dpmeps as the smallest floating point number such that */
05314 /* 1 + dpmeps differs from 1. */
05315
05316 /* This subroutine is based on the subroutine machar described in */
05317
05318 /* W. J. Cody, */
05319 /* MACHAR: A subroutine to dynamically determine machine parameters, */
05320 /* ACM Transactions on Mathematical Software, 14, 1988, pages 303-311. */
05321
05322 /* The subroutine statement is: */
05323
05324 /* subroutine dpeps(dpmeps) */
05325
05326 /* where */
05327
05328 /* dpmeps is a double precision variable. */
05329 /* On entry dpmeps need not be specified. */
05330 /* On exit dpmeps is the machine precision. */
05331
05332 /* MINPACK-2 Project. February 1991. */
05333 /* Argonne National Laboratory and University of Minnesota. */
05334 /* Brett M. Averick. */
05335
05336 /* ******* */
05337 /* determine ibeta, beta ala malcolm. */
05338 a = one;
05339 b = one;
05340 L10:
05341 a += a;
05342 temp = a + one;
05343 temp1 = temp - a;
05344 if (temp1 - one == zero) {
05345 goto L10;
05346 }
05347 L20:
05348 b += b;
05349 temp = a + b;
05350 itemp = (long int) (temp - a);
05351 if (itemp == 0) {
05352 goto L20;
05353 }
05354 ibeta = itemp;
05355 beta = (double) ibeta;
05356 /* determine it, irnd. */
05357 it = 0;
05358 b = one;
05359 L30:
05360 ++it;
05361 b *= beta;
05362 temp = b + one;
05363 temp1 = temp - b;
05364 if (temp1 - one == zero) {
05365 goto L30;
05366 }
05367 irnd = 0;
05368 betah = beta / two;
05369 temp = a + betah;
05370 if (temp - a != zero) {
05371 irnd = 1;
05372 }
05373 tempa = a + beta;
05374 temp = tempa + betah;
05375 if (irnd == 0 && temp - tempa != zero) {
05376 irnd = 2;
05377 }
05378 /* determine dpmeps. */
05379 negep = it + 3;
05380 betain = one / beta;
05381 a = one;
05382 i__1 = negep;
05383 for (i__ = 1; i__ <= i__1; ++i__) {
05384 a *= betain;
05385 /* L40: */
05386 }
05387 L50:
05388 temp = one + a;
05389 if (temp - one != zero) {
05390 goto L60;
05391 }
05392 a *= beta;
05393 goto L50;
05394 L60:
05395 ret_val = a;
05396 if (ibeta == 2 || irnd == 0) {
05397 goto L70;
05398 }
05399 a = a * (one + a) / two;
05400 temp = one + a;
05401 if (temp - one != zero) {
05402 ret_val = a;
05403 }
05404 L70:
05405 return ret_val;
05406 } /* dpmeps_ */
|
|
||||||||||||||||||||
|
Definition at line 5671 of file lbfgsb.cpp. References c__1, ddot_(), sqrt(), and t. Referenced by formk_(), and formt_(). 05674 {
05675 /* System generated locals */
05676 long int a_dim1, a_offset, i__1, i__2, i__3;
05677
05678 /* Builtin functions */
05679 //double sqrt();
05680
05681 /* Local variables */
05682 //extern double ddot_();
05683 static long int j, k;
05684 static double s, t;
05685 static long int jm1;
05686
05687
05688 /* dpofa factors a double precision symmetric positive definite */
05689 /* matrix. */
05690
05691 /* dpofa is usually called by dpoco, but it can be called */
05692 /* directly with a saving in time if rcond is not needed. */
05693 /* (time for dpoco) = (1 + 18/n)*(time for dpofa) . */
05694
05695 /* on entry */
05696
05697 /* a double precision(lda, n) */
05698 /* the symmetric matrix to be factored. only the */
05699 /* diagonal and upper triangle are used. */
05700
05701 /* lda long int */
05702 /* the leading dimension of the array a . */
05703
05704 /* n long int */
05705 /* the order of the matrix a . */
05706
05707 /* on return */
05708
05709 /* a an upper triangular matrix r so that a = trans(r)*r */
05710 /* where trans(r) is the transpose. */
05711 /* the strict lower triangle is unaltered. */
05712 /* if info .ne. 0 , the factorization is not complete. */
05713
05714 /* info long int */
05715 /* = 0 for normal return. */
05716 /* = k signals an error condition. the leading minor */
05717 /* of order k is not positive definite. */
05718
05719 /* linpack. this version dated 08/14/78 . */
05720 /* cleve moler, university of new mexico, argonne national lab. */
05721
05722 /* subroutines and functions */
05723
05724 /* blas ddot */
05725 /* fortran sqrt */
05726
05727 /* internal variables */
05728
05729 /* begin block with ...exits to 40 */
05730
05731
05732 /* Parameter adjustments */
05733 a_dim1 = *lda;
05734 a_offset = 1 + a_dim1 * 1;
05735 a -= a_offset;
05736
05737 /* Function Body */
05738 i__1 = *n;
05739 for (j = 1; j <= i__1; ++j) {
05740 *info = j;
05741 s = 0.;
05742 jm1 = j - 1;
05743 if (jm1 < 1) {
05744 goto L20;
05745 }
05746 i__2 = jm1;
05747 for (k = 1; k <= i__2; ++k) {
05748 i__3 = k - 1;
05749 t = a[k + j * a_dim1] - ddot_(&i__3, &a[k * a_dim1 + 1], &c__1, &
05750 a[j * a_dim1 + 1], &c__1);
05751 t /= a[k + k * a_dim1];
05752 a[k + j * a_dim1] = t;
05753 s += t * t;
05754 /* L10: */
05755 }
05756 L20:
05757 s = a[j + j * a_dim1] - s;
05758 /* ......exit */
05759 if (s <= 0.) {
05760 goto L40;
05761 }
05762 a[j + j * a_dim1] = sqrt(s);
05763 /* L30: */
05764 }
05765 *info = 0;
05766 L40:
05767 return 0;
05768 } /* dpofa_ */
|
|
||||||||||||||||||||
|
Definition at line 5771 of file lbfgsb.cpp. Referenced by cauchy_(), and mainlb_(). 05775 {
05776 /* System generated locals */
05777 long int i__1, i__2;
05778
05779 /* Local variables */
05780 static long int i__, m, nincx, mp1;
05781
05782
05783 /* scales a vector by a constant. */
05784 /* uses unrolled loops for increment equal to one. */
05785 /* jack dongarra, linpack, 3/11/78. */
05786 /* modified 3/93 to return if incx .le. 0. */
05787
05788
05789 /* Parameter adjustments */
05790 --dx;
05791
05792 /* Function Body */
05793 if (*n <= 0 || *incx <= 0) {
05794 return 0;
05795 }
05796 if (*incx == 1) {
05797 goto L20;
05798 }
05799
05800 /* code for increment not equal to 1 */
05801
05802 nincx = *n * *incx;
05803 i__1 = nincx;
05804 i__2 = *incx;
05805 for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
05806 dx[i__] = *da * dx[i__];
05807 /* L10: */
05808 }
05809 return 0;
05810
05811 /* code for increment equal to 1 */
05812
05813
05814 /* clean-up loop */
05815
05816 L20:
05817 m = *n % 5;
05818 if (m == 0) {
05819 goto L40;
05820 }
05821 i__2 = m;
05822 for (i__ = 1; i__ <= i__2; ++i__) {
05823 dx[i__] = *da * dx[i__];
05824 /* L30: */
05825 }
05826 if (*n < 5) {
05827 return 0;
05828 }
05829 L40:
05830 mp1 = m + 1;
05831 i__2 = *n;
05832 for (i__ = mp1; i__ <= i__2; i__ += 5) {
05833 dx[i__] = *da * dx[i__];
05834 dx[i__ + 1] = *da * dx[i__ + 1];
05835 dx[i__ + 2] = *da * dx[i__ + 2];
05836 dx[i__ + 3] = *da * dx[i__ + 3];
05837 dx[i__ + 4] = *da * dx[i__ + 4];
05838 /* L50: */
05839 }
05840 return 0;
05841 } /* dscal_ */
|
|
||||||||||||||||||||||||||||
|
Definition at line 5844 of file lbfgsb.cpp. References b, c__1, daxpy_(), ddot_(), and t. Referenced by bmv_(), formk_(), and subsm_(). 05849 {
05850 /* System generated locals */
05851 long int t_dim1, t_offset, i__1, i__2;
05852
05853 /* Local variables */
05854 static long int case__;
05855 //extern double ddot_(void);
05856 static double temp;
05857 static long int j;
05858 //extern /* Subroutine */ int daxpy_(double *a, long int *n, double *c, long int *d, long int *e);
05859 extern int daxpy_(long int *n, double *da, double *dx, long int *incx, double *dy, long int *incy);
05860 static long int jj;
05861
05862
05863
05864 /* dtrsl solves systems of the form */
05865
05866 /* t * x = b */
05867 /* or */
05868 /* trans(t) * x = b */
05869
05870 /* where t is a triangular matrix of order n. here trans(t) */
05871 /* denotes the transpose of the matrix t. */
05872
05873 /* on entry */
05874
05875 /* t double precision(ldt,n) */
05876 /* t contains the matrix of the system. the zero */
05877 /* elements of the matrix are not referenced, and */
05878 /* the corresponding elements of the array can be */
05879 /* used to store other information. */
05880
05881 /* ldt long int */
05882 /* ldt is the leading dimension of the array t. */
05883
05884 /* n long int */
05885 /* n is the order of the system. */
05886
05887 /* b double precision(n). */
05888 /* b contains the right hand side of the system. */
05889
05890 /* job long int */
05891 /* job specifies what kind of system is to be solved. */
05892 /* if job is */
05893
05894 /* 00 solve t*x=b, t lower triangular, */
05895 /* 01 solve t*x=b, t upper triangular, */
05896 /* 10 solve trans(t)*x=b, t lower triangular, */
05897 /* 11 solve trans(t)*x=b, t upper triangular. */
05898
05899 /* on return */
05900
05901 /* b b contains the solution, if info .eq. 0. */
05902 /* otherwise b is unaltered. */
05903
05904 /* info long int */
05905 /* info contains zero if the system is nonsingular. */
05906 /* otherwise info contains the index of */
05907 /* the first zero diagonal element of t. */
05908
05909 /* linpack. this version dated 08/14/78 . */
05910 /* g. w. stewart, university of maryland, argonne national lab. */
05911
05912 /* subroutines and functions */
05913
05914 /* blas daxpy,ddot */
05915 /* fortran mod */
05916
05917 /* internal variables */
05918
05919
05920 /* begin block permitting ...exits to 150 */
05921
05922 /* check for zero diagonal elements. */
05923
05924 /* Parameter adjustments */
05925 t_dim1 = *ldt;
05926 t_offset = 1 + t_dim1 * 1;
05927 t -= t_offset;
05928 --b;
05929
05930 /* Function Body */
05931 i__1 = *n;
05932 for (*info = 1; *info <= i__1; ++(*info)) {
05933 /* ......exit */
05934 if (t[*info + *info * t_dim1] == 0.) {
05935 goto L150;
05936 }
05937 /* L10: */
05938 }
05939 *info = 0;
05940
05941 /* determine the task and go to it. */
05942
05943 case__ = 1;
05944 if (*job % 10 != 0) {
05945 case__ = 2;
05946 }
05947 if (*job % 100 / 10 != 0) {
05948 case__ += 2;
05949 }
05950 switch ((int)case__) {
05951 case 1: goto L20;
05952 case 2: goto L50;
05953 case 3: goto L80;
05954 case 4: goto L110;
05955 }
05956
05957 /* solve t*x=b for t lower triangular */
05958
05959 L20:
05960 b[1] /= t[t_dim1 + 1];
05961 if (*n < 2) {
05962 goto L40;
05963 }
05964 i__1 = *n;
05965 for (j = 2; j <= i__1; ++j) {
05966 temp = -b[j - 1];
05967 i__2 = *n - j + 1;
05968 daxpy_(&i__2, &temp, &t[j + (j - 1) * t_dim1], &c__1, &b[j], &c__1);
05969 b[j] /= t[j + j * t_dim1];
05970 /* L30: */
05971 }
05972 L40:
05973 goto L140;
05974
05975 /* solve t*x=b for t upper triangular. */
05976
05977 L50:
05978 b[*n] /= t[*n + *n * t_dim1];
05979 if (*n < 2) {
05980 goto L70;
05981 }
05982 i__1 = *n;
05983 for (jj = 2; jj <= i__1; ++jj) {
05984 j = *n - jj + 1;
05985 temp = -b[j + 1];
05986 daxpy_(&j, &temp, &t[(j + 1) * t_dim1 + 1], &c__1, &b[1], &c__1);
05987 b[j] /= t[j + j * t_dim1];
05988 /* L60: */
05989 }
05990 L70:
05991 goto L140;
05992
05993 /* solve trans(t)*x=b for t lower triangular. */
05994
05995 L80:
05996 b[*n] /= t[*n + *n * t_dim1];
05997 if (*n < 2) {
05998 goto L100;
05999 }
06000 i__1 = *n;
06001 for (jj = 2; jj <= i__1; ++jj) {
06002 j = *n - jj + 1;
06003 i__2 = jj - 1;
06004 b[j] -= ddot_(&i__2, &t[j + 1 + j * t_dim1], &c__1, &b[j + 1], &c__1);
06005 b[j] /= t[j + j * t_dim1];
06006 /* L90: */
06007 }
06008 L100:
06009 goto L140;
06010
06011 /* solve trans(t)*x=b for t upper triangular. */
06012
06013 L110:
06014 b[1] /= t[t_dim1 + 1];
06015 if (*n < 2) {
06016 goto L130;
06017 }
06018 i__1 = *n;
06019 for (j = 2; j <= i__1; ++j) {
06020 i__2 = j - 1;
06021 b[j] -= ddot_(&i__2, &t[j * t_dim1 + 1], &c__1, &b[1], &c__1);
06022 b[j] /= t[j + j * t_dim1];
06023 /* L120: */
06024 }
06025 L130:
06026 L140:
06027 L150:
06028 return 0;
06029 } /* dtrsl_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2501 of file lbfgsb.cpp. References s_copy(). Referenced by mainlb_(). 02509 {
02510 /* System generated locals */
02511 long int i__1;
02512
02513 /* Builtin functions */
02514 /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
02515
02516 /* Local variables */
02517 static long int i__;
02518
02519 /* ************ */
02520
02521 /* Subroutine errclb */
02522
02523 /* This subroutine checks the validity of the input data. */
02524
02525
02526 /* * * * */
02527
02528 /* NEOS, November 1994. (Latest revision June 1996.) */
02529 /* Optimization Technology Center. */
02530 /* Argonne National Laboratory and Northwestern University. */
02531 /* Written by */
02532 /* Ciyou Zhu */
02533 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02534
02535
02536 /* ************ */
02537 /* Check the input arguments for errors. */
02538 /* Parameter adjustments */
02539 --nbd;
02540 --u;
02541 --l;
02542
02543 /* Function Body */
02544 if (*n <= 0) {
02545 s_copy(task, "ERROR: N .LE. 0", (long int)60, (long int)15);
02546 }
02547 if (*m <= 0) {
02548 s_copy(task, "ERROR: M .LE. 0", (long int)60, (long int)15);
02549 }
02550 if (*factr < 0.) {
02551 s_copy(task, "ERROR: FACTR .LT. 0", (long int)60, (long int)19);
02552 }
02553 /* Check the validity of the arrays nbd(i), u(i), and l(i). */
02554 i__1 = *n;
02555 for (i__ = 1; i__ <= i__1; ++i__) {
02556 if (nbd[i__] < 0 || nbd[i__] > 3) {
02557 /* return */
02558 s_copy(task, "ERROR: INVALID NBD", (long int)60, (long int)18);
02559 *info = -6;
02560 *k = i__;
02561 }
02562 if (nbd[i__] == 2) {
02563 if (l[i__] > u[i__]) {
02564 /* return */
02565 s_copy(task, "ERROR: NO FEASIBLE SOLUTION", (long int)60, (
02566 long int)27);
02567 *info = -7;
02568 *k = i__;
02569 }
02570 }
02571 /* L10: */
02572 }
02573 return 0;
02574 } /* errclb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 2577 of file lbfgsb.cpp. References b, c__1, c__11, dcopy_(), ddot_(), dpofa_(), dtrsl_(), and t. Referenced by mainlb_(). 02586 {
02587 /* System generated locals */
02588 long int wn_dim1, wn_offset, wn1_dim1, wn1_offset, ws_dim1, ws_offset,
02589 wy_dim1, wy_offset, sy_dim1, sy_offset, i__1, i__2, i__3;
02590
02591 /* Local variables */
02592 static long int dend, pend;
02593 extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
02594 static long int upcl;
02595 static double temp1, temp2, temp3, temp4;
02596 static long int i__, k;
02597 extern /* Subroutine */ int dpofa_(double *a, long int *lda, long int *n, long int *info),
02598 dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy),
02599 dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
02600 static long int ipntr, jpntr, k1, m2, dbegin, is, js, iy, jy, pbegin, is1,
02601 js1, col2;
02602
02603 /* ************ */
02604
02605 /* Subroutine formk */
02606
02607 /* This subroutine forms the LEL^T factorization of the indefinite */
02608
02609 /* matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
02610 /* [L_a -R_z theta*S'AA'S ] */
02611 /* where E = [-I 0] */
02612 /* [ 0 I] */
02613 /* The matrix K can be shown to be equal to the matrix M^[-1]N */
02614 /* occurring in section 5.1 of [1], as well as to the matrix */
02615 /* Mbar^[-1] Nbar in section 5.3. */
02616
02617 /* n is an long int variable. */
02618 /* On entry n is the dimension of the problem. */
02619 /* On exit n is unchanged. */
02620
02621 /* nsub is an long int variable */
02622 /* On entry nsub is the number of subspace variables in free set. */
02623 /* On exit nsub is not changed. */
02624
02625 /* ind is an long int array of dimension nsub. */
02626 /* On entry ind specifies the indices of subspace variables. */
02627 /* On exit ind is unchanged. */
02628
02629 /* nenter is an long int variable. */
02630 /* On entry nenter is the number of variables entering the */
02631 /* free set. */
02632 /* On exit nenter is unchanged. */
02633
02634 /* ileave is an long int variable. */
02635 /* On entry indx2(ileave),...,indx2(n) are the variables leaving */
02636 /* the free set. */
02637 /* On exit ileave is unchanged. */
02638
02639 /* indx2 is an long int array of dimension n. */
02640 /* On entry indx2(1),...,indx2(nenter) are the variables entering */
02641 /* the free set, while indx2(ileave),...,indx2(n) are the */
02642 /* variables leaving the free set. */
02643 /* On exit indx2 is unchanged. */
02644
02645 /* iupdat is an long int variable. */
02646 /* On entry iupdat is the total number of BFGS updates made so far. */
02647 /* On exit iupdat is unchanged. */
02648
02649 /* updatd is a long int variable. */
02650 /* On entry 'updatd' is true if the L-BFGS matrix is updatd. */
02651 /* On exit 'updatd' is unchanged. */
02652
02653 /* wn is a double precision array of dimension 2m x 2m. */
02654 /* On entry wn is unspecified. */
02655 /* On exit the upper triangle of wn stores the LEL^T factorization */
02656 /* of the 2*col x 2*col indefinite matrix */
02657 /* [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
02658 /* [L_a -R_z theta*S'AA'S ] */
02659
02660 /* wn1 is a double precision array of dimension 2m x 2m. */
02661 /* On entry wn1 stores the lower triangular part of */
02662 /* [Y' ZZ'Y L_a'+R_z'] */
02663 /* [L_a+R_z S'AA'S ] */
02664 /* in the previous iteration. */
02665 /* On exit wn1 stores the corresponding updated matrices. */
02666 /* The purpose of wn1 is just to store these inner products */
02667 /* so they can be easily updated and inserted into wn. */
02668
02669 /* m is an long int variable. */
02670 /* On entry m is the maximum number of variable metric corrections */
02671 /* used to define the limited memory matrix. */
02672 /* On exit m is unchanged. */
02673
02674 /* ws, wy, sy, and wtyy are double precision arrays; */
02675 /* theta is a double precision variable; */
02676 /* col is an long int variable; */
02677 /* head is an long int variable. */
02678 /* On entry they store the information defining the */
02679 /* limited memory BFGS matrix: */
02680 /* ws(n,m) stores S, a set of s-vectors; */
02681 /* wy(n,m) stores Y, a set of y-vectors; */
02682 /* sy(m,m) stores S'Y; */
02683 /* wtyy(m,m) stores the Cholesky factorization */
02684 /* of (theta*S'S+LD^(-1)L') */
02685 /* theta is the scaling factor specifying B_0 = theta I; */
02686 /* col is the number of variable metric corrections stored; */
02687 /* head is the location of the 1st s- (or y-) vector in S (or Y). */
02688 /* On exit they are unchanged. */
02689
02690 /* info is an long int variable. */
02691 /* On entry info is unspecified. */
02692 /* On exit info = 0 for normal return; */
02693 /* = -1 when the 1st Cholesky factorization failed; */
02694 /* = -2 when the 2st Cholesky factorization failed. */
02695
02696 /* Subprograms called: */
02697
02698 /* Linpack ... dcopy, dpofa, dtrsl. */
02699
02700
02701 /* References: */
02702 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
02703 /* memory algorithm for bound constrained optimization'', */
02704 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
02705
02706 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
02707 /* limited memory FORTRAN code for solving bound constrained */
02708 /* optimization problems'', Tech. Report, NAM-11, EECS Department, */
02709 /* Northwestern University, 1994. */
02710
02711 /* (Postscript files of these papers are available via anonymous */
02712 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
02713
02714 /* * * * */
02715
02716 /* NEOS, November 1994. (Latest revision June 1996.) */
02717 /* Optimization Technology Center. */
02718 /* Argonne National Laboratory and Northwestern University. */
02719 /* Written by */
02720 /* Ciyou Zhu */
02721 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
02722
02723
02724 /* ************ */
02725 /* Form the lower triangular part of */
02726 /* WN1 = [Y' ZZ'Y L_a'+R_z'] */
02727 /* [L_a+R_z S'AA'S ] */
02728 /* where L_a is the strictly lower triangular part of S'AA'Y */
02729 /* R_z is the upper triangular part of S'ZZ'Y. */
02730 /* Parameter adjustments */
02731 --indx2;
02732 --ind;
02733 sy_dim1 = *m;
02734 sy_offset = 1 + sy_dim1 * 1;
02735 sy -= sy_offset;
02736 wy_dim1 = *n;
02737 wy_offset = 1 + wy_dim1 * 1;
02738 wy -= wy_offset;
02739 ws_dim1 = *n;
02740 ws_offset = 1 + ws_dim1 * 1;
02741 ws -= ws_offset;
02742 wn1_dim1 = 2 * *m;
02743 wn1_offset = 1 + wn1_dim1 * 1;
02744 wn1 -= wn1_offset;
02745 wn_dim1 = 2 * *m;
02746 wn_offset = 1 + wn_dim1 * 1;
02747 wn -= wn_offset;
02748
02749 /* Function Body */
02750 if (*updatd) {
02751 if (*iupdat > *m) {
02752 /* shift old part of WN1. */
02753 i__1 = *m - 1;
02754 for (jy = 1; jy <= i__1; ++jy) {
02755 js = *m + jy;
02756 i__2 = *m - jy;
02757 dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[
02758 jy + jy * wn1_dim1], &c__1);
02759 i__2 = *m - jy;
02760 dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[
02761 js + js * wn1_dim1], &c__1);
02762 i__2 = *m - 1;
02763 dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[
02764 *m + 1 + jy * wn1_dim1], &c__1);
02765 /* L10: */
02766 }
02767 }
02768 /* put new rows in blocks (1,1), (2,1) and (2,2). */
02769 pbegin = 1;
02770 pend = *nsub;
02771 dbegin = *nsub + 1;
02772 dend = *n;
02773 iy = *col;
02774 is = *m + *col;
02775 ipntr = *head + *col - 1;
02776 if (ipntr > *m) {
02777 ipntr -= *m;
02778 }
02779 jpntr = *head;
02780 i__1 = *col;
02781 for (jy = 1; jy <= i__1; ++jy) {
02782 js = *m + jy;
02783 temp1 = 0.;
02784 temp2 = 0.;
02785 temp3 = 0.;
02786 /* compute element jy of row 'col' of Y'ZZ'Y */
02787 i__2 = pend;
02788 for (k = pbegin; k <= i__2; ++k) {
02789 k1 = ind[k];
02790 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02791 /* L15: */
02792 }
02793 /* compute elements jy of row 'col' of L_a and S'AA'S */
02794 i__2 = dend;
02795 for (k = dbegin; k <= i__2; ++k) {
02796 k1 = ind[k];
02797 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02798 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02799 /* L16: */
02800 }
02801 wn1[iy + jy * wn1_dim1] = temp1;
02802 wn1[is + js * wn1_dim1] = temp2;
02803 wn1[is + jy * wn1_dim1] = temp3;
02804 jpntr = jpntr % *m + 1;
02805 /* L20: */
02806 }
02807 /* put new column in block (2,1). */
02808 jy = *col;
02809 jpntr = *head + *col - 1;
02810 if (jpntr > *m) {
02811 jpntr -= *m;
02812 }
02813 ipntr = *head;
02814 i__1 = *col;
02815 for (i__ = 1; i__ <= i__1; ++i__) {
02816 is = *m + i__;
02817 temp3 = 0.;
02818 /* compute element i of column 'col' of R_z */
02819 i__2 = pend;
02820 for (k = pbegin; k <= i__2; ++k) {
02821 k1 = ind[k];
02822 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02823 /* L25: */
02824 }
02825 ipntr = ipntr % *m + 1;
02826 wn1[is + jy * wn1_dim1] = temp3;
02827 /* L30: */
02828 }
02829 upcl = *col - 1;
02830 } else {
02831 upcl = *col;
02832 }
02833 /* modify the old parts in blocks (1,1) and (2,2) due to changes */
02834 /* in the set of free variables. */
02835 ipntr = *head;
02836 i__1 = upcl;
02837 for (iy = 1; iy <= i__1; ++iy) {
02838 is = *m + iy;
02839 jpntr = *head;
02840 i__2 = iy;
02841 for (jy = 1; jy <= i__2; ++jy) {
02842 js = *m + jy;
02843 temp1 = 0.;
02844 temp2 = 0.;
02845 temp3 = 0.;
02846 temp4 = 0.;
02847 i__3 = *nenter;
02848 for (k = 1; k <= i__3; ++k) {
02849 k1 = indx2[k];
02850 temp1 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02851 temp2 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02852 /* L35: */
02853 }
02854 i__3 = *n;
02855 for (k = *ileave; k <= i__3; ++k) {
02856 k1 = indx2[k];
02857 temp3 += wy[k1 + ipntr * wy_dim1] * wy[k1 + jpntr * wy_dim1];
02858 temp4 += ws[k1 + ipntr * ws_dim1] * ws[k1 + jpntr * ws_dim1];
02859 /* L36: */
02860 }
02861 wn1[iy + jy * wn1_dim1] = wn1[iy + jy * wn1_dim1] + temp1 - temp3;
02862 wn1[is + js * wn1_dim1] = wn1[is + js * wn1_dim1] - temp2 + temp4;
02863 jpntr = jpntr % *m + 1;
02864 /* L40: */
02865 }
02866 ipntr = ipntr % *m + 1;
02867 /* L45: */
02868 }
02869 /* modify the old parts in block (2,1). */
02870 ipntr = *head;
02871 i__1 = *m + upcl;
02872 for (is = *m + 1; is <= i__1; ++is) {
02873 jpntr = *head;
02874 i__2 = upcl;
02875 for (jy = 1; jy <= i__2; ++jy) {
02876 temp1 = 0.;
02877 temp3 = 0.;
02878 i__3 = *nenter;
02879 for (k = 1; k <= i__3; ++k) {
02880 k1 = indx2[k];
02881 temp1 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02882 /* L50: */
02883 }
02884 i__3 = *n;
02885 for (k = *ileave; k <= i__3; ++k) {
02886 k1 = indx2[k];
02887 temp3 += ws[k1 + ipntr * ws_dim1] * wy[k1 + jpntr * wy_dim1];
02888 /* L51: */
02889 }
02890 if (is <= jy + *m) {
02891 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] + temp1 -
02892 temp3;
02893 } else {
02894 wn1[is + jy * wn1_dim1] = wn1[is + jy * wn1_dim1] - temp1 +
02895 temp3;
02896 }
02897 jpntr = jpntr % *m + 1;
02898 /* L55: */
02899 }
02900 ipntr = ipntr % *m + 1;
02901 /* L60: */
02902 }
02903 /* Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ] */
02904 /* [-L_a +R_z S'AA'S*theta] */
02905 m2 = *m << 1;
02906 i__1 = *col;
02907 for (iy = 1; iy <= i__1; ++iy) {
02908 is = *col + iy;
02909 is1 = *m + iy;
02910 i__2 = iy;
02911 for (jy = 1; jy <= i__2; ++jy) {
02912 js = *col + jy;
02913 js1 = *m + jy;
02914 wn[jy + iy * wn_dim1] = wn1[iy + jy * wn1_dim1] / *theta;
02915 wn[js + is * wn_dim1] = wn1[is1 + js1 * wn1_dim1] * *theta;
02916 /* L65: */
02917 }
02918 i__2 = iy - 1;
02919 for (jy = 1; jy <= i__2; ++jy) {
02920 wn[jy + is * wn_dim1] = -wn1[is1 + jy * wn1_dim1];
02921 /* L66: */
02922 }
02923 i__2 = *col;
02924 for (jy = iy; jy <= i__2; ++jy) {
02925 wn[jy + is * wn_dim1] = wn1[is1 + jy * wn1_dim1];
02926 /* L67: */
02927 }
02928 wn[iy + iy * wn_dim1] += sy[iy + iy * sy_dim1];
02929 /* L70: */
02930 }
02931 /* Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] */
02932 /* [(-L_a +R_z)L'^-1 S'AA'S*theta ] */
02933 /* first Cholesky factor (1,1) block of wn to get LL' */
02934 /* with L' stored in the upper triangle of wn. */
02935 dpofa_(&wn[wn_offset], &m2, col, info);
02936 if (*info != 0) {
02937 *info = -1;
02938 return 0;
02939 }
02940 /* then form L^-1(-L_a'+R_z') in the (1,2) block. */
02941 col2 = *col << 1;
02942 i__1 = col2;
02943 for (js = *col + 1; js <= i__1; ++js) {
02944 dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info);
02945 /* L71: */
02946 }
02947 /* Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the */
02948 /* upper triangle of (2,2) block of wn. */
02949 i__1 = col2;
02950 for (is = *col + 1; is <= i__1; ++is) {
02951 i__2 = col2;
02952 for (js = is; js <= i__2; ++js) {
02953 wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], &c__1,
02954 &wn[js * wn_dim1 + 1], &c__1);
02955 /* L74: */
02956 }
02957 /* L72: */
02958 }
02959 /* Cholesky factorization of (2,2) block of wn. */
02960 dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info);
02961 if (*info != 0) {
02962 *info = -2;
02963 return 0;
02964 }
02965 return 0;
02966 } /* formk_ */
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 2969 of file lbfgsb.cpp. References dpofa_(), min_, and theta. Referenced by mainlb_(). 02975 {
02976 /* System generated locals */
02977 long int wt_dim1, wt_offset, sy_dim1, sy_offset, ss_dim1, ss_offset, i__1,
02978 i__2, i__3;
02979
02980 /* Local variables */
02981 static double ddum;
02982 static long int i__, j, k;
02983 extern /* Subroutine */ int dpofa_(double *a, long int *lda, long int *n, long int *info);
02984 static long int k1;
02985
02986 /* ************ */
02987
02988 /* Subroutine formt */
02989
02990 /* This subroutine forms the upper half of the pos. def. and symm. */
02991 /* T = theta*SS + L*D^(-1)*L', stores T in the upper triangle */
02992 /* of the array wt, and performs the Cholesky factorization of T */
02993 /* to produce J*J', with J' stored in the upper triangle of wt. */
02994
02995 /* Subprograms called: */
02996
02997 /* Linpack ... dpofa. */
02998
02999
03000 /* * * * */
03001
03002 /* NEOS, November 1994. (Latest revision June 1996.) */
03003 /* Optimization Technology Center. */
03004 /* Argonne National Laboratory and Northwestern University. */
03005 /* Written by */
03006 /* Ciyou Zhu */
03007 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03008
03009
03010 /* ************ */
03011 /* Form the upper half of T = theta*SS + L*D^(-1)*L', */
03012 /* store T in the upper triangle of the array wt. */
03013 /* Parameter adjustments */
03014 ss_dim1 = *m;
03015 ss_offset = 1 + ss_dim1 * 1;
03016 ss -= ss_offset;
03017 sy_dim1 = *m;
03018 sy_offset = 1 + sy_dim1 * 1;
03019 sy -= sy_offset;
03020 wt_dim1 = *m;
03021 wt_offset = 1 + wt_dim1 * 1;
03022 wt -= wt_offset;
03023
03024 /* Function Body */
03025 i__1 = *col;
03026 for (j = 1; j <= i__1; ++j) {
03027 wt[j * wt_dim1 + 1] = *theta * ss[j * ss_dim1 + 1];
03028 /* L52: */
03029 }
03030 i__1 = *col;
03031 for (i__ = 2; i__ <= i__1; ++i__) {
03032 i__2 = *col;
03033 for (j = i__; j <= i__2; ++j) {
03034 k1 = min_(i__,j) - 1;
03035 ddum = 0.;
03036 i__3 = k1;
03037 for (k = 1; k <= i__3; ++k) {
03038 ddum += sy[i__ + k * sy_dim1] * sy[j + k * sy_dim1] / sy[k +
03039 k * sy_dim1];
03040 /* L53: */
03041 }
03042 wt[i__ + j * wt_dim1] = ddum + *theta * ss[i__ + j * ss_dim1];
03043 /* L54: */
03044 }
03045 /* L55: */
03046 }
03047 /* Cholesky factorize T to J*J' with */
03048 /* J' stored in the upper triangle of wt. */
03049 dpofa_(&wt[wt_offset], m, col, info);
03050 if (*info != 0) {
03051 *info = -3;
03052 }
03053 return 0;
03054 } /* formt_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 3057 of file lbfgsb.cpp. Referenced by mainlb_(). 03062 {
03063 /* System generated locals */
03064 long int i__1;
03065
03066 /* Builtin functions */
03067 // long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03068
03069 /* Local variables */
03070 static long int iact, i__, k;
03071
03072 /* Fortran I/O blocks */
03073 /* static cilist io___168 = { 0, 6, 0, 0, 0 };
03074 static cilist io___169 = { 0, 6, 0, 0, 0 };
03075 static cilist io___170 = { 0, 6, 0, 0, 0 };
03076 static cilist io___172 = { 0, 6, 0, 0, 0 };*/
03077
03078
03079 /* ************ */
03080
03081 /* Subroutine freev */
03082
03083 /* This subroutine counts the entering and leaving variables when */
03084 /* iter > 0, and finds the index set of free and active variables */
03085 /* at the GCP. */
03086
03087 /* cnstnd is a long int variable indicating whether bounds are present */
03088
03089 /* index is an long int array of dimension n */
03090 /* for i=1,...,nfree, index(i) are the indices of free variables */
03091 /* for i=nfree+1,...,n, index(i) are the indices of bound variables */
03092 /* On entry after the first iteration, index gives */
03093 /* the free variables at the previous iteration. */
03094 /* On exit it gives the free variables based on the determination */
03095 /* in cauchy using the array iwhere. */
03096
03097 /* indx2 is an long int array of dimension n */
03098 /* On entry indx2 is unspecified. */
03099 /* On exit with iter>0, indx2 indicates which variables */
03100 /* have changed status since the previous iteration. */
03101 /* For i= 1,...,nenter, indx2(i) have changed from bound to free. */
03102 /* For i= ileave+1,...,n, indx2(i) have changed from free to bound. */
03103
03104
03105 /* * * * */
03106
03107 /* NEOS, November 1994. (Latest revision June 1996.) */
03108 /* Optimization Technology Center. */
03109 /* Argonne National Laboratory and Northwestern University. */
03110 /* Written by */
03111 /* Ciyou Zhu */
03112 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03113
03114
03115 /* ************ */
03116 /* Parameter adjustments */
03117 --iwhere;
03118 --indx2;
03119 --index;
03120
03121 /* Function Body */
03122 *nenter = 0;
03123 *ileave = *n + 1;
03124 if (*iter > 0 && *cnstnd) {
03125 /* count the entering and leaving variables. */
03126 i__1 = *nfree;
03127 for (i__ = 1; i__ <= i__1; ++i__) {
03128 k = index[i__];
03129 if (iwhere[k] > 0) {
03130 --(*ileave);
03131 indx2[*ileave] = k;
03132 if (*iprint >= 100) {
03133 /* s_wsle(&io___168);
03134 do_lio(&c__9, &c__1, "Variable ", (long int)9);
03135 do_lio(&c__3, &c__1, (char *)&k, (long int)sizeof(long int));
03136 do_lio(&c__9, &c__1, " leaves the set of free variables",
03137 (long int)33);
03138 e_wsle();*/
03139 }
03140 }
03141 /* L20: */
03142 }
03143 i__1 = *n;
03144 for (i__ = *nfree + 1; i__ <= i__1; ++i__) {
03145 k = index[i__];
03146 if (iwhere[k] <= 0) {
03147 ++(*nenter);
03148 indx2[*nenter] = k;
03149 if (*iprint >= 100) {
03150 /* s_wsle(&io___169);
03151 do_lio(&c__9, &c__1, "Variable ", (long int)9);
03152 do_lio(&c__3, &c__1, (char *)&k, (long int)sizeof(long int));
03153 do_lio(&c__9, &c__1, " enters the set of free variables",
03154 (long int)33);
03155 e_wsle();*/
03156 }
03157 }
03158 /* L22: */
03159 }
03160 if (*iprint >= 99) {
03161 /* s_wsle(&io___170);
03162 i__1 = *n + 1 - *ileave;
03163 do_lio(&c__3, &c__1, (char *)&i__1, (long int)sizeof(long int));
03164 do_lio(&c__9, &c__1, " variables leave; ", (long int)18);
03165 do_lio(&c__3, &c__1, (char *)&(*nenter), (long int)sizeof(long int));
03166 do_lio(&c__9, &c__1, " variables enter", (long int)16);
03167 e_wsle();*/
03168 }
03169 }
03170 *wrk = *ileave < *n + 1 || *nenter > 0 || *updatd;
03171 /* Find the index set of free and active variables at the GCP. */
03172 *nfree = 0;
03173 iact = *n + 1;
03174 i__1 = *n;
03175 for (i__ = 1; i__ <= i__1; ++i__) {
03176 if (iwhere[i__] <= 0) {
03177 ++(*nfree);
03178 index[*nfree] = i__;
03179 } else {
03180 --iact;
03181 index[iact] = i__;
03182 }
03183 /* L24: */
03184 }
03185 if (*iprint >= 99) {
03186 /* s_wsle(&io___172);
03187 do_lio(&c__3, &c__1, (char *)&(*nfree), (long int)sizeof(long int));
03188 do_lio(&c__9, &c__1, " variables are free at GCP ", (long int)27);
03189 i__1 = *iter + 1;
03190 do_lio(&c__3, &c__1, (char *)&i__1, (long int)sizeof(long int));
03191 e_wsle();*/
03192 }
03193 return 0;
03194 } /* freev_ */
|
|
||||||||||||||||||||
|
Definition at line 3197 of file lbfgsb.cpp. References t. Referenced by cauchy_(). 03201 {
03202 /* System generated locals */
03203 long int i__1;
03204
03205 /* Local variables */
03206 static double ddum;
03207 static long int i__, j, k, indxin, indxou;
03208 static double out;
03209
03210 /* ************ */
03211
03212 /* Subroutine hpsolb */
03213
03214 /* This subroutine sorts out the least element of t, and puts the */
03215 /* remaining elements of t in a heap. */
03216
03217 /* n is an long int variable. */
03218 /* On entry n is the dimension of the arrays t and iorder. */
03219 /* On exit n is unchanged. */
03220
03221 /* t is a double precision array of dimension n. */
03222 /* On entry t stores the elements to be sorted, */
03223 /* On exit t(n) stores the least elements of t, and t(1) to t(n-1) */
03224 /* stores the remaining elements in the form of a heap. */
03225
03226 /* iorder is an long int array of dimension n. */
03227 /* On entry iorder(i) is the index of t(i). */
03228 /* On exit iorder(i) is still the index of t(i), but iorder may be */
03229 /* permuted in accordance with t. */
03230
03231 /* iheap is an long int variable specifying the task. */
03232 /* On entry iheap should be set as follows: */
03233 /* iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, */
03234 /* iheap .ne. 0 if otherwise. */
03235 /* On exit iheap is unchanged. */
03236
03237
03238 /* References: */
03239 /* Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT. */
03240
03241 /* * * * */
03242
03243 /* NEOS, November 1994. (Latest revision June 1996.) */
03244 /* Optimization Technology Center. */
03245 /* Argonne National Laboratory and Northwestern University. */
03246 /* Written by */
03247 /* Ciyou Zhu */
03248 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03249
03250 /* ************ */
03251 /* Parameter adjustments */
03252 --iorder;
03253 --t;
03254
03255 /* Function Body */
03256 if (*iheap == 0) {
03257 /* Rearrange the elements t(1) to t(n) to form a heap. */
03258 i__1 = *n;
03259 for (k = 2; k <= i__1; ++k) {
03260 ddum = t[k];
03261 indxin = iorder[k];
03262 /* Add ddum to the heap. */
03263 i__ = k;
03264 L10:
03265 if (i__ > 1) {
03266 j = i__ / 2;
03267 if (ddum < t[j]) {
03268 t[i__] = t[j];
03269 iorder[i__] = iorder[j];
03270 i__ = j;
03271 goto L10;
03272 }
03273 }
03274 t[i__] = ddum;
03275 iorder[i__] = indxin;
03276 /* L20: */
03277 }
03278 }
03279 /* Assign to 'out' the value of t(1), the least member of the heap, */
03280 /* and rearrange the remaining members to form a heap as */
03281 /* elements 1 to n-1 of t. */
03282 if (*n > 1) {
03283 i__ = 1;
03284 out = t[1];
03285 indxou = iorder[1];
03286 ddum = t[*n];
03287 indxin = iorder[*n];
03288 /* Restore the heap */
03289 L30:
03290 j = i__ + i__;
03291 if (j <= *n - 1) {
03292 if (t[j + 1] < t[j]) {
03293 ++j;
03294 }
03295 if (t[j] < ddum) {
03296 t[i__] = t[j];
03297 iorder[i__] = iorder[j];
03298 i__ = j;
03299 goto L30;
03300 }
03301 }
03302 t[i__] = ddum;
03303 iorder[i__] = indxin;
03304 /* Put the least member in t(n). */
03305 t[*n] = out;
03306 iorder[*n] = indxou;
03307 }
03308 return 0;
03309 } /* hpsolb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 3312 of file lbfgsb.cpp. References c__1, c_b275, c_b276, c_b277, c_b9, dcopy_(), dcsrch_(), ddot_(), min_, s_cmp(), s_copy(), sqrt(), t, and x. Referenced by mainlb_(). 03330 {
03331 /* System generated locals */
03332 long int i__1;
03333 double d__1;
03334
03335 /* Builtin functions */
03336 long int s_cmp(char *, const char *const, long int, long int);
03337 //double sqrt();
03338 /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
03339
03340 /* Local variables */
03341 extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03342 static long int i__;
03343 extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03344 static double a1, a2;
03345 extern /* Subroutine */ int dcsrch_(double *f, double *g, double *stp, double *ftol, double *gtol, double *xtol,
03346 double *stpmin, double *stpmax, char *task, long int *isave, double *dsave, long int task_len);
03347
03348 /* ********** */
03349
03350 /* Subroutine lnsrlb */
03351
03352 /* This subroutine calls subroutine dcsrch from the Minpack2 library */
03353 /* to perform the line search. Subroutine dscrch is safeguarded so */
03354 /* that all trial points lie within the feasible region. */
03355
03356 /* Subprograms called: */
03357
03358 /* Minpack2 Library ... dcsrch. */
03359
03360 /* Linpack ... dtrsl, ddot. */
03361
03362
03363 /* * * * */
03364
03365 /* NEOS, November 1994. (Latest revision June 1996.) */
03366 /* Optimization Technology Center. */
03367 /* Argonne National Laboratory and Northwestern University. */
03368 /* Written by */
03369 /* Ciyou Zhu */
03370 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03371
03372
03373 /* ********** */
03374 /* Parameter adjustments */
03375 --z__;
03376 --t;
03377 --r__;
03378 --d__;
03379 --g;
03380 --x;
03381 --nbd;
03382 --u;
03383 --l;
03384 --isave;
03385 --dsave;
03386
03387 /* Function Body */
03388 if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
03389 goto L556;
03390 }
03391 *dtd = ddot_(n, &d__[1], &c__1, &d__[1], &c__1);
03392 *dnorm = sqrt(*dtd);
03393 /* Determine the maximum step length. */
03394 *stpmx = 1e10;
03395 if (*cnstnd) {
03396 if (*iter == 0) {
03397 *stpmx = 1.;
03398 } else {
03399 i__1 = *n;
03400 for (i__ = 1; i__ <= i__1; ++i__) {
03401 a1 = d__[i__];
03402 if (nbd[i__] != 0) {
03403 if (a1 < 0. && nbd[i__] <= 2) {
03404 a2 = l[i__] - x[i__];
03405 if (a2 >= 0.) {
03406 *stpmx = 0.;
03407 } else if (a1 * *stpmx < a2) {
03408 *stpmx = a2 / a1;
03409 }
03410 } else if (a1 > 0. && nbd[i__] >= 2) {
03411 a2 = u[i__] - x[i__];
03412 if (a2 <= 0.) {
03413 *stpmx = 0.;
03414 } else if (a1 * *stpmx > a2) {
03415 *stpmx = a2 / a1;
03416 }
03417 }
03418 }
03419 /* L43: */
03420 }
03421 }
03422 }
03423 if (*iter == 0 && ! (*boxed)) {
03424 /* Computing MIN */
03425 d__1 = 1. / *dnorm;
03426 *stp = min_(d__1,*stpmx);
03427 } else {
03428 *stp = 1.;
03429 }
03430 dcopy_(n, &x[1], &c__1, &t[1], &c__1);
03431 dcopy_(n, &g[1], &c__1, &r__[1], &c__1);
03432 *fold = *f;
03433 *ifun = 0;
03434 *iback = 0;
03435 s_copy(csave, "START", (long int)60, (long int)5);
03436 L556:
03437 *gd = ddot_(n, &g[1], &c__1, &d__[1], &c__1);
03438 if (*ifun == 0) {
03439 *gdold = *gd;
03440 if (*gd >= 0.) {
03441 /* the directional derivative >=0. */
03442 /* Line search is impossible. */
03443 *info = -4;
03444 return 0;
03445 }
03446 }
03447 dcsrch_(f, gd, stp, &c_b275, &c_b276, &c_b277, &c_b9, stpmx, csave, &
03448 isave[1], &dsave[1], (long int)60);
03449 *xstep = *stp * *dnorm;
03450 if (s_cmp(csave, "CONV", (long int)4, (long int)4) != 0 && s_cmp(csave, "WARN"
03451 , (long int)4, (long int)4) != 0) {
03452 s_copy(task, "FG_LNSRCH", (long int)60, (long int)9);
03453 ++(*ifun);
03454 ++(*nfgv);
03455 *iback = *ifun - 1;
03456 if (*stp == 1.) {
03457 dcopy_(n, &z__[1], &c__1, &x[1], &c__1);
03458 } else {
03459 i__1 = *n;
03460 for (i__ = 1; i__ <= i__1; ++i__) {
03461 x[i__] = *stp * d__[i__] + t[i__];
03462 /* L41: */
03463 }
03464 }
03465 } else {
03466 s_copy(task, "NEW_X", (long int)60, (long int)5);
03467 }
03468 return 0;
03469 } /* lnsrlb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 650 of file lbfgsb.cpp. References abs_, active_(), c__1, c_b9, cauchy_(), cmprlb_(), dcopy_(), ddot_(), dpmeps_(), dscal_(), errclb_(), formk_(), formt_(), freev_(), lnsrlb_(), matupd_(), max_, prn1lb_(), prn2lb_(), prn3lb_(), projgr_(), s_cmp(), s_copy(), subsm_(), t, theta, timer_(), v, and x. Referenced by setulb_(). 00671 {
00672 /* Format strings *//*
00673 static char fmt_1002[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
00674 .5,4x,\002|proj g|= \002,1p,d12.5)";
00675 static char fmt_1003[] = "(2(1x,i4),5x,\002-\002,5x,\002-\002,3x,\002\
00676 -\002,5x,\002-\002,5x,\002-\002,8x,\002-\002,3x,1p,2(1x,d10.3))";
00677 static char fmt_1001[] = "(//,\002ITERATION \002,i5)";
00678 static char fmt_1005[] = "(/,\002 Singular triangular system detected\
00679 ;\002,/,\002 refresh the lbfgs memory and restart the iteration.\002)";
00680 static char fmt_1006[] = "(/,\002 Nonpositive definiteness in Cholesky f\
00681 actorization in formk;\002,/,\002 refresh the lbfgs memory and restart the\
00682 iteration.\002)";
00683 static char fmt_1008[] = "(/,\002 Bad direction in the line search;\002,\
00684 /,\002 refresh the lbfgs memory and restart the iteration.\002)";
00685 static char fmt_1004[] = "(\002 ys=\002,1p,e10.3,\002 -gs=\002,1p,e10.\
00686 3,\002 BFGS update SKIPPED\002)";
00687 static char fmt_1007[] = "(/,\002 Nonpositive definiteness in Cholesky f\
00688 actorization in formt;\002,/,\002 refresh the lbfgs memory and restart the\
00689 iteration.\002)";*/
00690
00691 /* System generated locals */
00692 long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
00693 ss_dim1, ss_offset, yy_dim1, yy_offset, wt_dim1, wt_offset,
00694 wn_dim1, wn_offset, snd_dim1, snd_offset, i__1;
00695 double d__1, d__2;
00696 // olist o__1;
00697
00698 /* Builtin functions */
00699 long int s_cmp(char *, const char *const, long int, long int);
00700 /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
00701 //long int f_open(olist *), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
00702
00703 /* Local variables */
00704 static long int head;
00705 static double fold;
00706 static long int nact;
00707 static double ddum;
00708 extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
00709 static long int info;
00710 static double time;
00711 static long int nfgv, ifun, iter, nint;
00712 static char word[3];
00713 static double time1, time2;
00714 static long int i__, iback, k;
00715 extern /* Subroutine */ int dscal_(long int *n, double *da, double *dx, long int *incx);
00716 static double gdold;
00717 static long int nfree;
00718 static long int boxed;
00719 static long int itail;
00720 static double theta;
00721 extern /* Subroutine */ int freev_(long int *n, long int *nfree, long int *index, long int *nenter, long int *ileave, long int *indx2, long int *iwhere,
00722 long int *wrk, long int *updatd, long int *cnstnd, long int *iprint, long int *iter),
00723 dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
00724 static double dnorm;
00725 extern /* Subroutine */ int timer_(double *ttime),
00726 formk_(long int *n, long int *nsub, long int *ind, long int *nenter, long int *ileave, long int *indx2, long int *iupdat,
00727 long int *updatd, double *wn, double *wn1, long int *m, double *ws, double *wy, double *sy,
00728 double *theta, long int *col, long int *head, long int *info);
00729 static long int nskip, iword;
00730 extern /* Subroutine */ int formt_(long int *m, double *wt, double *sy, double *ss, long int *col, double *theta, long int *info),
00731 subsm_(long int *n, long int *m, long int *nsub, long int *ind, double *l, double *u, long int *nbd,
00732 double *x, double *d__, double *ws, double *wy, double *theta, long int *col,
00733 long int *head, long int *iword, double *wv, double *wn, long int *iprint, long int *info);
00734 static double xstep, stpmx;
00735 extern /* Subroutine */ int prn1lb_(long int*, long int*, double*, double*, double*, long int*, long int*, double*),
00736 prn2lb_(long int *n, double *x, double *f, double *g, long int *iprint, long int *itfile,
00737 long int *iter, long int *nfgv, long int *nact, double *sbgnrm, long int *nint, char *word,
00738 long int *iword, long int *iback, double *stp, double *xstep, long int word_len),
00739 prn3lb_(long int *n, double *x, double *f, char *task, long int *iprint, long int *info,
00740 long int *itfile, long int *iter, long int *nfgv, long int *nintol, long int *nskip, long int *nact,
00741 double *sbgnrm, double *time, long int *nint, char *word, long int *iback, double *stp,
00742 double *xstep, long int *k, double *cachyt, double *sbtime, double *lnscht,
00743 long int task_len, long int word_len);
00744 static double gd, dr, rr;
00745 static long int ileave;
00746 extern /* Subroutine */ int errclb_(long int *n, long int *m, double *factr, double *l, double *u, long int *nbd, char *task,
00747 long int *info, long int *k, long int task_len);
00748 static long int itfile;
00749 static double cachyt, epsmch;
00750 static long int updatd;
00751 static double sbtime;
00752 extern /* Subroutine */ int active_(long int *n, double *l, double *u, long int *nbd, double *x, long int *iwhere,
00753 long int *iprint, long int *prjctd, long int *cnstnd, long int *boxed);
00754 static long int prjctd;
00755 static long int iupdat;
00756 extern double dpmeps_();
00757 static long int cnstnd;
00758 static double sbgnrm;
00759 static long int nenter;
00760 static double lnscht;
00761 extern /* Subroutine */ int cauchy_(long int *n, double *x, double *l, double *u, long int *nbd, double *g, long int *iorder,
00762 long int *iwhere, double *t, double *d__, double *xcp, long int *m, double *wy, double *ws,
00763 double *sy, double *wt, double *theta, long int *col, long int *head, double *p, double *c__,
00764 double *wbp, double *v, long int *nint, double *sg, double *yg, long int *iprint, double *sbgnrm, long int *info),
00765 cmprlb_(long int *n, long int *m, double *x, double *g, double *ws, double *wy, double *sy,
00766 double *wt, double *z__, double *r__, double *wa, long int *index, double *theta,
00767 long int *col, long int *head, long int *nfree, long int *cnstnd, long int *info),
00768 lnsrlb_(long int *n, double *l, double *u, long int *nbd, double *x, double *f, double *fold,
00769 double *gd, double *gdold, double *g, double *d__, double *r__, double *t,
00770 double *z__, double *stp, double *dnorm, double *dtd, double *xstep, double *stpmx,
00771 long int *iter, long int *ifun, long int *iback, long int *nfgv, long int *info, char *task, long int *boxed,
00772 long int *cnstnd, char *csave, long int *isave, double *dsave, long int task_len, long int csave_len),
00773 matupd_(long int *n, long int *m, double *ws, double *wy, double *sy, double *ss,
00774 double *d__, double *r__, long int *itail, long int *iupdat, long int *col, long int *head,
00775 double *theta, double *rr, double *dr, double *stp, double *dtd);
00776 static long int nintol;
00777 extern /* Subroutine */ int projgr_(long int *n, double *l, double *u, long int *nbd, double *x, double *g, double *sbgnrm);
00778 static double dtd;
00779 static long int col;
00780 static double tol;
00781 static long int wrk;
00782 static double stp, cpu1, cpu2;
00783
00784 /* Fortran I/O blocks */
00785 /* static cilist io___62 = { 0, 6, 0, fmt_1002, 0 };
00786 static cilist io___63 = { 0, 0, 0, fmt_1003, 0 };
00787 static cilist io___64 = { 0, 6, 0, fmt_1001, 0 };
00788 static cilist io___66 = { 0, 6, 0, fmt_1005, 0 };
00789 static cilist io___68 = { 0, 6, 0, fmt_1006, 0 };
00790 static cilist io___69 = { 0, 6, 0, fmt_1005, 0 };
00791 static cilist io___71 = { 0, 6, 0, fmt_1008, 0 };
00792 static cilist io___75 = { 0, 6, 0, fmt_1004, 0 };
00793 static cilist io___76 = { 0, 6, 0, fmt_1007, 0 };*/
00794
00795
00796 /* ************ */
00797
00798 /* Subroutine mainlb */
00799
00800 /* This subroutine solves bound constrained optimization problems by */
00801 /* using the compact formula of the limited memory BFGS updates. */
00802
00803 /* n is an long int variable. */
00804 /* On entry n is the number of variables. */
00805 /* On exit n is unchanged. */
00806
00807 /* m is an long int variable. */
00808 /* On entry m is the maximum number of variable metric */
00809 /* corrections allowed in the limited memory matrix. */
00810 /* On exit m is unchanged. */
00811
00812 /* x is a double precision array of dimension n. */
00813 /* On entry x is an approximation to the solution. */
00814 /* On exit x is the current approximation. */
00815
00816 /* l is a double precision array of dimension n. */
00817 /* On entry l is the lower bound of x. */
00818 /* On exit l is unchanged. */
00819
00820 /* u is a double precision array of dimension n. */
00821 /* On entry u is the upper bound of x. */
00822 /* On exit u is unchanged. */
00823
00824 /* nbd is an long int array of dimension n. */
00825 /* On entry nbd represents the type of bounds imposed on the */
00826 /* variables, and must be specified as follows: */
00827 /* nbd(i)=0 if x(i) is unbounded, */
00828 /* 1 if x(i) has only a lower bound, */
00829 /* 2 if x(i) has both lower and upper bounds, */
00830 /* 3 if x(i) has only an upper bound. */
00831 /* On exit nbd is unchanged. */
00832
00833 /* f is a double precision variable. */
00834 /* On first entry f is unspecified. */
00835 /* On final exit f is the value of the function at x. */
00836
00837 /* g is a double precision array of dimension n. */
00838 /* On first entry g is unspecified. */
00839 /* On final exit g is the value of the gradient at x. */
00840
00841 /* factr is a double precision variable. */
00842 /* On entry factr >= 0 is specified by the user. The iteration */
00843 /* will stop when */
00844
00845 /* (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
00846
00847 /* where epsmch is the machine precision, which is automatically */
00848 /* generated by the code. */
00849 /* On exit factr is unchanged. */
00850
00851 /* pgtol is a double precision variable. */
00852 /* On entry pgtol >= 0 is specified by the user. The iteration */
00853 /* will stop when */
00854
00855 /* max{|proj g_i | i = 1, ..., n} <= pgtol */
00856
00857 /* where pg_i is the ith component of the projected gradient. */
00858 /* On exit pgtol is unchanged. */
00859
00860 /* ws, wy, sy, and wt are double precision working arrays used to */
00861 /* store the following information defining the limited memory */
00862 /* BFGS matrix: */
00863 /* ws, of dimension n x m, stores S, the matrix of s-vectors; */
00864 /* wy, of dimension n x m, stores Y, the matrix of y-vectors; */
00865 /* sy, of dimension m x m, stores S'Y; */
00866 /* ss, of dimension m x m, stores S'S; */
00867 /* yy, of dimension m x m, stores Y'Y; */
00868 /* wt, of dimension m x m, stores the Cholesky factorization */
00869 /* of (theta*S'S+LD^(-1)L'); see eq. */
00870 /* (2.26) in [3]. */
00871
00872 /* wn is a double precision working array of dimension 2m x 2m */
00873 /* used to store the LEL^T factorization of the indefinite matrix */
00874 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
00875 /* [L_a -R_z theta*S'AA'S ] */
00876
00877 /* where E = [-I 0] */
00878 /* [ 0 I] */
00879
00880 /* snd is a double precision working array of dimension 2m x 2m */
00881 /* used to store the lower triangular part of */
00882 /* N = [Y' ZZ'Y L_a'+R_z'] */
00883 /* [L_a +R_z S'AA'S ] */
00884
00885 /* z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays. */
00886 /* z is used at different times to store the Cauchy point and */
00887 /* the Newton point. */
00888
00889 /* sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays. */
00890
00891 /* index is an long int working array of dimension n. */
00892 /* In subroutine freev, index is used to store the free and fixed */
00893 /* variables at the Generalized Cauchy Point (GCP). */
00894
00895 /* iwhere is an long int working array of dimension n used to record */
00896 /* the status of the vector x for GCP computation. */
00897 /* iwhere(i)=0 or -3 if x(i) is free and has bounds, */
00898 /* 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) */
00899 /* 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) */
00900 /* 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) */
00901 /* -1 if x(i) is always free, i.e., no bounds on it. */
00902
00903 /* indx2 is an long int working array of dimension n. */
00904 /* Within subroutine cauchy, indx2 corresponds to the array iorder. */
00905 /* In subroutine freev, a list of variables entering and leaving */
00906 /* the free set is stored in indx2, and it is passed on to */
00907 /* subroutine formk with this information. */
00908
00909 /* task is a working string of characters of length 60 indicating */
00910 /* the current job when entering and leaving this subroutine. */
00911
00912 /* iprint is an long int variable that must be set by the user. */
00913 /* It controls the frequency and type of output generated: */
00914 /* iprint<0 no output is generated; */
00915 /* iprint=0 print only one line at the last iteration; */
00916 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
00917 /* iprint=99 print details of every iteration except n-vectors; */
00918 /* iprint=100 print also the changes of active set and final x; */
00919 /* iprint>100 print details of every iteration including x and g; */
00920 /* When iprint > 0, the file iterate.dat will be created to */
00921 /* summarize the iteration. */
00922
00923 /* csave is a working string of characters of length 60. */
00924
00925 /* lsave is a long int working array of dimension 4. */
00926
00927 /* isave is an long int working array of dimension 23. */
00928
00929 /* dsave is a double precision working array of dimension 29. */
00930
00931
00932 /* Subprograms called */
00933
00934 /* L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk, */
00935
00936 /* errclb, prn1lb, prn2lb, prn3lb, active, projgr, */
00937
00938 /* freev, cmprlb, matupd, formt. */
00939
00940 /* Minpack2 Library ... timer, dpmeps. */
00941
00942 /* Linpack Library ... dcopy, ddot. */
00943
00944
00945 /* References: */
00946
00947 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
00948 /* memory algorithm for bound constrained optimization'', */
00949 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
00950
00951 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN */
00952 /* Subroutines for Large Scale Bound Constrained Optimization'' */
00953 /* Tech. Report, NAM-11, EECS Department, Northwestern University, */
00954 /* 1994. */
00955
00956 /* [3] R. Byrd, J. Nocedal and R. Schnabel "Representations of */
00957 /* Quasi-Newton Matrices and their use in Limited Memory Methods'', */
00958 /* Mathematical Programming 63 (1994), no. 4, pp. 129-156. */
00959
00960 /* (Postscript files of these papers are available via anonymous */
00961 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
00962
00963 /* * * * */
00964
00965 /* NEOS, November 1994. (Latest revision June 1996.) */
00966 /* Optimization Technology Center. */
00967 /* Argonne National Laboratory and Northwestern University. */
00968 /* Written by */
00969 /* Ciyou Zhu */
00970 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
00971
00972
00973 /* ************ */
00974 /* Parameter adjustments */
00975 --indx2;
00976 --iwhere;
00977 --index;
00978 --t;
00979 --d__;
00980 --r__;
00981 --z__;
00982 --g;
00983 --nbd;
00984 --u;
00985 --l;
00986 --x;
00987 --ygo;
00988 --yg;
00989 --sgo;
00990 --sg;
00991 --wa;
00992 snd_dim1 = 2 * *m;
00993 snd_offset = 1 + snd_dim1 * 1;
00994 snd -= snd_offset;
00995 wn_dim1 = 2 * *m;
00996 wn_offset = 1 + wn_dim1 * 1;
00997 wn -= wn_offset;
00998 wt_dim1 = *m;
00999 wt_offset = 1 + wt_dim1 * 1;
01000 wt -= wt_offset;
01001 yy_dim1 = *m;
01002 yy_offset = 1 + yy_dim1 * 1;
01003 yy -= yy_offset;
01004 ss_dim1 = *m;
01005 ss_offset = 1 + ss_dim1 * 1;
01006 ss -= ss_offset;
01007 sy_dim1 = *m;
01008 sy_offset = 1 + sy_dim1 * 1;
01009 sy -= sy_offset;
01010 wy_dim1 = *n;
01011 wy_offset = 1 + wy_dim1 * 1;
01012 wy -= wy_offset;
01013 ws_dim1 = *n;
01014 ws_offset = 1 + ws_dim1 * 1;
01015 ws -= ws_offset;
01016 --lsave;
01017 --isave;
01018 --dsave;
01019
01020
01021 /* Function Body */
01022 if (s_cmp(task, "START", (long int)60, (long int)5) == 0) {
01023 timer_(&time1);
01024 /* Generate the current machine precision. */
01025 epsmch = dpmeps_();
01026 /* Initialize counters and scalars when task='START'. */
01027 /* for the limited memory BFGS matrices: */
01028 col = 0;
01029 head = 1;
01030 theta = 1.;
01031 iupdat = 0;
01032 updatd = FALSE_;
01033 /* for operation counts: */
01034 iter = 0;
01035 nfgv = 0;
01036 nint = 0;
01037 nintol = 0;
01038 nskip = 0;
01039 nfree = *n;
01040 /* for stopping tolerance: */
01041 tol = *factr * epsmch;
01042 /* for measuring running time: */
01043 cachyt = 0.;
01044 sbtime = 0.;
01045 lnscht = 0.;
01046 /* 'word' records the status of subspace solutions. */
01047 s_copy(word, "---", (long int)3, (long int)3);
01048 /* 'info' records the termination information. */
01049 info = 0;
01050 if (*iprint >= 1) {
01051 /* open a summary file 'iterate.dat' */
01052 /* o__1.oerr = 0;
01053 o__1.ounit = 8;
01054 o__1.ofnmlen = 11;
01055 o__1.ofnm = "iterate.dat";
01056 o__1.orl = 0;
01057 o__1.osta = "unknown";
01058 o__1.oacc = 0;
01059 o__1.ofm = 0;
01060 o__1.oblnk = 0;
01061 f_open(&o__1); */
01062 itfile = 8;
01063 }
01064 /* Check the input arguments for errors. */
01065 errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (long int)60);
01066 if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
01067 prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &
01068 nintol, &nskip, &nact, &sbgnrm, &c_b9, &nint, word, &
01069 iback, &stp, &xstep, &k, &cachyt, &sbtime, &lnscht, (
01070 long int)60, (long int)3);
01071 return 0;
01072 }
01073 prn1lb_(n, m, &l[1], &u[1], &x[1], iprint, &itfile, &epsmch);
01074 /* Initialize iwhere & project x onto the feasible set. */
01075 active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd,
01076 &cnstnd, &boxed);
01077 /* The end of the initialization. */
01078 } else {
01079 /* restore local variables. */
01080 prjctd = lsave[1];
01081 cnstnd = lsave[2];
01082 boxed = lsave[3];
01083 updatd = lsave[4];
01084 nintol = isave[1];
01085 itfile = isave[3];
01086 iback = isave[4];
01087 nskip = isave[5];
01088 head = isave[6];
01089 col = isave[7];
01090 itail = isave[8];
01091 iter = isave[9];
01092 iupdat = isave[10];
01093 nint = isave[12];
01094 nfgv = isave[13];
01095 info = isave[14];
01096 ifun = isave[15];
01097 iword = isave[16];
01098 nfree = isave[17];
01099 nact = isave[18];
01100 ileave = isave[19];
01101 nenter = isave[20];
01102 theta = dsave[1];
01103 fold = dsave[2];
01104 tol = dsave[3];
01105 dnorm = dsave[4];
01106 epsmch = dsave[5];
01107 cpu1 = dsave[6];
01108 cachyt = dsave[7];
01109 sbtime = dsave[8];
01110 lnscht = dsave[9];
01111 time1 = dsave[10];
01112 gd = dsave[11];
01113 stpmx = dsave[12];
01114 sbgnrm = dsave[13];
01115 stp = dsave[14];
01116 gdold = dsave[15];
01117 dtd = dsave[16];
01118 /* After returning from the driver go to the point where execution */
01119 /* is to resume. */
01120 if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
01121 goto L666;
01122 }
01123 if (s_cmp(task, "NEW_X", (long int)5, (long int)5) == 0) {
01124 goto L777;
01125 }
01126 if (s_cmp(task, "FG_ST", (long int)5, (long int)5) == 0) {
01127 goto L111;
01128 }
01129 if (s_cmp(task, "STOP", (long int)4, (long int)4) == 0) {
01130 if (s_cmp(task + 6, "CPU", (long int)3, (long int)3) == 0) {
01131 /* restore the previous iterate. */
01132 dcopy_(n, &t[1], &c__1, &x[1], &c__1);
01133 dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
01134 *f = fold;
01135 }
01136 goto L999;
01137 }
01138 }
01139 /* Compute f0 and g0. */
01140 s_copy(task, "FG_START", (long int)60, (long int)8);
01141 /* return to the driver to calculate f and g; reenter at 111. */
01142 goto L1000;
01143 L111:
01144 nfgv = 1;
01145 /* Compute the infinity norm of the (-) projected gradient. */
01146 projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
01147 if (*iprint >= 1) {
01148 /* s_wsfe(&io___62);
01149 do_fio(&c__1, (char *)&iter, (long int)sizeof(long int));
01150 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
01151 do_fio(&c__1, (char *)&sbgnrm, (long int)sizeof(double));
01152 e_wsfe();
01153 io___63.ciunit = itfile;
01154 s_wsfe(&io___63);
01155 do_fio(&c__1, (char *)&iter, (long int)sizeof(long int));
01156 do_fio(&c__1, (char *)&nfgv, (long int)sizeof(long int));
01157 do_fio(&c__1, (char *)&sbgnrm, (long int)sizeof(double));
01158 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
01159 e_wsfe();*/
01160 }
01161 if (sbgnrm <= *pgtol) {
01162 /* terminate the algorithm. */
01163 s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
01164 long int)60, (long int)48);
01165 goto L999;
01166 }
01167 /* ----------------- the beginning of the loop -------------------------- */
01168 L222:
01169 if (*iprint >= 99) {
01170 /* s_wsfe(&io___64);
01171 i__1 = iter + 1;
01172 do_fio(&c__1, (char *)&i__1, (long int)sizeof(long int));
01173 e_wsfe();*/
01174 }
01175 iword = -1;
01176
01177 if (! cnstnd && col > 0) {
01178 /* skip the search for GCP. */
01179 dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
01180 wrk = updatd;
01181 nint = 0;
01182 goto L333;
01183 }
01184 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01185
01186 /* Compute the Generalized Cauchy Point (GCP). */
01187
01188 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01189 timer_(&cpu1);
01190 cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[
01191 1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[
01192 sy_offset], &wt[wt_offset], &theta, &col, &head, &wa[1], &wa[(*m
01193 << 1) + 1], &wa[(*m << 2) + 1], &wa[*m * 6 + 1], &nint, &sg[1], &
01194 yg[1], iprint, &sbgnrm, &info);
01195 if (info != 0) {
01196 /* singular triangular system detected; refresh the lbfgs memory. */
01197 if (*iprint >= 1) {
01198 /* s_wsfe(&io___66);
01199 e_wsfe();*/
01200 }
01201 info = 0;
01202 col = 0;
01203 head = 1;
01204 theta = 1.;
01205 iupdat = 0;
01206 updatd = FALSE_;
01207 timer_(&cpu2);
01208 cachyt = cachyt + cpu2 - cpu1;
01209 goto L222;
01210 }
01211 timer_(&cpu2);
01212 cachyt = cachyt + cpu2 - cpu1;
01213 nintol += nint;
01214 /* Count the entering and leaving variables for iter > 0; */
01215 /* find the index set of free and active variables at the GCP. */
01216 freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], &
01217 wrk, &updatd, &cnstnd, iprint, &iter);
01218 nact = *n - nfree;
01219 L333:
01220 /* If there are no free variables or B=I, then */
01221 /* skip the subspace minimization. */
01222 if (nfree == 0 || col == 0) {
01223 goto L555;
01224 }
01225 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01226
01227 /* Subspace minimization. */
01228
01229 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01230 timer_(&cpu1);
01231 /* Form the LEL^T factorization of the indefinite */
01232 /* matrix K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
01233 /* [L_a -R_z theta*S'AA'S ] */
01234 /* where E = [-I 0] */
01235 /* [ 0 I] */
01236 if (wrk) {
01237 formk_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iupdat, &
01238 updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], &
01239 wy[wy_offset], &sy[sy_offset], &theta, &col, &head, &info);
01240 }
01241 if (info != 0) {
01242 /* nonpositive definiteness in Cholesky factorization; */
01243 /* refresh the lbfgs memory and restart the iteration. */
01244 if (*iprint >= 1) {
01245 /* s_wsfe(&io___68);
01246 e_wsfe();*/
01247 }
01248 info = 0;
01249 col = 0;
01250 head = 1;
01251 theta = 1.;
01252 iupdat = 0;
01253 updatd = FALSE_;
01254 timer_(&cpu2);
01255 sbtime = sbtime + cpu2 - cpu1;
01256 goto L222;
01257 }
01258 /* compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */
01259 /* from 'cauchy'). */
01260 cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset]
01261 , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], &theta, &
01262 col, &head, &nfree, &cnstnd, &info);
01263 if (info != 0) {
01264 goto L444;
01265 }
01266 /* call the direct method. */
01267 subsm_(n, m, &nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], &
01268 ws[ws_offset], &wy[wy_offset], &theta, &col, &head, &iword, &wa[1]
01269 , &wn[wn_offset], iprint, &info);
01270 L444:
01271 if (info != 0) {
01272 /* singular triangular system detected; */
01273 /* refresh the lbfgs memory and restart the iteration. */
01274 if (*iprint >= 1) {
01275 /* s_wsfe(&io___69);
01276 e_wsfe();*/
01277 }
01278 info = 0;
01279 col = 0;
01280 head = 1;
01281 theta = 1.;
01282 iupdat = 0;
01283 updatd = FALSE_;
01284 timer_(&cpu2);
01285 sbtime = sbtime + cpu2 - cpu1;
01286 goto L222;
01287 }
01288 timer_(&cpu2);
01289 sbtime = sbtime + cpu2 - cpu1;
01290 L555:
01291 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01292
01293 /* Line search and optimality tests. */
01294
01295 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01296 /* Generate the search direction d:=z-x. */
01297 i__1 = *n;
01298 for (i__ = 1; i__ <= i__1; ++i__) {
01299 d__[i__] = z__[i__] - x[i__];
01300 /* L40: */
01301 }
01302 timer_(&cpu1);
01303 L666:
01304 lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], &
01305 d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, &
01306 stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd,
01307 csave, &isave[22], &dsave[17], (long int)60, (long int)60);
01308 if (info != 0 || iback >= 20) {
01309 /* restore the previous iterate. */
01310 dcopy_(n, &t[1], &c__1, &x[1], &c__1);
01311 dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
01312 *f = fold;
01313 if (col == 0) {
01314 /* abnormal termination. */
01315 if (info == 0) {
01316 info = -9;
01317 /* restore the actual number of f and g evaluations etc. */
01318 --nfgv;
01319 --ifun;
01320 --iback;
01321 }
01322 s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (long int)60, (
01323 long int)30);
01324 ++iter;
01325 goto L999;
01326 } else {
01327 /* refresh the lbfgs memory and restart the iteration. */
01328 if (*iprint >= 1) {
01329 /*s_wsfe(&io___71);
01330 e_wsfe();*/
01331 }
01332 if (info == 0) {
01333 --nfgv;
01334 }
01335 info = 0;
01336 col = 0;
01337 head = 1;
01338 theta = 1.;
01339 iupdat = 0;
01340 updatd = FALSE_;
01341 s_copy(task, "RESTART_FROM_LNSRCH", (long int)60, (long int)19);
01342 timer_(&cpu2);
01343 lnscht = lnscht + cpu2 - cpu1;
01344 goto L222;
01345 }
01346 } else if (s_cmp(task, "FG_LN", (long int)5, (long int)5) == 0) {
01347 /* return to the driver for calculating f and g; reenter at 666. */
01348 goto L1000;
01349 } else {
01350 /* calculate and print out the quantities related to the new X. */
01351 timer_(&cpu2);
01352 lnscht = lnscht + cpu2 - cpu1;
01353 ++iter;
01354 /* Compute the infinity norm of the projected (-)gradient. */
01355 projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
01356 /* Print iteration information. */
01357 prn2lb_(n, &x[1], f, &g[1], iprint, &itfile, &iter, &nfgv, &nact, &
01358 sbgnrm, &nint, word, &iword, &iback, &stp, &xstep, (long int)3);
01359 goto L1000;
01360 }
01361 L777:
01362 /* Test for termination. */
01363 if (sbgnrm <= *pgtol) {
01364 /* terminate the algorithm. */
01365 s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
01366 long int)60, (long int)48);
01367 goto L999;
01368 }
01369 /* Computing MAX */
01370 d__1 = abs_(fold), d__2 = abs_(*f), d__1 = max_(d__1,d__2);
01371 ddum = max_(d__1,1.);
01372 if (fold - *f <= tol * ddum) {
01373 /* terminate the algorithm. */
01374 s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", (
01375 long int)60, (long int)47);
01376 if (iback >= 10) {
01377 info = -5;
01378 }
01379 /* i.e., to issue a warning if iback>10 in the line search. */
01380 goto L999;
01381 }
01382 /* Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
01383 i__1 = *n;
01384 for (i__ = 1; i__ <= i__1; ++i__) {
01385 r__[i__] = g[i__] - r__[i__];
01386 /* L42: */
01387 }
01388 rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
01389 if (stp == 1.) {
01390 dr = gd - gdold;
01391 ddum = -gdold;
01392 } else {
01393 dr = (gd - gdold) * stp;
01394 dscal_(n, &stp, &d__[1], &c__1);
01395 ddum = -gdold * stp;
01396 }
01397 if (dr <= epsmch * ddum) {
01398 /* skip the L-BFGS update. */
01399 ++nskip;
01400 updatd = FALSE_;
01401 if (*iprint >= 1) {
01402 /* s_wsfe(&io___75);
01403 do_fio(&c__1, (char *)&dr, (long int)sizeof(double));
01404 do_fio(&c__1, (char *)&ddum, (long int)sizeof(double));
01405 e_wsfe();*/
01406 }
01407 goto L888;
01408 }
01409 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01410
01411 /* Update the L-BFGS matrix. */
01412
01413 /* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
01414 updatd = TRUE_;
01415 ++iupdat;
01416 /* Update matrices WS and WY and form the middle matrix in B. */
01417 matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[
01418 ss_offset], &d__[1], &r__[1], &itail, &iupdat, &col, &head, &
01419 theta, &rr, &dr, &stp, &dtd);
01420 /* Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */
01421 /* Store T in the upper triangular of the array wt; */
01422 /* Cholesky factorize T to J*J' with */
01423 /* J' stored in the upper triangular of wt. */
01424 formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], &col, &theta, &
01425 info);
01426 if (info != 0) {
01427 /* nonpositive definiteness in Cholesky factorization; */
01428 /* refresh the lbfgs memory and restart the iteration. */
01429 if (*iprint >= 1) {
01430 /*s_wsfe(&io___76);
01431 e_wsfe();*/
01432 }
01433 info = 0;
01434 col = 0;
01435 head = 1;
01436 theta = 1.;
01437 iupdat = 0;
01438 updatd = FALSE_;
01439 goto L222;
01440 }
01441 /* Now the inverse of the middle matrix in B is */
01442 /* [ D^(1/2) O ] [ -D^(1/2) D^(-1/2)*L' ] */
01443 /* [ -L*D^(-1/2) J ] [ 0 J' ] */
01444 L888:
01445 /* -------------------- the end of the loop ----------------------------- */
01446 goto L222;
01447 L999:
01448 timer_(&time2);
01449 time = time2 - time1;
01450 prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &nintol,
01451 &nskip, &nact, &sbgnrm, &time, &nint, word, &iback, &stp, &xstep,
01452 &k, &cachyt, &sbtime, &lnscht, (long int)60, (long int)3);
01453 L1000:
01454 /* Save local variables. */
01455 lsave[1] = prjctd;
01456 lsave[2] = cnstnd;
01457 lsave[3] = boxed;
01458 lsave[4] = updatd;
01459 isave[1] = nintol;
01460 isave[3] = itfile;
01461 isave[4] = iback;
01462 isave[5] = nskip;
01463 isave[6] = head;
01464 isave[7] = col;
01465 isave[8] = itail;
01466 isave[9] = iter;
01467 isave[10] = iupdat;
01468 isave[12] = nint;
01469 isave[13] = nfgv;
01470 isave[14] = info;
01471 isave[15] = ifun;
01472 isave[16] = iword;
01473 isave[17] = nfree;
01474 isave[18] = nact;
01475 isave[19] = ileave;
01476 isave[20] = nenter;
01477 dsave[1] = theta;
01478 dsave[2] = fold;
01479 dsave[3] = tol;
01480 dsave[4] = dnorm;
01481 dsave[5] = epsmch;
01482 dsave[6] = cpu1;
01483 dsave[7] = cachyt;
01484 dsave[8] = sbtime;
01485 dsave[9] = lnscht;
01486 dsave[10] = time1;
01487 dsave[11] = gd;
01488 dsave[12] = stpmx;
01489 dsave[13] = sbgnrm;
01490 dsave[14] = stp;
01491 dsave[15] = gdold;
01492 dsave[16] = dtd;
01493 return 0;
01494 } /* mainlb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 3472 of file lbfgsb.cpp. References c__1, dcopy_(), ddot_(), and theta. Referenced by mainlb_(). 03479 {
03480 /* System generated locals */
03481 long int ws_dim1, ws_offset, wy_dim1, wy_offset, sy_dim1, sy_offset,
03482 ss_dim1, ss_offset, i__1, i__2;
03483
03484 /* Local variables */
03485 extern double ddot_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03486 static long int j;
03487 extern /* Subroutine */ int dcopy_(long int *n, double *dx, long int *incx, double *dy, long int *incy);
03488 static long int pointr;
03489
03490 /* ************ */
03491
03492 /* Subroutine matupd */
03493
03494 /* This subroutine updates matrices WS and WY, and forms the */
03495 /* middle matrix in B. */
03496
03497 /* Subprograms called: */
03498
03499 /* Linpack ... dcopy, ddot. */
03500
03501
03502 /* * * * */
03503
03504 /* NEOS, November 1994. (Latest revision June 1996.) */
03505 /* Optimization Technology Center. */
03506 /* Argonne National Laboratory and Northwestern University. */
03507 /* Written by */
03508 /* Ciyou Zhu */
03509 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03510
03511
03512 /* ************ */
03513 /* Set pointers for matrices WS and WY. */
03514 /* Parameter adjustments */
03515 --r__;
03516 --d__;
03517 ss_dim1 = *m;
03518 ss_offset = 1 + ss_dim1 * 1;
03519 ss -= ss_offset;
03520 sy_dim1 = *m;
03521 sy_offset = 1 + sy_dim1 * 1;
03522 sy -= sy_offset;
03523 wy_dim1 = *n;
03524 wy_offset = 1 + wy_dim1 * 1;
03525 wy -= wy_offset;
03526 ws_dim1 = *n;
03527 ws_offset = 1 + ws_dim1 * 1;
03528 ws -= ws_offset;
03529
03530 /* Function Body */
03531 if (*iupdat <= *m) {
03532 *col = *iupdat;
03533 *itail = (*head + *iupdat - 2) % *m + 1;
03534 } else {
03535 *itail = *itail % *m + 1;
03536 *head = *head % *m + 1;
03537 }
03538 /* Update matrices WS and WY. */
03539 dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1);
03540 dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1);
03541 /* Set theta=yy/ys. */
03542 *theta = *rr / *dr;
03543 /* Form the middle matrix in B. */
03544 /* update the upper triangle of SS, */
03545 /* and the lower triangle of SY: */
03546 if (*iupdat > *m) {
03547 /* move old information */
03548 i__1 = *col - 1;
03549 for (j = 1; j <= i__1; ++j) {
03550 dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1]
03551 , &c__1);
03552 i__2 = *col - j;
03553 dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j *
03554 sy_dim1], &c__1);
03555 /* L50: */
03556 }
03557 }
03558 /* add new information: the last row of SY */
03559 /* and the last column of SS: */
03560 pointr = *head;
03561 i__1 = *col - 1;
03562 for (j = 1; j <= i__1; ++j) {
03563 sy[*col + j * sy_dim1] = ddot_(n, &d__[1], &c__1, &wy[pointr *
03564 wy_dim1 + 1], &c__1);
03565 ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, &
03566 d__[1], &c__1);
03567 pointr = pointr % *m + 1;
03568 /* L51: */
03569 }
03570 if (*stp == 1.) {
03571 ss[*col + *col * ss_dim1] = *dtd;
03572 } else {
03573 ss[*col + *col * ss_dim1] = *stp * *stp * *dtd;
03574 }
03575 sy[*col + *col * sy_dim1] = *dr;
03576 return 0;
03577 } /* matupd_ */
|
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 3580 of file lbfgsb.cpp. Referenced by mainlb_(). 03585 {
03586 /* Format strings *//*
03587 static char fmt_7001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002 \
03588 * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
03589 static char fmt_2001[] = "(\002RUNNING THE L-BFGS-B CODE\002,/,/,\002it \
03590 = iteration number\002,/,\002nf = number of function evaluations\002,/,\
03591 \002nint = number of segments explored during the Cauchy search\002,/,\002n\
03592 act = number of active bounds at the generalized Cauchy point\002,/,\002sub\
03593 = manner in which the subspace minimization terminated:\002,/,\002 \
03594 con = converged, bnd = a bound was reached\002,/,\002itls = number of iter\
03595 ations performed in the line search\002,/,\002stepl = step length used\002,/,\
03596 \002tstep = norm of the displacement (total step)\002,/,\002projg = norm of \
03597 the projected gradient\002,/,\002f = function value\002,/,/,\002 \
03598 * * *\002,/,/,\002Machine precision =\002,1p,d10.3)";
03599 static char fmt_9001[] = "(/,3x,\002it\002,3x,\002nf\002,2x,\002nint\002\
03600 ,2x,\002nact\002,2x,\002sub\002,2x,\002itls\002,2x,\002stepl\002,4x,\002tstep\
03601 \002,5x,\002projg\002,8x,\002f\002)";
03602 static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";*/
03603
03604 /* System generated locals */
03605 // long int i__1;
03606
03607 /* Builtin functions */
03608 // long int s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe(), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03609
03610 /* Local variables */
03611 // static long int i__;
03612
03613 /* Fortran I/O blocks */
03614 /* static cilist io___185 = { 0, 6, 0, fmt_7001, 0 };
03615 static cilist io___186 = { 0, 6, 0, 0, 0 };
03616 static cilist io___187 = { 0, 0, 0, fmt_2001, 0 };
03617 static cilist io___188 = { 0, 0, 0, 0, 0 };
03618 static cilist io___189 = { 0, 0, 0, fmt_9001, 0 };
03619 static cilist io___190 = { 0, 6, 0, fmt_1004, 0 };
03620 static cilist io___192 = { 0, 6, 0, fmt_1004, 0 };
03621 static cilist io___193 = { 0, 6, 0, fmt_1004, 0 }; */
03622
03623
03624 /* ************ */
03625
03626 /* Subroutine prn1lb */
03627
03628 /* This subroutine prints the input data, initial point, upper and */
03629 /* lower bounds of each variable, machine precision, as well as */
03630 /* the headings of the output. */
03631
03632
03633 /* * * * */
03634
03635 /* NEOS, November 1994. (Latest revision June 1996.) */
03636 /* Optimization Technology Center. */
03637 /* Argonne National Laboratory and Northwestern University. */
03638 /* Written by */
03639 /* Ciyou Zhu */
03640 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03641
03642
03643 /* ************ */
03644 /* Parameter adjustments */
03645 --x;
03646 --u;
03647 --l;
03648
03649 /* Function Body */
03650 if (*iprint >= 0) {
03651 /* s_wsfe(&io___185);
03652 do_fio(&c__1, (char *)&(*epsmch), (long int)sizeof(double));
03653 e_wsfe();
03654 s_wsle(&io___186);
03655 do_lio(&c__9, &c__1, "N = ", (long int)4);
03656 do_lio(&c__3, &c__1, (char *)&(*n), (long int)sizeof(long int));
03657 do_lio(&c__9, &c__1, " M = ", (long int)8);
03658 do_lio(&c__3, &c__1, (char *)&(*m), (long int)sizeof(long int));
03659 e_wsle();
03660 if (*iprint >= 1) {
03661 io___187.ciunit = *itfile;
03662 s_wsfe(&io___187);
03663 do_fio(&c__1, (char *)&(*epsmch), (long int)sizeof(double));
03664 e_wsfe();
03665 io___188.ciunit = *itfile;
03666 s_wsle(&io___188);
03667 do_lio(&c__9, &c__1, "N = ", (long int)4);
03668 do_lio(&c__3, &c__1, (char *)&(*n), (long int)sizeof(long int));
03669 do_lio(&c__9, &c__1, " M = ", (long int)8);
03670 do_lio(&c__3, &c__1, (char *)&(*m), (long int)sizeof(long int));
03671 e_wsle();
03672 io___189.ciunit = *itfile;
03673 s_wsfe(&io___189);
03674 e_wsfe();
03675 if (*iprint > 100) {
03676 s_wsfe(&io___190);
03677 do_fio(&c__1, "L =", (long int)3);
03678 i__1 = *n;
03679 for (i__ = 1; i__ <= i__1; ++i__) {
03680 do_fio(&c__1, (char *)&l[i__], (long int)sizeof(double))
03681 ;
03682 }
03683 e_wsfe();
03684 s_wsfe(&io___192);
03685 do_fio(&c__1, "X0 =", (long int)4);
03686 i__1 = *n;
03687 for (i__ = 1; i__ <= i__1; ++i__) {
03688 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double))
03689 ;
03690 }
03691 e_wsfe();
03692 s_wsfe(&io___193);
03693 do_fio(&c__1, "U =", (long int)3);
03694 i__1 = *n;
03695 for (i__ = 1; i__ <= i__1; ++i__) {
03696 do_fio(&c__1, (char *)&u[i__], (long int)sizeof(double))
03697 ;
03698 }
03699 e_wsfe();
03700 }
03701 }*/
03702 }
03703 return 0;
03704 } /* prn1lb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 3707 of file lbfgsb.cpp. References s_copy(). Referenced by mainlb_(). 03719 {
03720 /* Format strings *//*
03721 static char fmt_2001[] = "(/,\002At iterate\002,i5,4x,\002f= \002,1p,d12\
03722 .5,4x,\002|proj g|= \002,1p,d12.5)";
03723 static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
03724 static char fmt_3001[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),1\
03725 p,2(1x,d10.3))";*/
03726
03727 /* System generated locals */
03728 // long int i__1;
03729
03730 /* Builtin functions */
03731 /* Subroutine */ int s_copy(char *, const char *const, long int, long int);
03732 // long int s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle(), s_wsfe(cilist *), do_fio(long int*, char*, long int), e_wsfe();
03733
03734 /* Local variables */
03735 // static long int imod, i__;
03736
03737 /* Fortran I/O blocks */
03738 /* static cilist io___194 = { 0, 6, 0, 0, 0 };
03739 static cilist io___195 = { 0, 6, 0, fmt_2001, 0 };
03740 static cilist io___196 = { 0, 6, 0, fmt_1004, 0 };
03741 static cilist io___198 = { 0, 6, 0, fmt_1004, 0 };
03742 static cilist io___200 = { 0, 6, 0, fmt_2001, 0 };
03743 static cilist io___201 = { 0, 0, 0, fmt_3001, 0 };*/
03744
03745
03746 /* ************ */
03747
03748 /* Subroutine prn2lb */
03749
03750 /* This subroutine prints out new information after a successful */
03751 /* line search. */
03752
03753
03754 /* * * * */
03755
03756 /* NEOS, November 1994. (Latest revision June 1996.) */
03757 /* Optimization Technology Center. */
03758 /* Argonne National Laboratory and Northwestern University. */
03759 /* Written by */
03760 /* Ciyou Zhu */
03761 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03762
03763
03764 /* ************ */
03765 /* 'word' records the status of subspace solutions. */
03766 /* Parameter adjustments */
03767 --g;
03768 --x;
03769
03770 /* Function Body */
03771 if (*iword == 0) {
03772 /* the subspace minimization converged. */
03773 s_copy(word, "con", (long int)3, (long int)3);
03774 } else if (*iword == 1) {
03775 /* the subspace minimization stopped at a bound. */
03776 s_copy(word, "bnd", (long int)3, (long int)3);
03777 } else if (*iword == 5) {
03778 /* the truncated Newton step has been used. */
03779 s_copy(word, "TNT", (long int)3, (long int)3);
03780 } else {
03781 s_copy(word, "---", (long int)3, (long int)3);
03782 }
03783 if (*iprint >= 99) {
03784 /* s_wsle(&io___194);
03785 do_lio(&c__9, &c__1, "LINE SEARCH", (long int)11);
03786 do_lio(&c__3, &c__1, (char *)&(*iback), (long int)sizeof(long int));
03787 do_lio(&c__9, &c__1, " times; norm of step = ", (long int)23);
03788 do_lio(&c__5, &c__1, (char *)&(*xstep), (long int)sizeof(double));
03789 e_wsle();
03790 s_wsfe(&io___195);
03791 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03792 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03793 do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03794 e_wsfe();
03795 if (*iprint > 100) {
03796 s_wsfe(&io___196);
03797 do_fio(&c__1, "X =", (long int)3);
03798 i__1 = *n;
03799 for (i__ = 1; i__ <= i__1; ++i__) {
03800 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
03801 }
03802 e_wsfe();
03803 s_wsfe(&io___198);
03804 do_fio(&c__1, "G =", (long int)3);
03805 i__1 = *n;
03806 for (i__ = 1; i__ <= i__1; ++i__) {
03807 do_fio(&c__1, (char *)&g[i__], (long int)sizeof(double));
03808 }
03809 e_wsfe();
03810 }*/
03811 } else if (*iprint > 0) {
03812 /* imod = *iter % *iprint;
03813 if (imod == 0) {
03814 s_wsfe(&io___200);
03815 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03816 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03817 do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03818 e_wsfe();
03819 }*/
03820 }
03821 if (*iprint >= 1) {
03822 /* io___201.ciunit = *itfile;
03823 s_wsfe(&io___201);
03824 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03825 do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
03826 do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
03827 do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
03828 do_fio(&c__1, word, (long int)3);
03829 do_fio(&c__1, (char *)&(*iback), (long int)sizeof(long int));
03830 do_fio(&c__1, (char *)&(*stp), (long int)sizeof(double));
03831 do_fio(&c__1, (char *)&(*xstep), (long int)sizeof(double));
03832 do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03833 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03834 e_wsfe();*/
03835 }
03836 return 0;
03837 } /* prn2lb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 3840 of file lbfgsb.cpp. References s_cmp(). Referenced by mainlb_(). 03858 {
03859 /* Format strings *//*
03860 static char fmt_3003[] = "(/,\002 * * *\002,/,/,\002Tit = to\
03861 tal number of iterations\002,/,\002Tnf = total number of function evaluati\
03862 ons\002,/,\002Tnint = total number of segments explored during\002,\002 Cauc\
03863 hy searches\002,/,\002Skip = number of BFGS updates skipped\002,/,\002Nact \
03864 = number of active bounds at final generalized\002,\002 Cauchy point\002,/\
03865 ,\002Projg = norm of the final projected gradient\002,/,\002F = final fu\
03866 nction value\002,/,/,\002 * * *\002)";
03867 static char fmt_3004[] = "(/,3x,\002N\002,3x,\002Tit\002,2x,\002Tnf\002,\
03868 2x,\002Tnint\002,2x,\002Skip\002,2x,\002Nact\002,5x,\002Projg\002,8x,\002\
03869 F\002)";
03870 static char fmt_3005[] = "(i5,2(1x,i4),(1x,i6),(2x,i4),(1x,i5),1p,2(2x,d\
03871 10.3))";
03872 static char fmt_1004[] = "(/,a4,1p,6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))";
03873 static char fmt_3009[] = "(/,a60)";
03874 static char fmt_9011[] = "(/,\002 Matrix in 1st Cholesky factorization i\
03875 n formk is not Pos. Def.\002)";
03876 static char fmt_9012[] = "(/,\002 Matrix in 2st Cholesky factorization i\
03877 n formk is not Pos. Def.\002)";
03878 static char fmt_9013[] = "(/,\002 Matrix in the Cholesky factorization i\
03879 n formt is not Pos. Def.\002)";
03880 static char fmt_9014[] = "(/,\002 Derivative >= 0, backtracking line sea\
03881 rch impossible.\002,/,\002 Previous x, f and g restored.\002,/,\002 Possib\
03882 le causes: 1 error in function or gradient evaluation;\002,/,\002 \
03883 2 rounding errors dominate computation.\002)";
03884 static char fmt_9015[] = "(/,\002 Warning: more than 10 function and gr\
03885 adient\002,/,\002 evaluations in the last line search. Termination\002,/\
03886 ,\002 may possibly be caused by a bad search direction.\002)";
03887 static char fmt_9018[] = "(/,\002 The triangular system is singular.\002)"
03888 ;
03889 static char fmt_9019[] = "(/,\002 Line search cannot locate an adequate \
03890 point after 20 function\002,/,\002 and gradient evaluations. Previous x, f\
03891 and g restored.\002,/,\002 Possible causes: 1 error in function or gradient\
03892 evaluation;\002,/,\002 2 rounding error dominate computati\
03893 on.\002)";
03894 static char fmt_3007[] = "(/,\002 Cauchy time\002,1p,e10.\
03895 3,\002 seconds.\002,/\002 Subspace minimization time\002,1p,e10.3,\002 secon\
03896 ds.\002,/\002 Line search time\002,1p,e10.3,\002 seconds.\002)";
03897 static char fmt_3008[] = "(/,\002 Total User time\002,1p,e10.3,\002 seco\
03898 nds.\002,/)";
03899 static char fmt_3002[] = "(2(1x,i4),2(1x,i5),2x,a3,1x,i4,1p,2(2x,d7.1),6\
03900 x,\002-\002,10x,\002-\002)";*/
03901
03902 /* System generated locals */
03903 // long int i__1;
03904
03905 /* Builtin functions */
03906 long int s_cmp(char *, const char *const, long int, long int);
03907 // long int s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
03908
03909 /* Local variables */
03910 // static long int i__;
03911
03912 /* Fortran I/O blocks */
03913 /* static cilist io___202 = { 0, 6, 0, fmt_3003, 0 };
03914 static cilist io___203 = { 0, 6, 0, fmt_3004, 0 };
03915 static cilist io___204 = { 0, 6, 0, fmt_3005, 0 };
03916 static cilist io___205 = { 0, 6, 0, fmt_1004, 0 };
03917 static cilist io___207 = { 0, 6, 0, 0, 0 };
03918 static cilist io___208 = { 0, 6, 0, fmt_3009, 0 };
03919 static cilist io___209 = { 0, 6, 0, fmt_9011, 0 };
03920 static cilist io___210 = { 0, 6, 0, fmt_9012, 0 };
03921 static cilist io___211 = { 0, 6, 0, fmt_9013, 0 };
03922 static cilist io___212 = { 0, 6, 0, fmt_9014, 0 };
03923 static cilist io___213 = { 0, 6, 0, fmt_9015, 0 };
03924 static cilist io___214 = { 0, 6, 0, 0, 0 };
03925 static cilist io___215 = { 0, 6, 0, 0, 0 };
03926 static cilist io___216 = { 0, 6, 0, fmt_9018, 0 };
03927 static cilist io___217 = { 0, 6, 0, fmt_9019, 0 };
03928 static cilist io___218 = { 0, 6, 0, fmt_3007, 0 };
03929 static cilist io___219 = { 0, 6, 0, fmt_3008, 0 };
03930 static cilist io___220 = { 0, 0, 0, fmt_3002, 0 };
03931 static cilist io___221 = { 0, 0, 0, fmt_3009, 0 };
03932 static cilist io___222 = { 0, 0, 0, fmt_9011, 0 };
03933 static cilist io___223 = { 0, 0, 0, fmt_9012, 0 };
03934 static cilist io___224 = { 0, 0, 0, fmt_9013, 0 };
03935 static cilist io___225 = { 0, 0, 0, fmt_9014, 0 };
03936 static cilist io___226 = { 0, 0, 0, fmt_9015, 0 };
03937 static cilist io___227 = { 0, 0, 0, fmt_9018, 0 };
03938 static cilist io___228 = { 0, 0, 0, fmt_9019, 0 };
03939 static cilist io___229 = { 0, 0, 0, fmt_3008, 0 };*/
03940
03941
03942 /* ************ */
03943
03944 /* Subroutine prn3lb */
03945
03946 /* This subroutine prints out information when either a built-in */
03947 /* convergence test is satisfied or when an error message is */
03948 /* generated. */
03949
03950
03951 /* * * * */
03952
03953 /* NEOS, November 1994. (Latest revision June 1996.) */
03954 /* Optimization Technology Center. */
03955 /* Argonne National Laboratory and Northwestern University. */
03956 /* Written by */
03957 /* Ciyou Zhu */
03958 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
03959
03960
03961 /* ************ */
03962 /* Parameter adjustments */
03963 --x;
03964
03965 /* Function Body */
03966 if (s_cmp(task, "ERROR", (long int)5, (long int)5) == 0) {
03967 goto L999;
03968 }
03969 if (*iprint >= 0) {
03970 /* s_wsfe(&io___202);
03971 e_wsfe();
03972 s_wsfe(&io___203);
03973 e_wsfe();
03974 s_wsfe(&io___204);
03975 do_fio(&c__1, (char *)&(*n), (long int)sizeof(long int));
03976 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
03977 do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
03978 do_fio(&c__1, (char *)&(*nintol), (long int)sizeof(long int));
03979 do_fio(&c__1, (char *)&(*nskip), (long int)sizeof(long int));
03980 do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
03981 do_fio(&c__1, (char *)&(*sbgnrm), (long int)sizeof(double));
03982 do_fio(&c__1, (char *)&(*f), (long int)sizeof(double));
03983 e_wsfe();
03984 if (*iprint >= 100) {
03985 s_wsfe(&io___205);
03986 do_fio(&c__1, "X =", (long int)3);
03987 i__1 = *n;
03988 for (i__ = 1; i__ <= i__1; ++i__) {
03989 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
03990 }
03991 e_wsfe();
03992 }
03993 if (*iprint >= 1) {
03994 s_wsle(&io___207);
03995 do_lio(&c__9, &c__1, " F =", (long int)4);
03996 do_lio(&c__5, &c__1, (char *)&(*f), (long int)sizeof(double));
03997 e_wsle();
03998 }*/
03999 }
04000 L999:
04001 if (*iprint >= 0) {
04002 /* s_wsfe(&io___208);
04003 do_fio(&c__1, task, (long int)60);
04004 e_wsfe();
04005 if (*info != 0) {
04006 if (*info == -1) {
04007 s_wsfe(&io___209);
04008 e_wsfe();
04009 }
04010 if (*info == -2) {
04011 s_wsfe(&io___210);
04012 e_wsfe();
04013 }
04014 if (*info == -3) {
04015 s_wsfe(&io___211);
04016 e_wsfe();
04017 }
04018 if (*info == -4) {
04019 s_wsfe(&io___212);
04020 e_wsfe();
04021 }
04022 if (*info == -5) {
04023 s_wsfe(&io___213);
04024 e_wsfe();
04025 }
04026 if (*info == -6) {
04027 s_wsle(&io___214);
04028 do_lio(&c__9, &c__1, " Input nbd(", (long int)11);
04029 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04030 do_lio(&c__9, &c__1, ") is invalid.", (long int)13);
04031 e_wsle();
04032 }
04033 if (*info == -7) {
04034 s_wsle(&io___215);
04035 do_lio(&c__9, &c__1, " l(", (long int)3);
04036 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04037 do_lio(&c__9, &c__1, ") > u(", (long int)6);
04038 do_lio(&c__3, &c__1, (char *)&(*k), (long int)sizeof(long int));
04039 do_lio(&c__9, &c__1, "). No feasible solution.", (long int)25);
04040 e_wsle();
04041 }
04042 if (*info == -8) {
04043 s_wsfe(&io___216);
04044 e_wsfe();
04045 }
04046 if (*info == -9) {
04047 s_wsfe(&io___217);
04048 e_wsfe();
04049 }
04050 }
04051 if (*iprint >= 1) {
04052 s_wsfe(&io___218);
04053 do_fio(&c__1, (char *)&(*cachyt), (long int)sizeof(double));
04054 do_fio(&c__1, (char *)&(*sbtime), (long int)sizeof(double));
04055 do_fio(&c__1, (char *)&(*lnscht), (long int)sizeof(double));
04056 e_wsfe();
04057 }
04058 s_wsfe(&io___219);
04059 do_fio(&c__1, (char *)&(*time), (long int)sizeof(double));
04060 e_wsfe();
04061 if (*iprint >= 1) {
04062 if (*info == -4 || *info == -9) {
04063 io___220.ciunit = *itfile;
04064 s_wsfe(&io___220);
04065 do_fio(&c__1, (char *)&(*iter), (long int)sizeof(long int));
04066 do_fio(&c__1, (char *)&(*nfgv), (long int)sizeof(long int));
04067 do_fio(&c__1, (char *)&(*nint), (long int)sizeof(long int));
04068 do_fio(&c__1, (char *)&(*nact), (long int)sizeof(long int));
04069 do_fio(&c__1, word, (long int)3);
04070 do_fio(&c__1, (char *)&(*iback), (long int)sizeof(long int));
04071 do_fio(&c__1, (char *)&(*stp), (long int)sizeof(double));
04072 do_fio(&c__1, (char *)&(*xstep), (long int)sizeof(double));
04073 e_wsfe();
04074 }
04075 io___221.ciunit = *itfile;
04076 s_wsfe(&io___221);
04077 do_fio(&c__1, task, (long int)60);
04078 e_wsfe();
04079 if (*info != 0) {
04080 if (*info == -1) {
04081 io___222.ciunit = *itfile;
04082 s_wsfe(&io___222);
04083 e_wsfe();
04084 }
04085 if (*info == -2) {
04086 io___223.ciunit = *itfile;
04087 s_wsfe(&io___223);
04088 e_wsfe();
04089 }
04090 if (*info == -3) {
04091 io___224.ciunit = *itfile;
04092 s_wsfe(&io___224);
04093 e_wsfe();
04094 }
04095 if (*info == -4) {
04096 io___225.ciunit = *itfile;
04097 s_wsfe(&io___225);
04098 e_wsfe();
04099 }
04100 if (*info == -5) {
04101 io___226.ciunit = *itfile;
04102 s_wsfe(&io___226);
04103 e_wsfe();
04104 }
04105 if (*info == -8) {
04106 io___227.ciunit = *itfile;
04107 s_wsfe(&io___227);
04108 e_wsfe();
04109 }
04110 if (*info == -9) {
04111 io___228.ciunit = *itfile;
04112 s_wsfe(&io___228);
04113 e_wsfe();
04114 }
04115 }
04116 io___229.ciunit = *itfile;
04117 s_wsfe(&io___229);
04118 do_fio(&c__1, (char *)&(*time), (long int)sizeof(double));
04119 e_wsfe();
04120 }*/
04121 }
04122 /* L3006: */
04123 return 0;
04124 } /* prn3lb_ */
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 4127 of file lbfgsb.cpp. References abs_, max_, min_, and x. Referenced by mainlb_(). 04132 {
04133 /* System generated locals */
04134 long int i__1;
04135 double d__1, d__2;
04136
04137 /* Local variables */
04138 static long int i__;
04139 static double gi;
04140
04141 /* ************ */
04142
04143 /* Subroutine projgr */
04144
04145 /* This subroutine computes the infinity norm of the projected */
04146 /* gradient. */
04147
04148
04149 /* * * * */
04150
04151 /* NEOS, November 1994. (Latest revision June 1996.) */
04152 /* Optimization Technology Center. */
04153 /* Argonne National Laboratory and Northwestern University. */
04154 /* Written by */
04155 /* Ciyou Zhu */
04156 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
04157
04158
04159 /* ************ */
04160 /* Parameter adjustments */
04161 --g;
04162 --x;
04163 --nbd;
04164 --u;
04165 --l;
04166
04167 /* Function Body */
04168 *sbgnrm = 0.;
04169 i__1 = *n;
04170 for (i__ = 1; i__ <= i__1; ++i__) {
04171 gi = g[i__];
04172 if (nbd[i__] != 0) {
04173 if (gi < 0.) {
04174 if (nbd[i__] >= 2) {
04175 /* Computing MAX */
04176 d__1 = x[i__] - u[i__];
04177 gi = max_(d__1,gi);
04178 }
04179 } else {
04180 if (nbd[i__] <= 2) {
04181 /* Computing MIN */
04182 d__1 = x[i__] - l[i__];
04183 gi = min_(d__1,gi);
04184 }
04185 }
04186 }
04187 /* Computing MAX */
04188 d__1 = *sbgnrm, d__2 = abs_(gi);
04189 *sbgnrm = max_(d__1,d__2);
04190 /* L15: */
04191 }
04192 return 0;
04193 } /* projgr_ */
|
|
||||||||||||||||||||
|
Definition at line 183 of file lbfgsb.cpp. Referenced by dcsrch_(), ilaenv_(), lnsrlb_(), mainlb_(), prn3lb_(), and setulb_(). 00183 {
00184 long int l;
00185 int i = 0;
00186 long int flag = 0;
00187
00188 l = min_(l1,l2);
00189
00190 while (i<l && flag==0) {
00191 if (str1[i] != str2[i]) { flag = 1; }
00192 i++;
00193 }
00194 return flag;
00195 }
|
|
||||||||||||||||||||
|
Definition at line 197 of file lbfgsb.cpp. References min_. Referenced by dcsrch_(), errclb_(), ilaenv_(), lnsrlb_(), mainlb_(), and prn2lb_(). 00197 {
00198 long int l;
00199 int i = 0;
00200 l = min_(l1,l2);
00201
00202 while (i<l) {
00203 str1[i]=str2[i];
00204 i++;
00205 }
00206 return 0;
00207 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 368 of file lbfgsb.cpp. References mainlb_(), s_cmp(), t, and x. Referenced by EMAN::Util::multi_align_error(), EMAN::Util::twoD_fine_ali(), EMAN::Util::twoD_fine_ali_G(), and EMAN::Util::twoD_to_3D_ali(). 00384 {
00385 /* System generated locals */
00386 long int i__1;
00387
00388 /* Builtin functions */
00389 long int s_cmp(char *, const char *const, long int, long int);
00390
00391 /* Local variables */
00392 static long int lsnd, lsgo, lygo, l1, l2, l3, ld, lr, lt;
00393 extern /* Subroutine */ int mainlb_(long int *n, long int *m, double *x, double *l, double *u, long int *nbd, double *f,
00394 double *g, double *factr, double *pgtol, double *ws, double *wy,
00395 double *sy, double *ss, double *yy, double *wt, double *wn, double *snd,
00396 double *z__, double *r__, double *d__, double *t, double *wa, double *sg,
00397 double *sgo, double *yg, double *ygo, long int *index, long int *iwhere, long int *indx2,
00398 char *task, long int *iprint, char *csave, long int *lsave, long int *isave, double *dsave,
00399 long int task_len, long int csave_len);
00400 static long int lz, lwa, lsg, lyg, lwn, lss, lws, lwt, lsy, lwy, lyy;
00401
00402 /* ************ */
00403
00404 /* Subroutine setulb */
00405
00406 /* This subroutine partitions the working arrays wa and iwa, and */
00407 /* then uses the limited memory BFGS method to solve the bound */
00408 /* constrained optimization problem by calling mainlb. */
00409 /* (The direct method will be used in the subspace minimization.) */
00410
00411 /* n is an long int variable. */
00412 /* On entry n is the dimension of the problem. */
00413 /* On exit n is unchanged. */
00414
00415 /* m is an long int variable. */
00416 /* On entry m is the maximum number of variable metric corrections */
00417 /* used to define the limited memory matrix. */
00418 /* On exit m is unchanged. */
00419
00420 /* x is a double precision array of dimension n. */
00421 /* On entry x is an approximation to the solution. */
00422 /* On exit x is the current approximation. */
00423
00424 /* l is a double precision array of dimension n. */
00425 /* On entry l is the lower bound on x. */
00426 /* On exit l is unchanged. */
00427
00428 /* u is a double precision array of dimension n. */
00429 /* On entry u is the upper bound on x. */
00430 /* On exit u is unchanged. */
00431
00432 /* nbd is an long int array of dimension n. */
00433 /* On entry nbd represents the type of bounds imposed on the */
00434 /* variables, and must be specified as follows: */
00435 /* nbd(i)=0 if x(i) is unbounded, */
00436 /* 1 if x(i) has only a lower bound, */
00437 /* 2 if x(i) has both lower and upper bounds, and */
00438 /* 3 if x(i) has only an upper bound. */
00439 /* On exit nbd is unchanged. */
00440
00441 /* f is a double precision variable. */
00442 /* On first entry f is unspecified. */
00443 /* On final exit f is the value of the function at x. */
00444
00445 /* g is a double precision array of dimension n. */
00446 /* On first entry g is unspecified. */
00447 /* On final exit g is the value of the gradient at x. */
00448
00449 /* factr is a double precision variable. */
00450 /* On entry factr >= 0 is specified by the user. The iteration */
00451 /* will stop when */
00452
00453 /* (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */
00454
00455 /* where epsmch is the machine precision, which is automatically */
00456 /* generated by the code. Typical values for factr: 1.d+12 for */
00457 /* low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely */
00458 /* high accuracy. */
00459 /* On exit factr is unchanged. */
00460
00461 /* pgtol is a double precision variable. */
00462 /* On entry pgtol >= 0 is specified by the user. The iteration */
00463 /* will stop when */
00464
00465 /* max{|proj g_i | i = 1, ..., n} <= pgtol */
00466
00467 /* where pg_i is the ith component of the projected gradient. */
00468 /* On exit pgtol is unchanged. */
00469
00470 /* wa is a double precision working array of length */
00471 /* (2mmax + 4)nmax + 12mmax^2 + 12mmax. */
00472
00473 /* iwa is an long int working array of length 3nmax. */
00474
00475 /* task is a working string of characters of length 60 indicating */
00476 /* the current job when entering and quitting this subroutine. */
00477
00478 /* iprint is an long int variable that must be set by the user. */
00479 /* It controls the frequency and type of output generated: */
00480 /* iprint<0 no output is generated; */
00481 /* iprint=0 print only one line at the last iteration; */
00482 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
00483 /* iprint=99 print details of every iteration except n-vectors; */
00484 /* iprint=100 print also the changes of active set and final x; */
00485 /* iprint>100 print details of every iteration including x and g; */
00486 /* When iprint > 0, the file iterate.dat will be created to */
00487 /* summarize the iteration. */
00488
00489 /* csave is a working string of characters of length 60. */
00490
00491 /* lsave is a long int working array of dimension 4. */
00492 /* On exit with 'task' = NEW_X, the following information is */
00493 /* available: */
00494 /* If lsave(1) = .true. then the initial X has been replaced by */
00495 /* its projection in the feasible set; */
00496 /* If lsave(2) = .true. then the problem is constrained; */
00497 /* If lsave(3) = .true. then each variable has upper and lower */
00498 /* bounds; */
00499
00500 /* isave is an long int working array of dimension 44. */
00501 /* On exit with 'task' = NEW_X, the following information is */
00502 /* available: */
00503 /* isave(22) = the total number of intervals explored in the */
00504 /* search of Cauchy points; */
00505 /* isave(26) = the total number of skipped BFGS updates before */
00506 /* the current iteration; */
00507 /* isave(30) = the number of current iteration; */
00508 /* isave(31) = the total number of BFGS updates prior the current */
00509 /* iteration; */
00510 /* isave(33) = the number of intervals explored in the search of */
00511 /* Cauchy point in the current iteration; */
00512 /* isave(34) = the total number of function and gradient */
00513 /* evaluations; */
00514 /* isave(36) = the number of function value or gradient */
00515 /* evaluations in the current iteration; */
00516 /* if isave(37) = 0 then the subspace argmin is within the box; */
00517 /* if isave(37) = 1 then the subspace argmin is beyond the box; */
00518 /* isave(38) = the number of free variables in the current */
00519 /* iteration; */
00520 /* isave(39) = the number of active constraints in the current */
00521 /* iteration; */
00522 /* n + 1 - isave(40) = the number of variables leaving the set of */
00523 /* active constraints in the current iteration; */
00524 /* isave(41) = the number of variables entering the set of active */
00525 /* constraints in the current iteration. */
00526
00527 /* dsave is a double precision working array of dimension 29. */
00528 /* On exit with 'task' = NEW_X, the following information is */
00529 /* available: */
00530 /* dsave(1) = current 'theta' in the BFGS matrix; */
00531 /* dsave(2) = f(x) in the previous iteration; */
00532 /* dsave(3) = factr*epsmch; */
00533 /* dsave(4) = 2-norm of the line search direction vector; */
00534 /* dsave(5) = the machine precision epsmch generated by the code; */
00535 /* dsave(7) = the accumulated time spent on searching for */
00536 /* Cauchy points; */
00537 /* dsave(8) = the accumulated time spent on */
00538 /* subspace minimization; */
00539 /* dsave(9) = the accumulated time spent on line search; */
00540 /* dsave(11) = the slope of the line search function at */
00541 /* the current point of line search; */
00542 /* dsave(12) = the maximum relative step length imposed in */
00543 /* line search; */
00544 /* dsave(13) = the infinity norm of the projected gradient; */
00545 /* dsave(14) = the relative step length in the line search; */
00546 /* dsave(15) = the slope of the line search function at */
00547 /* the starting point of the line search; */
00548 /* dsave(16) = the square of the 2-norm of the line search */
00549 /* direction vector. */
00550
00551 /* Subprograms called: */
00552
00553 /* L-BFGS-B Library ... mainlb. */
00554
00555
00556 /* References: */
00557
00558 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
00559 /* memory algorithm for bound constrained optimization'', */
00560 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
00561
00562 /* [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
00563 /* limited memory FORTRAN code for solving bound constrained */
00564 /* optimization problems'', Tech. Report, NAM-11, EECS Department, */
00565 /* Northwestern University, 1994. */
00566
00567 /* (Postscript files of these papers are available via anonymous */
00568 /* ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */
00569
00570 /* * * * */
00571
00572 /* NEOS, November 1994. (Latest revision June 1996.) */
00573 /* Optimization Technology Center. */
00574 /* Argonne National Laboratory and Northwestern University. */
00575 /* Written by */
00576 /* Ciyou Zhu */
00577 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
00578
00579
00580 /* ************ */
00581 /* Parameter adjustments */
00582 --iwa;
00583 --g;
00584 --nbd;
00585 --u;
00586 --l;
00587 --x;
00588 --wa;
00589 --lsave;
00590 --isave;
00591 --dsave;
00592
00593 /* Function Body */
00594 if (s_cmp(task, "START", (long int)60, (long int)5) == 0) {
00595 isave[1] = *m * *n;
00596 /* Computing 2nd power */
00597 i__1 = *m;
00598 isave[2] = i__1 * i__1;
00599 /* Computing 2nd power */
00600 i__1 = *m;
00601 isave[3] = i__1 * i__1 << 2;
00602 isave[4] = 1;
00603 isave[5] = isave[4] + isave[1];
00604 isave[6] = isave[5] + isave[1];
00605 isave[7] = isave[6] + isave[2];
00606 isave[8] = isave[7] + isave[2];
00607 isave[9] = isave[8] + isave[2];
00608 isave[10] = isave[9] + isave[2];
00609 isave[11] = isave[10] + isave[3];
00610 isave[12] = isave[11] + isave[3];
00611 isave[13] = isave[12] + *n;
00612 isave[14] = isave[13] + *n;
00613 isave[15] = isave[14] + *n;
00614 isave[16] = isave[15] + *n;
00615 isave[17] = isave[16] + (*m << 3);
00616 isave[18] = isave[17] + *m;
00617 isave[19] = isave[18] + *m;
00618 isave[20] = isave[19] + *m;
00619 }
00620 l1 = isave[1];
00621 l2 = isave[2];
00622 l3 = isave[3];
00623 lws = isave[4];
00624 lwy = isave[5];
00625 lsy = isave[6];
00626 lss = isave[7];
00627 lyy = isave[8];
00628 lwt = isave[9];
00629 lwn = isave[10];
00630 lsnd = isave[11];
00631 lz = isave[12];
00632 lr = isave[13];
00633 ld = isave[14];
00634 lt = isave[15];
00635 lwa = isave[16];
00636 lsg = isave[17];
00637 lsgo = isave[18];
00638 lyg = isave[19];
00639 lygo = isave[20];
00640 mainlb_(n, m, &x[1], &l[1], &u[1], &nbd[1], f, &g[1], factr, pgtol, &wa[
00641 lws], &wa[lwy], &wa[lsy], &wa[lss], &wa[lyy], &wa[lwt], &wa[lwn],
00642 &wa[lsnd], &wa[lz], &wa[lr], &wa[ld], &wa[lt], &wa[lwa], &wa[lsg],
00643 &wa[lsgo], &wa[lyg], &wa[lygo], &iwa[1], &iwa[*n + 1], &iwa[(*n
00644 << 1) + 1], task, iprint, csave, &lsave[1], &isave[22], &dsave[1],
00645 (long int)60, (long int)60);
00646 return 0;
00647 } /* setulb_ */
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 4196 of file lbfgsb.cpp. References b, c__1, c__11, dtrsl_(), t, theta, and x. Referenced by mainlb_(). 04206 {
04207 /* Format strings *//*
04208 static char fmt_1001[] = "(/,\002----------------SUBSM entered----------\
04209 -------\002,/)";
04210 static char fmt_1002[] = "(\002ALPHA = \002,f7.5,\002 backtrack to the B\
04211 OX\002)";
04212 static char fmt_1003[] = "(\002Subspace solution X = \002,/,(4x,1p,6(1x\
04213 ,d11.4)))";
04214 static char fmt_1004[] = "(/,\002----------------exit SUBSM ------------\
04215 --------\002,/)";*/
04216
04217 /* System generated locals */
04218 long int ws_dim1, ws_offset, wy_dim1, wy_offset, wn_dim1, wn_offset, i__1,
04219 i__2;
04220
04221 /* Builtin functions */
04222 // long int s_wsfe(cilist *), e_wsfe(), do_fio(long int*, char*, long int), s_wsle(cilist *), do_lio(long int *, long int *, char *, long int), e_wsle();
04223
04224 /* Local variables */
04225 static double temp1, temp2;
04226 static long int i__, j, k;
04227 static double alpha;
04228 //extern /* Subroutine */ int dtrsl_();
04229 extern int dtrsl_(double *t, long int *ldt, long int *n, double *b, long int *job, long int *info);
04230 static long int m2;
04231 static double dk;
04232 static long int js, jy, pointr, ibd, col2;
04233
04234 /* Fortran I/O blocks */
04235 /* static cilist io___232 = { 0, 6, 0, fmt_1001, 0 };
04236 static cilist io___246 = { 0, 6, 0, fmt_1002, 0 };
04237 static cilist io___247 = { 0, 6, 0, 0, 0 };
04238 static cilist io___248 = { 0, 6, 0, fmt_1003, 0 };
04239 static cilist io___249 = { 0, 6, 0, fmt_1004, 0 }; */
04240
04241
04242 /* ************ */
04243
04244 /* Subroutine subsm */
04245
04246 /* Given xcp, l, u, r, an index set that specifies */
04247 /* the active set at xcp, and an l-BFGS matrix B */
04248 /* (in terms of WY, WS, SY, WT, head, col, and theta), */
04249 /* this subroutine computes an approximate solution */
04250 /* of the subspace problem */
04251
04252 /* (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp) */
04253
04254 /* subject to l<=x<=u */
04255 /* x_i=xcp_i for all i in A(xcp) */
04256
04257 /* along the subspace unconstrained Newton direction */
04258
04259 /* d = -(Z'BZ)^(-1) r. */
04260
04261 /* The formula for the Newton direction, given the L-BFGS matrix */
04262 /* and the Sherman-Morrison formula, is */
04263
04264 /* d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r. */
04265
04266 /* where */
04267 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
04268 /* [L_a -R_z theta*S'AA'S ] */
04269
04270 /* Note that this procedure for computing d differs */
04271 /* from that described in [1]. One can show that the matrix K is */
04272 /* equal to the matrix M^[-1]N in that paper. */
04273
04274 /* n is an long int variable. */
04275 /* On entry n is the dimension of the problem. */
04276 /* On exit n is unchanged. */
04277
04278 /* m is an long int variable. */
04279 /* On entry m is the maximum number of variable metric corrections */
04280 /* used to define the limited memory matrix. */
04281 /* On exit m is unchanged. */
04282
04283 /* nsub is an long int variable. */
04284 /* On entry nsub is the number of free variables. */
04285 /* On exit nsub is unchanged. */
04286
04287 /* ind is an long int array of dimension nsub. */
04288 /* On entry ind specifies the coordinate indices of free variables. */
04289 /* On exit ind is unchanged. */
04290
04291 /* l is a double precision array of dimension n. */
04292 /* On entry l is the lower bound of x. */
04293 /* On exit l is unchanged. */
04294
04295 /* u is a double precision array of dimension n. */
04296 /* On entry u is the upper bound of x. */
04297 /* On exit u is unchanged. */
04298
04299 /* nbd is a long int array of dimension n. */
04300 /* On entry nbd represents the type of bounds imposed on the */
04301 /* variables, and must be specified as follows: */
04302 /* nbd(i)=0 if x(i) is unbounded, */
04303 /* 1 if x(i) has only a lower bound, */
04304 /* 2 if x(i) has both lower and upper bounds, and */
04305 /* 3 if x(i) has only an upper bound. */
04306 /* On exit nbd is unchanged. */
04307
04308 /* x is a double precision array of dimension n. */
04309 /* On entry x specifies the Cauchy point xcp. */
04310 /* On exit x(i) is the minimizer of Q over the subspace of */
04311 /* free variables. */
04312
04313 /* d is a double precision array of dimension n. */
04314 /* On entry d is the reduced gradient of Q at xcp. */
04315 /* On exit d is the Newton direction of Q. */
04316
04317 /* ws and wy are double precision arrays; */
04318 /* theta is a double precision variable; */
04319 /* col is an long int variable; */
04320 /* head is an long int variable. */
04321 /* On entry they store the information defining the */
04322 /* limited memory BFGS matrix: */
04323 /* ws(n,m) stores S, a set of s-vectors; */
04324 /* wy(n,m) stores Y, a set of y-vectors; */
04325 /* theta is the scaling factor specifying B_0 = theta I; */
04326 /* col is the number of variable metric corrections stored; */
04327 /* head is the location of the 1st s- (or y-) vector in S (or Y). */
04328 /* On exit they are unchanged. */
04329
04330 /* iword is an long int variable. */
04331 /* On entry iword is unspecified. */
04332 /* On exit iword specifies the status of the subspace solution. */
04333 /* iword = 0 if the solution is in the box, */
04334 /* 1 if some bound is encountered. */
04335
04336 /* wv is a double precision working array of dimension 2m. */
04337
04338 /* wn is a double precision array of dimension 2m x 2m. */
04339 /* On entry the upper triangle of wn stores the LEL^T factorization */
04340 /* of the indefinite matrix */
04341
04342 /* K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] */
04343 /* [L_a -R_z theta*S'AA'S ] */
04344 /* where E = [-I 0] */
04345 /* [ 0 I] */
04346 /* On exit wn is unchanged. */
04347
04348 /* iprint is an long int variable that must be set by the user. */
04349 /* It controls the frequency and type of output generated: */
04350 /* iprint<0 no output is generated; */
04351 /* iprint=0 print only one line at the last iteration; */
04352 /* 0<iprint<99 print also f and |proj g| every iprint iterations; */
04353 /* iprint=99 print details of every iteration except n-vectors; */
04354 /* iprint=100 print also the changes of active set and final x; */
04355 /* iprint>100 print details of every iteration including x and g; */
04356 /* When iprint > 0, the file iterate.dat will be created to */
04357 /* summarize the iteration. */
04358
04359 /* info is an long int variable. */
04360 /* On entry info is unspecified. */
04361 /* On exit info = 0 for normal return, */
04362 /* = nonzero for abnormal return */
04363 /* when the matrix K is ill-conditioned. */
04364
04365 /* Subprograms called: */
04366
04367 /* Linpack dtrsl. */
04368
04369
04370 /* References: */
04371
04372 /* [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
04373 /* memory algorithm for bound constrained optimization'', */
04374 /* SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */
04375
04376
04377
04378 /* * * * */
04379
04380 /* NEOS, November 1994. (Latest revision June 1996.) */
04381 /* Optimization Technology Center. */
04382 /* Argonne National Laboratory and Northwestern University. */
04383 /* Written by */
04384 /* Ciyou Zhu */
04385 /* in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */
04386
04387
04388 /* ************ */
04389 /* Parameter adjustments */
04390 --d__;
04391 --x;
04392 --nbd;
04393 --u;
04394 --l;
04395 wn_dim1 = 2 * *m;
04396 wn_offset = 1 + wn_dim1 * 1;
04397 wn -= wn_offset;
04398 --wv;
04399 wy_dim1 = *n;
04400 wy_offset = 1 + wy_dim1 * 1;
04401 wy -= wy_offset;
04402 ws_dim1 = *n;
04403 ws_offset = 1 + ws_dim1 * 1;
04404 ws -= ws_offset;
04405 --ind;
04406
04407 /* Function Body */
04408 if (*nsub <= 0) {
04409 return 0;
04410 }
04411 if (*iprint >= 99) {
04412 /* s_wsfe(&io___232);
04413 e_wsfe();*/
04414 }
04415 /* Compute wv = W'Zd. */
04416 pointr = *head;
04417 i__1 = *col;
04418 for (i__ = 1; i__ <= i__1; ++i__) {
04419 temp1 = 0.;
04420 temp2 = 0.;
04421 i__2 = *nsub;
04422 for (j = 1; j <= i__2; ++j) {
04423 k = ind[j];
04424 temp1 += wy[k + pointr * wy_dim1] * d__[j];
04425 temp2 += ws[k + pointr * ws_dim1] * d__[j];
04426 /* L10: */
04427 }
04428 wv[i__] = temp1;
04429 wv[*col + i__] = *theta * temp2;
04430 pointr = pointr % *m + 1;
04431 /* L20: */
04432 }
04433 /* Compute wv:=K^(-1)wv. */
04434 m2 = *m << 1;
04435 col2 = *col << 1;
04436 dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info);
04437 if (*info != 0) {
04438 return 0;
04439 }
04440 i__1 = *col;
04441 for (i__ = 1; i__ <= i__1; ++i__) {
04442 wv[i__] = -wv[i__];
04443 /* L25: */
04444 }
04445 dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info);
04446 if (*info != 0) {
04447 return 0;
04448 }
04449 /* Compute d = (1/theta)d + (1/theta**2)Z'W wv. */
04450 pointr = *head;
04451 i__1 = *col;
04452 for (jy = 1; jy <= i__1; ++jy) {
04453 js = *col + jy;
04454 i__2 = *nsub;
04455 for (i__ = 1; i__ <= i__2; ++i__) {
04456 k = ind[i__];
04457 d__[i__] = d__[i__] + wy[k + pointr * wy_dim1] * wv[jy] / *theta
04458 + ws[k + pointr * ws_dim1] * wv[js];
04459 /* L30: */
04460 }
04461 pointr = pointr % *m + 1;
04462 /* L40: */
04463 }
04464 i__1 = *nsub;
04465 for (i__ = 1; i__ <= i__1; ++i__) {
04466 d__[i__] /= *theta;
04467 /* L50: */
04468 }
04469 /* Backtrack to the feasible region. */
04470 alpha = 1.;
04471 temp1 = alpha;
04472 i__1 = *nsub;
04473 for (i__ = 1; i__ <= i__1; ++i__) {
04474 k = ind[i__];
04475 dk = d__[i__];
04476 if (nbd[k] != 0) {
04477 if (dk < 0. && nbd[k] <= 2) {
04478 temp2 = l[k] - x[k];
04479 if (temp2 >= 0.) {
04480 temp1 = 0.;
04481 } else if (dk * alpha < temp2) {
04482 temp1 = temp2 / dk;
04483 }
04484 } else if (dk > 0. && nbd[k] >= 2) {
04485 temp2 = u[k] - x[k];
04486 if (temp2 <= 0.) {
04487 temp1 = 0.;
04488 } else if (dk * alpha > temp2) {
04489 temp1 = temp2 / dk;
04490 }
04491 }
04492 if (temp1 < alpha) {
04493 alpha = temp1;
04494 ibd = i__;
04495 }
04496 }
04497 /* L60: */
04498 }
04499 if (alpha < 1.) {
04500 dk = d__[ibd];
04501 k = ind[ibd];
04502 if (dk > 0.) {
04503 x[k] = u[k];
04504 d__[ibd] = 0.;
04505 } else if (dk < 0.) {
04506 x[k] = l[k];
04507 d__[ibd] = 0.;
04508 }
04509 }
04510 i__1 = *nsub;
04511 for (i__ = 1; i__ <= i__1; ++i__) {
04512 k = ind[i__];
04513 x[k] += alpha * d__[i__];
04514 /* L70: */
04515 }
04516 if (*iprint >= 99) {
04517 /* if (alpha < 1.) {
04518 s_wsfe(&io___246);
04519 do_fio(&c__1, (char *)&alpha, (long int)sizeof(double));
04520 e_wsfe();
04521 } else {
04522 s_wsle(&io___247);
04523 do_lio(&c__9, &c__1, "SM solution inside the box", (long int)26);
04524 e_wsle();
04525 }
04526 if (*iprint > 100) {
04527 s_wsfe(&io___248);
04528 i__1 = *n;
04529 for (i__ = 1; i__ <= i__1; ++i__) {
04530 do_fio(&c__1, (char *)&x[i__], (long int)sizeof(double));
04531 }
04532 e_wsfe();
04533 }*/
04534 }
04535 if (alpha < 1.) {
04536 *iword = 1;
04537 } else {
04538 *iword = 0;
04539 }
04540 if (*iprint >= 99) {
04541 /* s_wsfe(&io___249);
04542 e_wsfe();*/
04543 }
04544 return 0;
04545 } /* subsm_ */
|
|
|
Definition at line 5172 of file lbfgsb.cpp. Referenced by mainlb_(). 05174 {
05175 // static float temp;
05176 // extern double etime_(float *);
05177 static float tarray[2];
05178
05179 /* ********* */
05180
05181 /* Subroutine timer */
05182
05183 /* This subroutine is used to determine user time. In a typical */
05184 /* application, the user time for a code segment requires calls */
05185 /* to subroutine timer to determine the initial and final time. */
05186
05187 /* The subroutine statement is */
05188
05189 /* subroutine timer(ttime) */
05190
05191 /* where */
05192
05193 /* ttime is an output variable which specifies the user time. */
05194
05195 /* Argonne National Laboratory and University of Minnesota. */
05196 /* MINPACK-2 Project. */
05197
05198 /* Modified October 1990 by Brett M. Averick. */
05199
05200 /* ********** */
05201 /* The first element of the array tarray specifies user time */
05202 //temp = etime_(tarray);
05203 tarray[0]=0.0;
05204 tarray[1]=1.0;
05205 *ttime = (double) tarray[0];
05206 return 0;
05207 } /* timer_ */
|
|
|
Definition at line 35 of file lbfgsb.cpp. Referenced by bmv_(), cauchy_(), dpofa_(), dtrsl_(), formk_(), lnsrlb_(), mainlb_(), matupd_(), and subsm_(). |
|
|
Definition at line 37 of file lbfgsb.cpp. |
|
|
Definition at line 39 of file lbfgsb.cpp. Referenced by lnsrlb_(). |
|
|
Definition at line 40 of file lbfgsb.cpp. Referenced by lnsrlb_(). |
|
|
Definition at line 41 of file lbfgsb.cpp. Referenced by lnsrlb_(). |
|
|
Definition at line 34 of file lbfgsb.cpp. |
|
|
Definition at line 44 of file lbfgsb.cpp. Referenced by EMAN::Util::multi_align_error(), EMAN::Util::twoD_fine_ali(), EMAN::Util::twoD_fine_ali_G(), and EMAN::Util::twoD_to_3D_ali(). |
1.3.9.1