This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Functions | |
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 *, long int *, double *, double *, double *, long int *, long int *, double *) |
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) |
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) |
|
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_ */
|