ext/lbfgsb/src/lbfgsb.c in lbfgsb-0.2.0 vs ext/lbfgsb/src/lbfgsb.c in lbfgsb-0.3.0
- old
+ new
@@ -671,12 +671,12 @@
goto L111;
}
if (strncmp(task, "STOP", 4) == 0) {
if (strncmp(task + 6, "CPU", 3) == 0) {
/* restore the previous iterate. */
- dcopy_(n, &t[1], &c__1, &x[1], &c__1);
- dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &t[1], &c__1, &x[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
*f = fold;
}
goto L999;
}
}
@@ -707,11 +707,11 @@
}
iword = -1;
if (! cnstnd && col > 0) {
/* skip the search for GCP. */
- dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
wrk = updatd;
nseg = 0;
goto L333;
}
/**
@@ -832,12 +832,12 @@
&d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep,
&stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd,
csave, &isave[22], &dsave[17]);
if (info != 0 || iback >= 20) {
/* restore the previous iterate. */
- dcopy_(n, &t[1], &c__1, &x[1], &c__1);
- dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &t[1], &c__1, &x[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
*f = fold;
if (col == 0) {
/* abnormal termination. */
if (info == 0) {
info = -9;
@@ -908,17 +908,17 @@
/* Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
r__[i__] = g[i__] - r__[i__];
}
- rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
+ rr = lbfgsb_rb_ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
if (stp == 1.) {
dr = gd - gdold;
ddum = -gdold;
} else {
dr = (gd - gdold) * stp;
- dscal_(n, &stp, &d__[1], &c__1);
+ lbfgsb_rb_dscal_(n, &stp, &d__[1], &c__1);
ddum = -gdold * stp;
}
if (dr <= epsmch * ddum) {
/* skip the L-BFGS update. */
++nskip;
@@ -1191,11 +1191,11 @@
sum += sy[i__ + k * sy_dim1] * v[k] / sy[k + k * sy_dim1];
}
p[i2] = v[i2] + sum;
}
/* Solve the triangular system */
- dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info);
+ lbfgsb_rb_dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__11, info);
if (*info != 0) {
return 0;
}
/* solve D^(1/2)p1=v1. */
i__1 = *col;
@@ -1203,11 +1203,11 @@
p[i__] = v[i__] / sqrt(sy[i__ + i__ * sy_dim1]);
}
/* PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ] */
/* [ 0 J' ] [ p2 ] [ p2 ]. */
/* solve J^Tp2=p2. */
- dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info);
+ lbfgsb_rb_dtrsl_(&wt[wt_offset], m, col, &p[*col + 1], &c__1, info);
if (*info != 0) {
return 0;
}
/* compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2) */
/* =-D^(-1/2)p1+D^(-1)L'p2. */
@@ -1461,11 +1461,11 @@
/* the derivative f1 and the vector p = W'd (for theta = 1). */
if (*sbgnrm <= 0.) {
if (*iprint >= 0) {
fprintf(stdout, " Subgnorm = 0. GCP = X.\n");
}
- dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
return 0;
}
bnded = TRUE_;
nfree = *n + 1;
nbreak = 0;
@@ -1560,14 +1560,14 @@
/* The indices of the nonzero components of d are now stored */
/* in iorder(1),...,iorder(nbreak) and iorder(nfree),...,iorder(n). */
/* The smallest of the nbreak breakpoints is in t(ibkmin)=bkmin. */
if (*theta != 1.) {
/* complete the initialization of p for theta not= one. */
- dscal_(col, theta, &p[*col + 1], &c__1);
+ lbfgsb_rb_dscal_(col, theta, &p[*col + 1], &c__1);
}
/* Initialize GCP xcp = x. */
- dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &x[1], &c__1, &xcp[1], &c__1);
if (nbreak == 0 && nfree == *n + 1) {
/* is a zero vector, return with the initial xcp as GCP. */
if (*iprint > 100) {
fprintf(stdout, "Cauchy X = \n");
fprintf(stdout, " ");
@@ -1594,11 +1594,11 @@
if (*col > 0) {
bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &p[1], &v[1], info);
if (*info != 0) {
return 0;
}
- f2 -= ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
+ f2 -= lbfgsb_rb_ddot_(&col2, &v[1], &c__1, &p[1], &c__1);
}
dtm = -f1 / f2;
tsum = 0.;
*nseg = 1;
if (*iprint >= 99) {
@@ -1683,11 +1683,11 @@
/* temporarily set f1 and f2 for col=0. */
f1 = f1 + dt * f2 + dibp2 - *theta * dibp * zibp;
f2 -= *theta * dibp2;
if (*col > 0) {
/* update c = c + dt*p. */
- daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
+ lbfgsb_rb_daxpy_(&col2, &dt, &p[1], &c__1, &c__[1], &c__1);
/* choose wbp, */
/* the row of W corresponding to the breakpoint encountered. */
pointr = *head;
i__1 = *col;
for (j = 1; j <= i__1; ++j) {
@@ -1698,16 +1698,16 @@
/* compute (wbp)Mc, (wbp)Mp, and (wbp)M(wbp)'. */
bmv_(m, &sy[sy_offset], &wt[wt_offset], col, &wbp[1], &v[1], info);
if (*info != 0) {
return 0;
}
- wmc = ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
- wmp = ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
- wmw = ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
+ wmc = lbfgsb_rb_ddot_(&col2, &c__[1], &c__1, &v[1], &c__1);
+ wmp = lbfgsb_rb_ddot_(&col2, &p[1], &c__1, &v[1], &c__1);
+ wmw = lbfgsb_rb_ddot_(&col2, &wbp[1], &c__1, &v[1], &c__1);
/* update p = p - dibp*wbp. */
d__1 = -dibp;
- daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
+ lbfgsb_rb_daxpy_(&col2, &d__1, &wbp[1], &c__1, &p[1], &c__1);
/* complete updating f1 and f2 while col > 0. */
f1 += dibp * wmc;
f2 = f2 + dibp * 2. * wmp - dibp2 * wmw;
}
d__1 = *epsmch * f2_org__;
@@ -1735,16 +1735,16 @@
dtm = 0.;
}
tsum += dtm;
/* Move free variables (i.e., the ones w/o breakpoints) and */
/* the variables whose breakpoints haven't been reached. */
- daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
+ lbfgsb_rb_daxpy_(n, &tsum, &d__[1], &c__1, &xcp[1], &c__1);
L999:
/* Update c = c + dtm*p = W'(x^c - x) */
/* which will be used in computing r = Z'(B(x^c - x) + g). */
if (*col > 0) {
- daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
+ lbfgsb_rb_daxpy_(&col2, &dtm, &p[1], &c__1, &c__[1], &c__1);
}
if (*iprint > 100) {
fprintf(stdout, "Cauchy X = \n");
fprintf(stdout, " ");
i__1 = *n;
@@ -2058,15 +2058,15 @@
/* shift old part of WN1. */
i__1 = *m - 1;
for (jy = 1; jy <= i__1; ++jy) {
js = *m + jy;
i__2 = *m - jy;
- dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[jy + jy * wn1_dim1], &c__1);
+ lbfgsb_rb_dcopy_(&i__2, &wn1[jy + 1 + (jy + 1) * wn1_dim1], &c__1, &wn1[jy + jy * wn1_dim1], &c__1);
i__2 = *m - jy;
- dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[js + js * wn1_dim1], &c__1);
+ lbfgsb_rb_dcopy_(&i__2, &wn1[js + 1 + (js + 1) * wn1_dim1], &c__1, &wn1[js + js * wn1_dim1], &c__1);
i__2 = *m - 1;
- dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[*m + 1 + jy * wn1_dim1], &c__1);
+ lbfgsb_rb_dcopy_(&i__2, &wn1[*m + 2 + (jy + 1) * wn1_dim1], &c__1, &wn1[*m + 1 + jy * wn1_dim1], &c__1);
}
}
/* put new rows in blocks (1,1), (2,1) and (2,2). */
pbegin = 1;
pend = *nsub;
@@ -2213,32 +2213,32 @@
}
/* Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')] */
/* [(-L_a +R_z)L'^-1 S'AA'S*theta ] */
/* first Cholesky factor (1,1) block of wn to get LL' */
/* with L' stored in the upper triangle of wn. */
- dpofa_(&wn[wn_offset], &m2, col, info);
+ lbfgsb_rb_dpofa_(&wn[wn_offset], &m2, col, info);
if (*info != 0) {
*info = -1;
return 0;
}
/* then form L^-1(-L_a'+R_z') in the (1,2) block. */
col2 = *col << 1;
i__1 = col2;
for (js = *col + 1; js <= i__1; ++js) {
- dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info);
+ lbfgsb_rb_dtrsl_(&wn[wn_offset], &m2, col, &wn[js * wn_dim1 + 1], &c__11, info);
}
/* Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the */
/* upper triangle of (2,2) block of wn. */
i__1 = col2;
for (is = *col + 1; is <= i__1; ++is) {
i__2 = col2;
for (js = is; js <= i__2; ++js) {
- wn[is + js * wn_dim1] += ddot_(col, &wn[is * wn_dim1 + 1], &c__1, &wn[js * wn_dim1 + 1], &c__1);
+ wn[is + js * wn_dim1] += lbfgsb_rb_ddot_(col, &wn[is * wn_dim1 + 1], &c__1, &wn[js * wn_dim1 + 1], &c__1);
}
}
/* Cholesky factorization of (2,2) block of wn. */
- dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info);
+ lbfgsb_rb_dpofa_(&wn[*col + 1 + (*col + 1) * wn_dim1], &m2, col, info);
if (*info != 0) {
*info = -2;
return 0;
}
return 0;
@@ -2301,11 +2301,11 @@
wt[i__ + j * wt_dim1] = ddum + *theta * ss[i__ + j * ss_dim1];
}
}
/* Cholesky factorize T to J*J' with */
/* J' stored in the upper triangle of wt. */
- dpofa_(&wt[wt_offset], m, col, info);
+ lbfgsb_rb_dpofa_(&wt[wt_offset], m, col, info);
if (*info != 0) {
*info = -3;
}
return 0;
}
@@ -2561,11 +2561,11 @@
--dsave;
if (strncmp(task, "FG_LN", 5) == 0) {
goto L556;
}
- *dtd = ddot_(n, &d__[1], &c__1, &d__[1], &c__1);
+ *dtd = lbfgsb_rb_ddot_(n, &d__[1], &c__1, &d__[1], &c__1);
*dnorm = sqrt(*dtd);
/* Determine the maximum step length. */
*stpmx = 1e10;
if (*cnstnd) {
if (*iter == 0) {
@@ -2598,18 +2598,18 @@
d__1 = 1. / *dnorm;
*stp = d__1 <= *stpmx ? d__1 : *stpmx;
} else {
*stp = 1.;
}
- dcopy_(n, &x[1], &c__1, &t[1], &c__1);
- dcopy_(n, &g[1], &c__1, &r__[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &x[1], &c__1, &t[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &g[1], &c__1, &r__[1], &c__1);
*fold = *f;
*ifun = 0;
*iback = 0;
strcpy(csave, "START");
L556:
- *gd = ddot_(n, &g[1], &c__1, &d__[1], &c__1);
+ *gd = lbfgsb_rb_ddot_(n, &g[1], &c__1, &d__[1], &c__1);
if (*ifun == 0) {
*gdold = *gd;
if (*gd >= 0.) {
/* the directional derivative >=0. */
/* Line search is impossible. */
@@ -2624,11 +2624,11 @@
strcpy(task, "FG_LNSRCH");
++(*ifun);
++(*nfgv);
*iback = *ifun - 1;
if (*stp == 1.) {
- dcopy_(n, &z__[1], &c__1, &x[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &z__[1], &c__1, &x[1], &c__1);
} else {
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
x[i__] = *stp * d__[i__] + t[i__];
}
@@ -2691,33 +2691,33 @@
} else {
*itail = *itail % *m + 1;
*head = *head % *m + 1;
}
/* Update matrices WS and WY. */
- dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1);
- dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1);
+ lbfgsb_rb_dcopy_(n, &d__[1], &c__1, &ws[*itail * ws_dim1 + 1], &c__1);
+ lbfgsb_rb_dcopy_(n, &r__[1], &c__1, &wy[*itail * wy_dim1 + 1], &c__1);
/* Set theta=yy/ys. */
*theta = *rr / *dr;
/* Form the middle matrix in B. */
/* update the upper triangle of SS, */
/* and the lower triangle of SY: */
if (*iupdat > *m) {
/* move old information */
i__1 = *col - 1;
for (j = 1; j <= i__1; ++j) {
- dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1], &c__1);
+ lbfgsb_rb_dcopy_(&j, &ss[(j + 1) * ss_dim1 + 2], &c__1, &ss[j * ss_dim1 + 1], &c__1);
i__2 = *col - j;
- dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j * sy_dim1], &c__1);
+ lbfgsb_rb_dcopy_(&i__2, &sy[j + 1 + (j + 1) * sy_dim1], &c__1, &sy[j + j * sy_dim1], &c__1);
}
}
/* add new information: the last row of SY */
/* and the last column of SS: */
pointr = *head;
i__1 = *col - 1;
for (j = 1; j <= i__1; ++j) {
- sy[*col + j * sy_dim1] = ddot_(n, &d__[1], &c__1, &wy[pointr * wy_dim1 + 1], &c__1);
- ss[j + *col * ss_dim1] = ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, &d__[1], &c__1);
+ sy[*col + j * sy_dim1] = lbfgsb_rb_ddot_(n, &d__[1], &c__1, &wy[pointr * wy_dim1 + 1], &c__1);
+ ss[j + *col * ss_dim1] = lbfgsb_rb_ddot_(n, &ws[pointr * ws_dim1 + 1], &c__1, &d__[1], &c__1);
pointr = pointr % *m + 1;
}
if (*stp == 1.) {
ss[*col + *col * ss_dim1] = *dtd;
} else {
@@ -3359,19 +3359,19 @@
pointr = pointr % *m + 1;
}
/* Compute wv:=K^(-1)wv. */
m2 = *m << 1;
col2 = *col << 1;
- dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info);
+ lbfgsb_rb_dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__11, info);
if (*info != 0) {
return 0;
}
i__1 = *col;
for (i__ = 1; i__ <= i__1; ++i__) {
wv[i__] = -wv[i__];
}
- dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info);
+ lbfgsb_rb_dtrsl_(&wn[wn_offset], &m2, &col2, &wv[1], &c__1, info);
if (*info != 0) {
return 0;
}
/* Compute d = (1/theta)d + (1/theta**2)Z'W wv. */
pointr = *head;
@@ -3385,16 +3385,16 @@
+ ws[k + pointr * ws_dim1] * wv[js];
}
pointr = pointr % *m + 1;
}
d__1 = 1. / *theta;
- dscal_(nsub, &d__1, &d__[1], &c__1);
+ lbfgsb_rb_dscal_(nsub, &d__1, &d__[1], &c__1);
/* ----------------------------------------------------------------- */
/* Let us try the projection, d is the Newton direction */
*iword = 0;
- dcopy_(n, &x[1], &c__1, &xp[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &x[1], &c__1, &xp[1], &c__1);
i__1 = *nsub;
for (i__ = 1; i__ <= i__1; ++i__) {
k = ind[i__];
dk = d__[i__];
@@ -3443,10 +3443,10 @@
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
dd_p__ += (x[i__] - xx[i__]) * gg[i__];
}
if (dd_p__ > 0.) {
- dcopy_(n, &xp[1], &c__1, &x[1], &c__1);
+ lbfgsb_rb_dcopy_(n, &xp[1], &c__1, &x[1], &c__1);
fprintf(stderr, " Positive dir derivative in projection\n");
fprintf(stderr, " Using the backtracking step\n");
} else {
goto L911;
}