| 565 | |
| 566 | |
| 567 | void LBFGS::lbfgs(int n, int m, double f, Eigen::VectorXd g, double eps, |
| 568 | double xtol) |
| 569 | { |
| 570 | bool execute_entire_while_loop = false; |
| 571 | if(iflag == 0) { |
| 572 | iter=0; |
| 573 | if ( n <= 0 || m <= 0 ) |
| 574 | { |
| 575 | iflag= -3; |
| 576 | } |
| 577 | |
| 578 | nfun= 1; |
| 579 | point = 0; |
| 580 | finish = false; |
| 581 | ispt = n + 2*m; |
| 582 | iypt = ispt + n*m; |
| 583 | npt = 0; |
| 584 | w.segment(ispt, n) = (-g).cwiseProduct(diag); |
| 585 | |
| 586 | stp1 = 1.0 / g.norm(); |
| 587 | ftol= 0.0001; |
| 588 | maxfev= 20; |
| 589 | execute_entire_while_loop = true; |
| 590 | } |
| 591 | while(true) { |
| 592 | if(execute_entire_while_loop) { |
| 593 | iter++; |
| 594 | info = 0; |
| 595 | bound=iter-1; |
| 596 | if (iter!=1) { |
| 597 | if (iter > m) bound = m; |
| 598 | double ys = w.segment(iypt + npt, n).dot(w.segment(ispt + npt, n)); |
| 599 | double yy = w.segment(iypt + npt, n).squaredNorm(); |
| 600 | diag.setConstant(ys / yy); |
| 601 | cp = point; |
| 602 | if (point ==0 ) cp =m; |
| 603 | w[n + cp-1] = 1.0 / ys; |
| 604 | w.head(n) = -g; |
| 605 | cp = point; |
| 606 | for (int i = 0; i < bound; i++) { |
| 607 | cp -= 1; |
| 608 | if (cp == -1) { |
| 609 | cp = m - 1; |
| 610 | } |
| 611 | sq = w.segment(ispt + cp *n,n).dot(w.head(n)); |
| 612 | inmc = n + m + cp; |
| 613 | iycn = iypt + cp * n; |
| 614 | w[inmc] = sq * w[n + cp]; |
| 615 | w.head(n) -= w[inmc] * w.segment(iycn, n); |
| 616 | } |
| 617 | w.head(n)=w.head(n).cwiseProduct(diag); |
| 618 | |
| 619 | for (int i = 0; i < bound; i++) { |
| 620 | yr = w.segment(iypt + cp * n, n).dot(w.head(n)); |
| 621 | inmc = n + m + cp; |
| 622 | beta = w[inmc] - w[n + cp] * yr; |
| 623 | iscn = ispt + cp * n; |
| 624 | w.head(n) += beta * w.segment(iscn, n); |