/* * R : A Computer Language for Statistical Data Analysis * Copyright (C) 1999-2023 The R Core Team * Copyright (C) 2003-2017 The R Foundation * Copyright (C) 1997-1999 Saikat DebRoy * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, a copy is available at * https://www.R-project.org/Licenses/ * USA */ /* ../appl/uncmin.f -- translated by f2c (version of 1 June 1993 23:00:00). -- and hand edited by Saikat DebRoy */ /*--- The Dennis + Schnabel Minimizer -- used by R's nlm() ---*/ #include #include /* DBL_MAX */ #include #include /* Rprintf */ #include /* printRealVector */ #include /* ddot dnrm2 dscal */ #include /* dtrsl */ #include /* fdhess */ #include // as in : #define Rexp10(x) pow(10.0, x) /* CC subroutines mvmlt[lsu] should be REPLACED by BLAS ones! * CC * CC--- choldc(nr,n,a,diagmx,tol,addmax) is ``cholesky + tolerance'' * CC ------ * CC it should make use of BLAS routines as [linkpack's dpofa!] */ void fdhess(int n, double *x, double fval, fcn_p fun, void *state, double *h, int nfd, double *step, double *f, int ndigit, double *typx) { /* calculates a numerical approximation to the upper triangular * portion of the second derivative matrix (the hessian). * Algorithm A5.6.2 from Dennis and Schnabel (1983), numerical methods * for unconstrained optimization and nonlinear equations, * prentice-hall, 321-322. * programmed by richard h. jones, january 11, 1989 * INPUT to subroutine * n the number of parameters * x vector of parameter values * fval double precision value of function at x * fun a function provided by the user which must be declared as * external in the calling program. its call must * be of the call fun(n,x,state,fval) where fval is the * computed value of the function * state information other than x and n that fun requires. * state is not modified in fdhess (but can be modified by fun). * nfd first dimension of h in the calling program * OUTPUT from subroutine * h an n by n matrix of the approximate hessian * Work space : * step a real array of length n * f a double precision array of length n */ int i, j; double tempi, tempj, fii, eta, fij; eta = Rexp10(-ndigit/3.0); for (i = 0; i < n; ++i) { step[i] = eta * fmax2(x[i], typx[i]); if (typx[i] < 0.) step[i] = -step[i]; tempi = x[i]; x[i] += step[i]; step[i] = x[i] - tempi; (*fun)(n, x, &f[i], state); x[i] = tempi; } for (i = 0; i < n; ++i) { tempi = x[i]; x[i] += step[i] * 2.; (*fun)(n, x, &fii, state); h[i + i * nfd] = (fval - f[i] + (fii - f[i]))/(step[i] * step[i]); x[i] = tempi + step[i]; for (j = i + 1; j < n; ++j) { tempj = x[j]; x[j] += step[j]; (*fun)(n, x, &fij, state); h[i + j * nfd] = (fval - f[i] + (fij - f[j]))/(step[i] * step[j]); x[j] = tempj; } x[i] = tempi; } } /* fdhess */ #if 0 static void d1fcn_dum(int n, double *x, double *g, void *state) { /* dummy routine to prevent unsatisfied external diagnostic * when specific analytic gradient function not supplied. */ } static void d2fcn_dum(int nr, int n, double *x, double *h, void *state) { /* dummy routine to prevent unsatisfied external diagnostic * when specific analytic hessian function not supplied. */ } #endif static void mvmltl(int nr, int n, double *a, double *x, double *y) { /* compute y = l x * where l is a lower triangular matrix stored in a * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) --> lower triangular (n*n) matrix * x(n) --> operand vector * y(n) <-- result vector * note * x and y cannot share storage */ int i, j; double sum; for (i = 0; i < n; ++i) { sum = 0.; for (j = 0; j <= i; ++j) sum += a[i + j * nr] * x[j]; y[i] = sum; } } /* mvmltl */ static void mvmltu(int nr, int n, double *a, double *x, double *y) { /* compute y = (L+) x * where L is a lower triangular matrix stored in a * (L-transpose (L+) is taken implicitly) * ARGUMENTS : * nr --> row dimension of matrix * n --> dimension of problem * a(nr,1) --> lower triangular (n*n) matrix * x(n) --> operand vector * y(n) <-- result vector * NOTE : x and y cannot share storage */ int i, length, one = 1; for (i = 0, length = n; i < n; --length, ++i) y[i] = F77_CALL(ddot)(&length, &a[i + i * nr], &one, &x[i], &one); } /* mvmltu */ static void mvmlts(int nr, int n, double *a, double *x, double *y) { /* compute y=ax * where "a" is a symmetric (n*n) matrix stored in its lower * triangular part and x,y are n-vectors * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) --> symmetric (n*n) matrix stored in * lower triangular part and diagonal * x(n) --> operand vector * y(n) <-- result vector * NOTE: x and y cannot share storage. */ int i, j; double sum; for (i = 0; i < n; ++i) { sum = 0.; for (j = 0; j <= i; ++j) sum += a[i + j * nr] * x[j]; for (j = i+1; j < n; ++j) sum += a[j + i * nr] * x[j]; y[i] = sum; } } /* mvmlts */ static void lltslv(int nr, int n, double *a, double *x, double *b) { /* solve ax=b where a has the form l(l-transpose) * but only the lower triangular part, l, is stored. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) --> matrix of form l(l-transpose). * on return a is unchanged. * x(n) <-- solution vector * b(n) --> right-hand side vector * note * if b is not required by calling program, then * b and x may share the same storage. */ int job = 0, info; if (x != b) Memcpy(x, b, n); F77_CALL(dtrsl)(a, &nr, &n, x, &job, &info); job = 10; F77_CALL(dtrsl)(a, &nr, &n, x, &job, &info); } /* lltslv */ static void choldc(int nr, int n, double *a, double diagmx, double tol, double *addmax) { /* Find the perturbed L(L-transpose) [written LL+] decomposition * of a+d, where d is a non-negative diagonal matrix added to a if * necessary to allow the cholesky decomposition to continue. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) <--> on entry: matrix for which to find perturbed * cholesky decomposition * on exit: contains L of LL+ decomposition * in lower triangular part and diagonal of "a" * diagmx --> maximum diagonal element of "a" * tol --> tolerance * addmax <-- maximum amount implicitly added to diagonal of "a" * in forming the cholesky decomposition of a+d * internal variables * aminl smallest element allowed on diagonal of L * amnlsq =aminl**2 * offmax maximum off-diagonal element in column of a * description * the normal cholesky decomposition is performed. however, if at any * point the algorithm would attempt to set L(i,i)=sqrt(temp) * with temp < tol*diagmx, then L(i,i) is set to sqrt(tol*diagmx) * instead. this is equivalent to adding tol*diagmx-temp to a(i,i) */ double tmp1, tmp2; int i, j, k; double aminl, offmax, amnlsq; double sum; *addmax = 0.0; aminl = sqrt(diagmx * tol); amnlsq = aminl * aminl; /* form row i of L */ for (i = 0; i < n; ++i) { // A[i,j] := * || find i,j element of lower triangular matrix L for (j = 0; j < i; ++j) { sum = 0.; for (k = 0; k < j; ++k) sum += a[i + k * nr] * a[j + k * nr]; a[i + j * nr] = (a[i + j * nr] - sum) / a[j + j * nr]; } // A[i,i] := * || find diagonal elements of L sum = 0.; for (k = 0; k < i; ++k) sum += a[i + k * nr] * a[i + k * nr]; tmp1 = a[i + i * nr] - sum; if (tmp1 >= amnlsq) { // normal Cholesky a[i + i * nr] = sqrt(tmp1); } else { // augment diagonal of L /* find maximum off-diagonal element in row */ offmax = 0.; for (j = 0; j < i; ++j) { if(offmax < (tmp2 = fabs(a[i + j * nr]))) offmax = tmp2; } if (offmax <= amnlsq) offmax = amnlsq; /* add to diagonal element to * allow cholesky decomposition to continue */ a[i + i * nr] = sqrt(offmax); if(*addmax < (tmp2 = offmax - tmp1)) *addmax = tmp2; } } } /* choldc */ static void qraux1(int nr, int n, double *r, int i) { /* Interchange rows i,i+1 of the upper hessenberg matrix r, columns i to n . * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of matrix * r[n*n] <--> upper hessenberg matrix * i --> index of row to interchange (i < n-1) */ double tmp; double *r1, *r2; /* pointer arithmetic : */ r1 = r + i + i * nr; r2 = r1 + 1; while(n-- > i) { tmp = *r1; *r1 = *r2; *r2 = tmp; r1 += nr; r2 += nr; } } /* qraux1 */ static void qraux2(int nr, int n, double *r, int i, double a, double b) { /* Pre-multiply r by the jacobi rotation j(i,i+1,a,b) . * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of matrix * r(n,n) <--> upper hessenberg matrix * i --> index of row * a --> scalar * b --> scalar */ double c, s; double y, z, den; double *r1, *r2; den = hypot(a,b); c = a / den; s = b / den; /* pointer arithmetic : */ r1 = r + i + i*nr; r2 = r1 + 1; while(n-- > i) { y = *r1; z = *r2; *r1 = c * y - s * z; *r2 = s * y + c * z; r1 += nr; r2 += nr; } } /* qraux2 */ static void qrupdt(int nr, int n, double *a, double *u, double *v) { /* Find an orthogonal (n*n) matrix (q*) and an upper triangular (n*n) * matrix (r*) such that (q*)(r*)=r+u(v+) * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) <--> on input: contains r * on output: contains (r*) * u(n) --> vector * v(n) --> vector */ int i, j, k; double t1, t2; int ii; /* determine last non-zero in u(.) */ for(k = n-1; k > 0 && u[k] == 0.0; k--) ; /* (k-1) jacobi rotations transform * r + u(v+) --> (r*) + (u(1)*e1)(v+) * which is upper hessenberg */ if (k > 0) { ii = k; while(ii > 0) { i = ii - 1; if (u[i] == 0.0) { qraux1(nr, n, a, i); u[i] = u[ii]; } else { qraux2(nr, n, a, i, u[i], -u[ii]); u[i] = hypot(u[i], u[ii]); } ii = i; } } /* r <-- r + (u(1)*e1)(v+) */ for (j = 0; j < n; ++j) a[j * nr] += u[0] * v[j]; /* (k-1) jacobi rotations transform upper hessenberg r * to upper triangular (r*) */ for (i = 0; i < k; ++i) { if (a[i + i * nr] == 0.) qraux1(nr, n, a, i); else { t1 = a[i + i * nr]; t2 = -a[i + 1 + i * nr]; qraux2(nr, n, a, i, t1, t2); } } } /* qrupdt */ static void tregup(int nr, int n, double *x, double f, double *g, double *a, fcn_p fcn, void *state, double *sc, double *sx, Rboolean nwtake, double stepmx, double steptl, double *dlt, int *iretcd, double *xplsp, double *fplsp, double *xpls, double *fpls, Rboolean *mxtake, int method, double *udiag) { /* TRust REGion UPdating * == == == * Decide whether to accept xpls = x+sc as the next iterate and * update the trust region radius dlt. * Used iff method == 2 or 3 * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> old iterate x[k-1] * f --> function value at old iterate, f(x) * g(n) --> gradient at old iterate, g(x), or approximate * a(n,n) --> cholesky decomposition of hessian in * lower triangular part and diagonal. * hessian or approx in upper triangular part * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in tregup (but can be * modified by fcn). * sc(n) --> current step * sx(n) --> diagonal scaling matrix for x * nwtake --> boolean, = TRUE if newton step taken * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * dlt <--> trust region radius * iretcd <--> return code * =0 xpls accepted as next iterate; * dlt trust region for next iteration. * =1 xpls unsatisfactory but accepted as next iterate * because xpls-x < smallest allowable step length. * =2 f(xpls) too large. continue current iteration * with new reduced dlt. * =3 f(xpls) sufficiently small, but quadratic model * predicts f(xpls) sufficiently well to continue * current iteration with new doubled dlt. * xplsp(n) <--> workspace [value needs to be retained between * succesive calls of k-th global step] * fplsp <--> [retain value between successive calls] * xpls(n) <-- new iterate x[k] * fpls <-- function value at new iterate, f(xpls) * mxtake <-- boolean flag indicating step of maximum length used * ipr --> device to which to send output * method --> algorithm to use to solve minimization problem * =1 line search * =2 double dogleg * =3 more-hebdon * udiag(n) --> diagonal of hessian in a(.,.) */ double dltf; double temp1; int i, j, one = 1; double dltfp, dltmp; double rln, slp; *mxtake = FALSE; for (i = 0; i < n; ++i) xpls[i] = x[i] + sc[i]; (*fcn)(n, xpls, fpls, state); dltf = *fpls - f; slp = F77_CALL(ddot)(&n, g, &one, sc, &one); /* next statement added for case of compilers which do not optimize evaluation of next "if" statement (in which case fplsp could be undefined). if (*iretcd == 4) { *fplsp = 0.; } */ if (*iretcd == 3 && (*fpls >= *fplsp || dltf > slp * 1e-4)) { /* reset xpls to xplsp and terminate global step */ *iretcd = 0; for (i = 0; i < n; ++i) xpls[i] = xplsp[i]; *fpls = *fplsp; *dlt *= .5; } else { /* fpls too large */ if (dltf > slp * 1e-4) { rln = 0.; for (i = 0; i < n; ++i) { temp1 = fabs(sc[i])/fmax2(fabs(xpls[i]), 1./sx[i]); if(rln < temp1) rln = temp1; } if (rln < steptl) { /* cannot find satisfactory xpls sufficiently distinct from x */ *iretcd = 1; } else { /* reduce trust region and continue global step */ *iretcd = 2; dltmp = -slp * *dlt / ((dltf - slp) * 2.); if (dltmp < *dlt * .1) *dlt *= .1; else *dlt = dltmp; } } else { /* fpls sufficiently small */ dltfp = 0.; if (method == 2) { for (i = 0; i < n; ++i) { temp1 = 0.; for (j = i; j < n; ++j) temp1 += a[j + i * nr] * sc[j]; dltfp += temp1 * temp1; } } else { /* method != 2 */ for (i = 0; i < n; ++i) { dltfp += udiag[i] * sc[i] * sc[i]; temp1 = 0.; for (j = i+1; j < n; ++j) temp1 += a[i + j * nr] * sc[i] * sc[j]; dltfp += temp1 * 2.; } } dltfp = slp + dltfp / 2.; if (*iretcd != 2 && fabs(dltfp - dltf) <= fabs(dltf) * 0.1 && nwtake && *dlt <= stepmx * .99) { /* double trust region and continue global step */ *iretcd = 3; for (i = 0; i < n; ++i) xplsp[i] = xpls[i]; *fplsp = *fpls; temp1 = *dlt * 2.0; *dlt = fmin2(temp1, stepmx); } else { /* accept xpls as next iterate. choose new trust region. */ *iretcd = 0; if (*dlt > stepmx * .99) *mxtake = TRUE; if (dltf >= dltfp * .1) { /* decrease trust region for next iteration */ *dlt *= .5; } else { /* check whether to increase trust region for next iteration */ if (dltf <= dltfp * .75) { temp1 = *dlt * 2.0; *dlt = fmin2(temp1, stepmx); } } } } } } /* tregup */ static void lnsrch(int n, double *x, double f, double *g, double *p, double *xpls, double *fpls, fcn_p fcn, void *state, Rboolean *mxtake, int *iretcd, double stepmx, double steptl, double *sx) { /* Find a next newton iterate by line search. (iff method == 1) * PARAMETERS : * n --> dimension of problem * x(n) --> old iterate: x[k-1] * f --> function value at old iterate, f(x) * g(n) --> gradient at old iterate, g(x), or approximate * p(n) --> non-zero newton step * xpls(n) <-- new iterate x[k] * fpls <-- function value at new iterate, f(xpls) * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in lnsrch (but can be * modified by fcn). * iretcd <-- return code * mxtake <-- boolean flag indicating step of maximum length used * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * sx(n) --> diagonal scaling matrix for x * internal variables * sln newton length * rln relative length of newton step */ int i, one = 1; Rboolean firstback = TRUE; double disc; double a3, b; double t1, t2, t3, lambda, tlmbda, rmnlmb; double scl, rln, sln, slp; double temp1; double pfpls = 0., plmbda = 0.; /* -Wall */ temp1 = 0.; for (i = 0; i < n; ++i) temp1 += sx[i] * sx[i] * p[i] * p[i]; sln = sqrt(temp1); if (sln > stepmx) { /* newton step longer than maximum allowed */ scl = stepmx / sln; F77_CALL(dscal)(&n, &scl, p, &one); sln = stepmx; } slp = F77_CALL(ddot)(&n, g, &one, p, &one); rln = 0.; for (i = 0; i < n; ++i) { temp1 = fabs(p[i])/ fmax2(fabs(x[i]), 1./sx[i]); if(rln < temp1) rln = temp1; } rmnlmb = steptl / rln; lambda = 1.0; /* check if new iterate satisfactory. generate new lambda if necessary. */ *mxtake = FALSE; *iretcd = 2; do { for (i = 0; i < n; ++i) xpls[i] = x[i] + lambda * p[i]; (*fcn)(n, xpls, fpls, state); if (*fpls <= f + slp * 1e-4 * lambda) { /* solution found */ *iretcd = 0; if (lambda == 1. && sln > stepmx * .99) *mxtake = TRUE; return; } /* else : solution not (yet) found */ /* First find a point with a finite value */ if (lambda < rmnlmb) { /* no satisfactory xpls found sufficiently distinct from x */ *iretcd = 1; return; } else { /* calculate new lambda */ /* modifications by BDR 2000/01/05 to cover non-finite values * ">=" instead of "==" : MM 2001/07/24 */ if (*fpls >= DBL_MAX) { lambda *= 0.1; firstback = TRUE; } else { if (firstback) { /* first backtrack: quadratic fit */ tlmbda = -lambda * slp / ((*fpls - f - slp) * 2.); firstback = FALSE; } else { /* all subsequent backtracks: cubic fit */ t1 = *fpls - f - lambda * slp; t2 = pfpls - f - plmbda * slp; t3 = 1. / (lambda - plmbda); a3 = 3. * t3 * (t1 / (lambda * lambda) - t2 / (plmbda * plmbda)); b = t3 * (t2 * lambda / (plmbda * plmbda) - t1 * plmbda / (lambda * lambda)); disc = b * b - a3 * slp; if (disc > b * b) /* only one positive critical point, must be minimum */ tlmbda = (-b + ((a3 < 0)? -sqrt(disc): sqrt(disc))) /a3; else /* both critical points positive, first is minimum */ tlmbda = (-b + ((a3 < 0)? sqrt(disc): -sqrt(disc))) /a3; if (tlmbda > lambda * .5) tlmbda = lambda * .5; } plmbda = lambda; pfpls = *fpls; if (tlmbda < lambda * .1) lambda *= .1; else lambda = tlmbda; } } } while(*iretcd > 1); } /* lnsrch */ static void dog_1step(int nr, int n, double *g, double *a, double *p, double *sx, double rnwtln, double *dlt, Rboolean *nwtake, Rboolean *fstdog, double *ssd, double *v, double *cln, double *eta, double *sc, double stepmx) { /* Find new step by double dogleg algorithm (iff method == 2); * repeatedly called by dogdrv() only. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * g(n) --> gradient at current iterate, g(x) * a(n,n) --> cholesky decomposition of hessian in * lower part and diagonal * p(n) --> newton step * sx(n) --> diagonal scaling matrix for x * rnwtln --> newton step length * dlt <--> trust region radius * nwtake <--> boolean, =.true. if newton step taken * fstdog <--> boolean, =.true. if on first leg of dogleg * ssd(n) <--> workspace [cauchy step to the minimum of the * quadratic model in the scaled steepest descent * direction] [retain value between successive calls] * v(n) <--> workspace [retain value between successive calls] * cln <--> cauchy length * [retain value between successive calls] * eta [retain value between successive calls] * sc(n) <-- current step * ipr --> device to which to send output * stepmx --> maximum allowable step size * internal variables * cln length of cauchy step */ int i, j, one = 1; double alam, bet, alpha, tmp, dot1, dot2; /* can we take newton step */ *nwtake = (rnwtln <= *dlt); if (*nwtake) { for (i = 0; i < n; ++i) sc[i] = p[i]; *dlt = rnwtln; return; } /* else *nwtake = FALSE : * newton step too long -- cauchy step is on double dogleg curve */ if (*fstdog) { /* calculate double dogleg curve (ssd) */ *fstdog = FALSE; alpha = 0.; for (i = 0; i < n; ++i) alpha += g[i] * g[i] / (sx[i] * sx[i]); bet = 0.; for (i = 0; i < n; ++i) { tmp = 0.; for (j = i; j < n; ++j) tmp += a[j + i * nr] * g[j] / (sx[j] * sx[j]); bet += tmp * tmp; } for (i = 0; i < n; ++i) ssd[i] = -(alpha / bet) * g[i] / sx[i]; *cln = alpha * sqrt(alpha) / bet; *eta = (.8 * alpha * alpha / (-bet * F77_CALL(ddot)(&n, g, &one, p, &one))) + .2; for (i = 0; i < n; ++i) v[i] = *eta * sx[i] * p[i] - ssd[i]; if (*dlt == -1.) *dlt = fmin2(*cln, stepmx); } if (*eta * rnwtln <= *dlt) { /* take partial step in newton direction */ for (i = 0; i < n; ++i) sc[i] = *dlt / rnwtln * p[i]; } else if (*cln >= *dlt) { /* take step in steepest descent direction */ for (i = 0; i < n; ++i) sc[i] = *dlt / *cln * ssd[i] / sx[i]; } else { /* calculate convex combination of ssd and eta*p which has scaled length dlt */ dot1 = F77_CALL(ddot)(&n, v, &one, ssd, &one); dot2 = F77_CALL(ddot)(&n, v, &one, v, &one); alam = (-dot1 + sqrt(dot1 * dot1 - dot2 * (*cln * *cln - *dlt * *dlt))) / dot2; for (i = 0; i < n; ++i) sc[i] = (ssd[i] + alam * v[i]) / sx[i]; } } /* dog_1step */ static void dogdrv(int nr, int n, double *x, double f, double *g, double *a, double *p, double *xpls, double *fpls, fcn_p fcn, void *state, double *sx, double stepmx, double steptl, double *dlt, int *iretcd, Rboolean *mxtake, double *sc, double *wrk1, double *wrk2, double *wrk3, int *itncnt) { /* Find a next newton iterate (xpls) by the double dogleg method * (iff method == 2 ). * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> old iterate x[k-1] * f --> function value at old iterate, f(x) * g(n) --> gradient at old iterate, g(x), or approximate * a(n,n) --> cholesky decomposition of hessian * in lower triangular part and diagonal * p(n) --> newton step * xpls(n) <-- new iterate x[k] * fpls <-- function value at new iterate, f(xpls) * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in dogdrv (but can be * modified by fcn). * sx(n) --> diagonal scaling matrix for x * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * dlt <--> trust region radius * [retain value between successive calls] * iretcd <-- return code * =0 satisfactory xpls found * =1 failed to find satisfactory xpls sufficiently * distinct from x * mxtake <-- boolean flag indicating step of maximum length used * sc(n) --> workspace [current step] * wrk1(n) --> workspace (and place holding argument to tregup) * wrk2(n) --> workspace * wrk3(n) --> workspace * ipr --> device to which to send output */ Rboolean fstdog, nwtake; int i; double fplsp = 0.0, rnwtln, eta = 0.0, cln = 0.0, tmp; /* -Wall */ tmp = 0.; for (i = 0; i < n; ++i) tmp += sx[i] * sx[i] * p[i] * p[i]; rnwtln = sqrt(tmp); *iretcd = 4; fstdog = TRUE; do { /* find new step by double dogleg algorithm */ dog_1step(nr, n, g, a, p, sx, rnwtln, dlt, &nwtake, &fstdog, wrk1, wrk2, &cln, &eta, sc, stepmx); /* check new point and update trust region */ tregup(nr, n, x, f, g, a, (fcn_p)fcn, state, sc, sx, nwtake, stepmx, steptl, dlt, iretcd, wrk3, &fplsp, xpls, fpls, mxtake, 2, wrk1); } while(*iretcd > 1); } /* dogdrv */ static void hook_1step(int nr, int n, double *g, double *a, double *udiag, double *p, double *sx, double rnwtln, double *dlt, double *amu, double dltp, double *phi, double *phip0, Rboolean *fstime, double *sc, Rboolean *nwtake, double *wrk0, double epsm) { /* Find new step by more-hebdon algorithm (iff method == 3); * repeatedly called by hookdrv() only. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * g(n) --> gradient at current iterate, g(x) * a(n,n) --> cholesky decomposition of hessian in * lower triangular part and diagonal. * hessian or approx in upper triangular part * udiag(n) --> diagonal of hessian in a(.,.) * p(n) --> newton step * sx(n) --> diagonal scaling matrix for n * rnwtln --> newton step length * dlt <--> trust region radius * amu <--> [retain value between successive calls] * dltp --> trust region radius at last exit from this routine * phi <--> [retain value between successive calls] * phip0 <--> [retain value between successive calls] * fstime <--> boolean. =.true. if first entry to this routine * during k-th iteration * sc(n) <-- current step * nwtake <-- boolean, =.true. if newton step taken * wrk0(n) --> workspace * epsm --> machine epsilon */ int one = 1, job = 0, info; int i, j; double phip; double amulo, amuup; double addmax, stepln; double temp1; const double hi = 1.5, alo = 0.75; /* hi and alo are constants used in this routine. */ /* change here if other values are to be substituted. */ /* shall we take newton step ? */ *nwtake = (rnwtln <= hi * *dlt); if (*nwtake) { /* take newton step */ for (i = 0; i < n; ++i) sc[i] = p[i]; *dlt = fmin2(*dlt, rnwtln); *amu = 0.; return; } /* else *nwtake = FALSE : newton step not taken */ if (*amu > 0.) *amu -= (*phi + dltp) * (dltp - *dlt + *phi) / (*dlt * *phip0); *phi = rnwtln - *dlt; if (*fstime) { for (i = 0; i < n; ++i) wrk0[i] = sx[i] * sx[i] * p[i]; /* solve l*y = (sx**2)*p */ F77_CALL(dtrsl)(a, &nr, &n, wrk0, &job, &info); /* Computing 2nd power */ temp1 = F77_CALL(dnrm2)(&n, wrk0, &one); *phip0 = -(temp1 * temp1) / rnwtln; *fstime = FALSE; } phip = *phip0; amulo = -(*phi) / phip; amuup = 0.; for (i = 0; i < n; ++i) amuup += g[i] * g[i] / (sx[i] * sx[i]); amuup = sqrt(amuup) / *dlt; while (1) { /* test value of amu; generate next amu if necessary */ if (*amu < amulo || *amu > amuup) { *amu = fmax2(sqrt(amulo * amuup), amuup * .001); } /* copy (h,udiag) to l */ /* where h <-- h+amu*(sx**2) [do not actually change (h,udiag)] */ /* The original code was http://people.sc.fsu.edu/~jburkardt/f77_src/uncmin/uncmin.f do 100 j=1,n a(j,j)=udiag(j) + amu*sx(j)*sx(j) if(j.eq.n) go to 100 jp1=j+1 do i=jp1,n a(i,j)=a(j,i) end do 100 continue */ for (i = 0; i < n; ++i) { a[i + i * nr] = udiag[i] + *amu * sx[i] * sx[i]; for (j = 0; j < i; ++j) // changed from ++i 2012-05-31, PR# a[i + j * nr] = a[j + i * nr]; } /* factor h=l(l+) */ temp1 = sqrt(epsm); choldc(nr, n, a, 0.0, temp1, &addmax); /* solve h*p = l(l+)*sc = -g */ for (i = 0; i < n; ++i) wrk0[i] = -g[i]; lltslv(nr, n, a, sc, wrk0); /* reset h. note since udiag has not been destroyed we need do */ /* nothing here. h is in the upper part and in udiag, still intact */ stepln = 0.; for (i = 0; i < n; ++i) stepln += sx[i] * sx[i] * sc[i] * sc[i]; stepln = sqrt(stepln); *phi = stepln - *dlt; for (i = 0; i < n; ++i) wrk0[i] = sx[i] * sx[i] * sc[i]; F77_CALL(dtrsl)(a, &nr, &n, wrk0, &job, &info); temp1 = F77_CALL(dnrm2)(&n, wrk0, &one); phip = -(temp1 * temp1) / stepln; if ((alo * *dlt <= stepln && stepln <= hi * *dlt) || (amuup - amulo > 0.)) { /* sc is acceptable hookstep */ break; } else { /* sc not acceptable hookstep. select new amu */ temp1 = (*amu - *phi)/phip; amulo = fmax2(amulo, temp1); if (*phi < 0.) amuup = fmin2(amuup,*amu); *amu -= stepln * *phi / (*dlt * phip); } } } /* hook_1step */ static void hookdrv(int nr, int n, double *x, double f, double *g, double *a, double *udiag, double *p, double *xpls, double *fpls, fcn_p fcn, void *state, double *sx, double stepmx, double steptl, double *dlt, int *iretcd, Rboolean *mxtake, double *amu, double *dltp, double *phi, double *phip0, double *sc, double *xplsp, double *wrk0, double epsm, int itncnt) { /* Find a next newton iterate (xpls) by the more-hebdon method. * (iff method == 3) * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> old iterate x[k-1] * f --> function value at old iterate, f(x) * g(n) --> gradient at old iterate, g(x), or approximate * a(n,n) --> cholesky decomposition of hessian in lower * triangular part and diagonal. * hessian in upper triangular part and udiag. * udiag(n) --> diagonal of hessian in a(.,.) * p(n) --> newton step * xpls(n) <-- new iterate x[k] * fpls <-- function value at new iterate, f(xpls) * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in hookdrv (but can be * modified by fcn). * sx(n) --> diagonal scaling matrix for x * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * dlt <--> trust region radius * iretcd <-- return code * =0 satisfactory xpls found * =1 failed to find satisfactory xpls sufficiently * distinct from x * mxtake <-- boolean flag indicating step of maximum length used * amu <--> [retain value between successive calls] * dltp <--> [retain value between successive calls] * phi <--> [retain value between successive calls] * phip0 <--> [retain value between successive calls] * sc(n) --> workspace * xplsp(n) --> workspace * wrk0(n) --> workspace * epsm --> machine epsilon * itncnt --> iteration count * ipr --> device to which to send output */ Rboolean fstime, nwtake; int i, j; double bet, alpha, fplsp = 0.0 /* -Wall */, rnwtln, tmp; tmp = 0.; for (i = 0; i < n; ++i) tmp += sx[i] * sx[i] * p[i] * p[i]; rnwtln = sqrt(tmp); if (itncnt == 1) { *amu = 0.; /* if first iteration and trust region not provided by user, compute initial trust region. */ if (*dlt == -1.) { alpha = 0.; for (i = 0; i < n; ++i) alpha += g[i] * g[i] / (sx[i] * sx[i]); bet = 0.; for (i = 0; i < n; ++i) { tmp = 0.; for (j = i; j < n; ++j) tmp += a[j + i * nr] * g[j] / (sx[j] * sx[j]); bet += tmp * tmp; } *dlt = alpha * sqrt(alpha) / bet; if(*dlt > stepmx) *dlt = stepmx; } } *iretcd = 4; fstime = TRUE; do { /* find new step by more-hebdon algorithm */ hook_1step(nr, n, g, a, udiag, p, sx, rnwtln, dlt, amu, *dltp, phi, phip0, &fstime, sc, &nwtake, wrk0, epsm); *dltp = *dlt; /* check new point and update trust region */ tregup(nr, n, x, f, g, a, (fcn_p)fcn, state, sc, sx, nwtake, stepmx, steptl, dlt, iretcd, xplsp, &fplsp, xpls, fpls, mxtake, 3, udiag); } while(*iretcd > 1); } /* hookdrv */ static void secunf(int nr, int n, double *x, double *g, double *a, double *udiag, double *xpls, double *gpls, double epsm, int itncnt, double rnf, int iagflg, Rboolean *noupdt, double *s, double *y, double *t) { /* Update hessian by the bfgs unfactored method (only when method == 3) * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> old iterate, x[k-1] * g(n) --> gradient or approximate at old iterate * a(n,n) <--> on entry: approximate hessian at old iterate * in upper triangular part (and udiag) * on exit: updated approx hessian at new iterate * in lower triangular part and diagonal * [lower triangular part of symmetric matrix] * udiag --> on entry: diagonal of hessian * xpls(n) --> new iterate, x[k] * gpls(n) --> gradient or approximate at new iterate * epsm --> machine epsilon * itncnt --> iteration count * rnf --> relative noise in optimization function fcn * iagflg --> =1 if analytic gradient supplied, =0 otherwise * noupdt <--> boolean: no update yet * [retain value between successive calls] * s(n) --> workspace * y(n) --> workspace * t(n) --> workspace */ double ynrm2, snorm2; int i, j, one = 1; Rboolean skpupd; double gam, tol, den1, den2; /* copy hessian in upper triangular part and udiag to lower triangular part and diagonal */ for (i = 0; i < n; ++i) { a[i + i * nr] = udiag[i]; for (j = 0; j < i; ++j) a[i + j * nr] = a[j + i * nr]; } *noupdt = (itncnt == 1); for (i = 0; i < n; ++i) { s[i] = xpls[i] - x[i]; y[i] = gpls[i] - g[i]; } den1 = F77_CALL(ddot)(&n, s, &one, y, &one); snorm2 = F77_CALL(dnrm2)(&n, s, &one); ynrm2 = F77_CALL(dnrm2)(&n, y, &one); if (den1 < sqrt(epsm) * snorm2 * ynrm2) return; mvmlts(nr, n, a, s, t); den2 = F77_CALL(ddot)(&n, s, &one, t, &one); if (*noupdt) { /* h <-- [(s+)y/(s+)hs]h */ gam = den1 / den2; den2 *= gam; for (j = 0; j < n; ++j) { t[j] *= gam; for (i = j; i < n; ++i) a[i + j * nr] *= gam; } *noupdt = FALSE; } skpupd = TRUE; /* check update condition on row i */ for (i = 0; i < n; ++i) { tol = rnf * fmax2(fabs(g[i]), fabs(gpls[i])); if (iagflg == 0) tol /= sqrt(rnf); if (fabs(y[i] - t[i]) >= tol) { skpupd = FALSE; break; } } if (skpupd) return; /* bfgs update */ for (j = 0; j < n; ++j) { for (i = j; i < n; ++i) a[i + j * nr] += y[i] * y[j] / den1 - t[i] * t[j] / den2; } } /* secunf */ static void secfac(int nr, int n, double *x, double *g, double *a, double *xpls, double *gpls, double epsm, int itncnt, double rnf, int iagflg, Rboolean *noupdt, double *s, double *y, double *u, double *w) { /* Update hessian by the bfgs factored method (only when method == 1 or 2) * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> old iterate, x[k-1] * g(n) --> gradient or approximate at old iterate * a(n,n) <--> on entry: cholesky decomposition of hessian in * lower part and diagonal. * on exit: updated cholesky decomposition of hessian * in lower triangular part and diagonal * xpls(n) --> new iterate, x[k] * gpls(n) --> gradient or approximate at new iterate * epsm --> machine epsilon * itncnt --> iteration count * rnf --> relative noise in optimization function fcn * iagflg --> =1 if analytic gradient supplied, =0 itherwise * noupdt <--> boolean: no update yet * [retain value between successive calls] * s(n) --> workspace * y(n) --> workspace * u(n) --> workspace * w(n) --> workspace */ double ynrm2; int i, j, one = 1; Rboolean skpupd; double snorm2, reltol; double alp, den1, den2; *noupdt = (itncnt == 1); for (i = 0; i < n; ++i) { s[i] = xpls[i] - x[i]; y[i] = gpls[i] - g[i]; } den1 = F77_CALL(ddot)(&n, s, &one, y, &one); snorm2 = F77_CALL(dnrm2)(&n, s, &one); ynrm2 = F77_CALL(dnrm2)(&n, y, &one); if (den1 < sqrt(epsm) * snorm2 * ynrm2) return; mvmltu(nr, n, a, s, u); den2 = F77_CALL(ddot)(&n, u, &one, u, &one); /* l <-- sqrt(den1/den2)*l */ alp = sqrt(den1 / den2); if (*noupdt) { for (j = 0; j < n; ++j) { u[j] = alp * u[j]; for (i = j; i < n; ++i) { a[i + j * nr] *= alp; } } *noupdt = FALSE; den2 = den1; alp = 1.; } /* w = l(l+)s = hs */ mvmltl(nr, n, a, u, w); if (iagflg == 0) reltol = sqrt(rnf); else reltol = rnf; skpupd = TRUE; for (i = 0; i < n; ++i) { skpupd = (fabs(y[i] - w[i]) < reltol * fmax2(fabs(g[i]), fabs(gpls[i]))); if(!skpupd) break; } if(skpupd) return; /* w = y-alp*l(l+)s */ for (i = 0; i < n; ++i) w[i] = y[i] - alp * w[i]; /* alp=1/sqrt(den1*den2) */ alp /= den1; /* u=(l+)/sqrt(den1*den2) = (l+)s/sqrt((y+)s * (s+)l(l+)s) */ for (i = 0; i < n; ++i) u[i] *= alp; /* copy l into upper triangular part. zero l. */ for (i = 1; i < n; ++i) { for (j = 0; j < i; ++j) { a[j + i * nr] = a[i + j * nr]; a[i + j * nr] = 0.; } } /* find q, (l+) such that q(l+) = (l+) + u(w+) */ qrupdt(nr, n, a, u, w); /* upper triangular part and diagonal of a[] now contain updated * cholesky decomposition of hessian. * copy back to lower triangular part. */ for (i = 1; i < n; ++i) for (j = 0; j < i; ++j) a[i + j * nr] = a[j + i * nr]; } /* secfac */ static void chlhsn(int nr, int n, double *a, double epsm, double *sx, double *udiag) { /* find the l(l-transpose) [written LL+] decomposition of the perturbed * model hessian matrix a+mu*i(where mu\0 and i is the identity matrix) * which is safely positive definite. if a is safely positive definite * upon entry, then mu=0. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) <--> on entry; "a" is model hessian (only lower * triangular part and diagonal stored) * on exit: a contains l of LL+ decomposition of * perturbed model hessian in lower triangular * part and diagonal and contains hessian in upper * triangular part and udiag * epsm --> machine epsilon * sx(n) --> diagonal scaling matrix for x * udiag(n) <-- on exit: contains diagonal of hessian * INTERNAL VARIABLES * tol tolerance * diagmn minimum element on diagonal of a * diagmx maximum element on diagonal of a * offmax maximum off-diagonal element of a * offrow sum of off-diagonal elements in a row of a * evmin minimum eigenvalue of a * evmax maximum eigenvalue of a * DESCRIPTION * 1. if "a" has any negative diagonal elements, then choose mu>0 * such that the diagonal of a:=a+mu*i is all positive * with the ratio of its smallest to largest element on the * order of sqrt(epsm). * 2. "a" undergoes a perturbed cholesky decomposition which * results in an LL+ decomposition of a+d, where d is a * non-negative diagonal matrix which is implicitly added to * "a" during the decomposition if "a" is not positive definite. * "a" is retained and not changed during this process by * copying l into the upper triangular part of "a" and the * diagonal into udiag. then the cholesky decomposition routine * is called. on return, addmax contains maximum element of d. * 3. if addmax=0, "a" was positive definite going into step 2 * and return is made to calling program. otherwise, * the minimum number sdd which must be added to the * diagonal of a to make it safely strictly diagonally dominant * is calculated. since a+addmax*i and a+sdd*i are safely * positive definite, choose mu=fmin2(addmax,sdd) and decompose * a+mu*i to obtain l. */ int i, j; double evmin, evmax; double addmax, diagmn, diagmx, offmax, offrow, posmax; double sdd, amu, tol, tmp; /* scale hessian */ /* pre- and post- multiply "a" by inv(sx) */ for (j = 0; j < n; ++j) for (i = j; i < n; ++i) a[i + j * nr] /= sx[i] * sx[j]; /* step1 * ----- * note: if a different tolerance is desired throughout this * algorithm, change tolerance here: */ tol = sqrt(epsm); diagmx = a[0]; diagmn = a[0]; if (n > 1) { for (i = 1; i < n; ++i) { tmp = a[i + i * nr]; if(diagmn > tmp) diagmn = tmp; if(diagmx < tmp) diagmx = tmp; } } posmax = fmax2(diagmx, 0.0); if (diagmn <= posmax * tol) { amu = tol * (posmax - diagmn) - diagmn; if (amu == 0.) { /* find largest off-diagonal element of a */ offmax = 0.; for (i = 1; i < n; ++i) { for (j = 0; j < i; ++j) if (offmax < (tmp = fabs(a[i + j * nr]))) offmax = tmp; } if (offmax == 0.) amu = 1.; else amu = offmax * (tol + 1.); } /* a=a + mu*i */ for (i = 0; i < n; ++i) a[i + i * nr] += amu; diagmx += amu; } /* copy lower triangular part of "a" to upper triangular part */ /* and diagonal of "a" to udiag */ for (i = 0; i < n; ++i) { udiag[i] = a[i + i * nr]; for (j = 0; j < i; ++j) a[j + i * nr] = a[i + j * nr]; } choldc(nr, n, a, diagmx, tol, &addmax); /* step3 * if addmax=0, "a" was positive definite going into step 2, * the LL+ decomposition has been done, and we return. * otherwise, addmax>0. perturb "a" so that it is safely * diagonally dominant and find LL+ decomposition */ if (addmax > 0.0) { /* restore original "a" (lower triangular part and diagonal) */ for (i = 0; i < n; ++i) { a[i + i * nr] = udiag[i]; for (j = 0; j < i; ++j) a[i + j * nr] = a[j + i * nr]; } /* find sdd such that a+sdd*i is safely positive definite */ /* note: evmin<0 since a is not positive definite; */ evmin = 0.; evmax = a[0]; for (i = 0; i < n; ++i) { offrow = 0.; for (j = 0; j < i; ++j) offrow += fabs(a[i + j * nr]); for (j = i+1; j < n; ++j) offrow += fabs(a[j + i * nr]); tmp = a[i + i * nr] - offrow; if(evmin > tmp) evmin = tmp; tmp = a[i + i * nr] + offrow; if(evmax < tmp) evmax = tmp; } sdd = tol * (evmax - evmin) - evmin; /* perturb "a" and decompose again */ amu = fmin2(sdd, addmax); for (i = 0; i < n; ++i) { a[i + i * nr] += amu; udiag[i] = a[i + i * nr]; } /* "a" now guaranteed safely positive definite */ choldc(nr, n, a, 0.0, tol, &addmax); } /* unscale hessian and cholesky decomposition matrix */ for (j = 0; j < n; ++j) { for (i = j; i < n; ++i) a[i + j * nr] *= sx[i]; for (i = 0; i < j; ++i) a[i + j * nr] *= sx[i] * sx[j]; udiag[j] *= sx[j] * sx[j]; } } /* chlhsn */ static void hsnint(int nr, int n, double *a, double *sx, int method) { /* Provide initial hessian when using secant updates . * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * a(n,n) <-- initial hessian (lower triangular matrix) * sx(n) --> diagonal scaling matrix for x * method --> algorithm to use to solve minimization problem * =1,2 factored secant method used * =3 unfactored secant method used */ int i, j; for (i = 0; i < n; ++i) { if (method == 3) a[i + i * nr] = sx[i] * sx[i]; else a[i + i * nr] = sx[i]; for (j = 0; j < i; ++j) a[i + j * nr] = 0.; } } /* hsnint */ static void fstofd(int nr, int m, int n, double *xpls, fcn_p fcn, void *state, const double *fpls, double *a, double *sx, double rnoise, double *fhat, int icase) { /* find first order forward finite difference approximation "a" to the * first derivative of the function defined by the subprogram "fname" * evaluated at the new iterate "xpls". * for optimization use this routine to estimate: * 1) the first derivative (gradient) of the optimization function "fcn * analytic user routine has been supplied; * 2) the second derivative (hessian) of the optimization function * if no analytic user routine has been supplied for the hessian but * one has been supplied for the gradient ("fcn") and if the * optimization function is inexpensive to evaluate * note * _m=1 (optimization) algorithm estimates the gradient of the function * (fcn). fcn(x) # f: r(n)-->r(1) * _m=n (systems) algorithm estimates the jacobian of the function * fcn(x) # f: r(n)-->r(n). * _m=n (optimization) algorithm estimates the hessian of the optimizatio * function, where the hessian is the first derivative of "fcn" * PARAMETERS : * nr --> row dimension of matrix * m --> number of rows in a * n --> number of columns in a; dimension of problem * xpls(n) --> new iterate: x[k] * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in fstofd (but can be * modified by fcn). * fpls(m) --> _m=1 (optimization) function value at new iterate: * fcn(xpls) * _m=n (optimization) value of first derivative * (gradient) given by user function fcn * _m=n (systems) function value of associated * minimization function * a(nr,n) <-- finite difference approximation (see note). only * lower triangular matrix and diagonal are returned * sx(n) --> diagonal scaling matrix for x * rnoise --> relative noise in fcn [f(x)] * fhat(m) --> workspace * icase --> =1 optimization (gradient) * =2 systems * =3 optimization (hessian) * internal variables * stepsz - stepsize in the j-th variable direction */ int i, j; double xtmpj, stepsz, temp1, temp2; /* find j-th column of a each column is derivative of f(fcn) with respect to xpls(j) */ for (j = 0; j < n; ++j) { temp1 = fabs(xpls[j]); temp2 = 1.0/sx[j]; stepsz = sqrt(rnoise) * fmax2(temp1, temp2); xtmpj = xpls[j]; xpls[j] = xtmpj + stepsz; (*fcn)(n, xpls, fhat, state); xpls[j] = xtmpj; for (i = 0; i < m; ++i) a[i + j * nr] = (fhat[i] - fpls[i]) / stepsz; } if (icase == 3 && n > 1) { /* if computing hessian, a must be symmetric */ for (i = 1; i < m; ++i) for (j = 0; j < i; ++j) a[i + j * nr] = (a[i + j * nr] + a[j + i * nr]) / 2.0; } } /* fstofd */ static void fstocd(int n, double *x, fcn_p fcn, void *state, double *sx, double rnoise, double *g) { /* Find central difference approximation g to the first derivative * (gradient) of the function defined by fcn at the point x. * PARAMETERS : * n --> dimension of problem * x --> point at which gradient is to be approximated. * fcn --> name of subroutine to evaluate function. * state <--> information other than x and n that fcn requires. * state is not modified in fstocd (but can be * modified by fcn). * sx --> diagonal scaling matrix for x. * rnoise --> relative noise in fcn [f(x)]. * g <-- central difference approximation to gradient. */ int i; double stepi, fplus, fminus, xtempi, temp1, temp2; /* find i th stepsize, evaluate two neighbors in direction of i th */ /* unit vector, and evaluate i th component of gradient. */ for (i = 0; i < n; ++i) { xtempi = x[i]; temp1 = fabs(xtempi); temp2 = 1.0/sx[i]; stepi = pow(rnoise, 1.0/3.0) * fmax2(temp1, temp2); x[i] = xtempi + stepi; (*fcn)(n, x, &fplus, state); x[i] = xtempi - stepi; (*fcn)(n, x, &fminus, state); x[i] = xtempi; g[i] = (fplus - fminus) / (stepi * 2.); } } /* fstocd */ static void sndofd(int nr, int n, double *xpls, fcn_p fcn, void *state, double fpls, double *a, double *sx, double rnoise, double *stepsz, double *anbr) { /* Find second order forward finite difference approximation "a" * to the second derivative (hessian) of the function defined by the subp * "fcn" evaluated at the new iterate "xpls" . * For optimization use this routine to estimate * 1) the second derivative (hessian) of the optimization function * if no analytical user function has been supplied for either * the gradient or the hessian and if the optimization function * "fcn" is inexpensive to evaluate. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * xpls(n) --> new iterate: x[k] * fcn --> name of subroutine to evaluate function * state <--> information other than x and n that fcn requires. * state is not modified in sndofd (but can be * modified by fcn). * fpls --> function value at new iterate, f(xpls) * a(n,n) <-- finite difference approximation to hessian * only lower triangular matrix and diagonal * are returned * sx(n) --> diagonal scaling matrix for x * rnoise --> relative noise in fname [f(x)] * stepsz(n) --> workspace (stepsize in i-th component direction) * anbr(n) --> workspace (neighbor in i-th direction) */ double fhat; int i, j; double xtmpi, xtmpj; /* find i-th stepsize and evaluate neighbor in direction of i-th unit vector. */ for (i = 0; i < n; ++i) { xtmpi = xpls[i]; stepsz[i] = pow(rnoise, 1.0/3.0) * fmax2(fabs(xtmpi), 1./sx[i]); xpls[i] = xtmpi + stepsz[i]; (*fcn)(n, xpls, &anbr[i], state); xpls[i] = xtmpi; } /* calculate row i of a */ for (i = 0; i < n; ++i) { xtmpi = xpls[i]; xpls[i] = xtmpi + stepsz[i] * 2.; (*fcn)(n, xpls, &fhat, state); a[i + i * nr] = ((fpls - anbr[i]) + (fhat - anbr[i])) / (stepsz[i] * stepsz[i]); /* calculate sub-diagonal elements of column */ if(i == 0) { xpls[i] = xtmpi; continue; } xpls[i] = xtmpi + stepsz[i]; for (j = 0; j < i; ++j) { xtmpj = xpls[j]; xpls[j] = xtmpj + stepsz[j]; (*fcn)(n, xpls, &fhat, state); a[i + j*nr] = ((fpls - anbr[i]) + (fhat - anbr[j])) / (stepsz[i]*stepsz[j]); xpls[j] = xtmpj; } xpls[i] = xtmpi; } } /* sndofd */ static void grdchk(int n, double *x, fcn_p fcn, void *state, double f, double *g, double *typsiz, double *sx, double fscale, double rnf, double analtl, double *wrk1, int *msg) { /* Check analytic gradient against estimated gradient * PARAMETERS : * n --> dimension of problem * x(n) --> estimate to a root of fcn * fcn --> name of subroutine to evaluate optimization function * must be declared external in calling routine * fcn: r(n) --> r(1) * state <--> information other than x and n that fcn requires. * state is not modified in grdchk (but can be * modified by fcn). * f --> function value: fcn(x) * g(n) --> gradient: g(x) * typsiz(n) --> typical size for each component of x * sx(n) --> diagonal scaling matrix: sx(i)=1./typsiz(i) * fscale --> estimate of scale of objective function fcn * rnf --> relative noise in optimization function fcn * analtl --> tolerance for comparison of estimated and * analytical gradients * wrk1(n) --> workspace * msg <-- message or error code * on output: =-21, probable coding error of gradient */ int i; double gs, wrk; /* compute first order finite difference gradient and compare to analytic gradient. */ fstofd(1, 1, n, x, (fcn_p)fcn, state, &f, wrk1, sx, rnf, &wrk, 1); for (i = 0; i < n; ++i) { gs = fmax2(fabs(f), fscale) / fmax2(fabs(x[i]), typsiz[i]); if (fabs(g[i] - wrk1[i]) > fmax2(fabs(g[i]), gs) * analtl) { *msg = -21; return; } } } /* grdchk */ static void heschk(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn, void *state, double f, double *g, double *a, double *typsiz, double *sx, double rnf, double analtl, int iagflg, double *udiag, double *wrk1, double *wrk2, int *msg) { /* Check analytic hessian against estimated hessian * (this may be done only if the user supplied analytic hessian * d2fcn fills only the lower triangular part and diagonal of a) * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> estimate to a root of fcn * fcn --> name of subroutine to evaluate optimization function * must be declared external in calling routine * fcn: r(n) --> r(1) * d1fcn --> name of subroutine to evaluate gradient of fcn. * must be declared external in calling routine * d2fcn --> name of subroutine to evaluate hessian of fcn. * must be declared external in calling routine * state <--> information other than x and n that fcn, * d1fcn and d2fcn requires. * state is not modified in heschk (but can be * modified by fcn, d1fcn or d2fcn). * f --> function value: fcn(x) * g(n) <-- gradient: g(x) * a(n,n) <-- on exit: hessian in lower triangular part and diag * typsiz(n) --> typical size for each component of x * sx(n) --> diagonal scaling matrix: sx(i)=1./typsiz(i) * rnf --> relative noise in optimization function fcn * analtl --> tolerance for comparison of estimated and * analytical gradients * iagflg --> =1 if analytic gradient supplied * udiag(n) --> workspace * wrk1(n) --> workspace * wrk2(n) --> workspace * msg <--> message or error code * on input : if = 1xx do not compare anal + est hess * on output: = -22, probable coding error of hessian */ int i, j; double hs, temp1, temp2; /* compute finite difference approximation a to the hessian. */ if (iagflg) fstofd(nr, n, n, x, (fcn_p)d1fcn, state, g, a, sx, rnf, wrk1, 3); else sndofd(nr, n, x, (fcn_p)fcn, state, f, a, sx, rnf, wrk1, wrk2); /* copy lower triangular part of "a" to upper triangular part and diagonal of "a" to udiag */ for (j = 0; j < n; ++j) { udiag[j] = a[j + j * nr]; for (i = j+1; i < n; ++i) a[j + i * nr] = a[i + j * nr]; } /* compute analytic hessian and compare to finite difference approximation. */ (*d2fcn)(nr, n, x, a, state); for (j = 0; j < n; ++j) { hs = fmax2(fabs(g[j]), 1.0) / fmax2(fabs(x[j]), typsiz[j]); if (fabs(a[j + j * nr] - udiag[j]) > fmax2(fabs(udiag[j]), hs) * analtl) { *msg = -22; return; } for (i = j+1; i < n; ++i) { temp1 = a[i + j * nr]; temp2 = fabs(temp1 - a[j + i * nr]); temp1 = fabs(temp1); if (temp2 > fmax2(temp1, hs) * analtl) { *msg = -22; return; } } } } /* heschk */ static int opt_stop(int n, double *xpls, double fpls, double *gpls, double *x, int itncnt, int *icscmx, double gradtl, double steptl, double *sx, double fscale, int itnlim, int iretcd, Rboolean mxtake, int *msg) { /* Unconstrained minimization stopping criteria : * Find whether the algorithm should terminate, due to any * of the following: * 1) problem solved within user tolerance * 2) convergence within user tolerance * 3) iteration limit reached * 4) divergence or too restrictive maximum step (stepmx) suspected * ARGUMENTS : * n --> dimension of problem * xpls(n) --> new iterate x[k] * fpls --> function value at new iterate f(xpls) * gpls(n) --> gradient at new iterate, g(xpls), or approximate * x(n) --> old iterate x[k-1] * itncnt --> current iteration k * icscmx <--> number consecutive steps >= stepmx * [retain value between successive calls] * gradtl --> tolerance at which relative gradient considered close * enough to zero to terminate algorithm * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * sx(n) --> diagonal scaling matrix for x * fscale --> estimate of scale of objective function * itnlim --> maximum number of allowable iterations * iretcd --> return code * mxtake --> boolean flag indicating step of maximum length used * msg --> if msg includes a term 8, suppress output * * VALUE : * `itrmcd' : termination code */ int i, jtrmcd; double d, relgrd, relstp, rgx, rsx; /* last global step failed to locate a point lower than x */ if (iretcd == 1) return 3; /* else : */ /* find direction in which relative gradient maximum. */ /* check whether within tolerance */ d = fmax2(fabs(fpls), fscale); rgx = 0.; for (i = 0; i < n; ++i) { relgrd = fabs(gpls[i]) * fmax2(fabs(xpls[i]), 1./sx[i]) / d; if(rgx < relgrd) rgx = relgrd; } jtrmcd = 1; if (rgx > gradtl) { if (itncnt == 0) return 0; /* find direction in which relative stepsize maximum */ /* check whether within tolerance. */ rsx = 0.; for (i = 0; i < n; ++i) { relstp = fabs(xpls[i] - x[i]) / fmax2(fabs(xpls[i]), 1./sx[i]); if(rsx < relstp) rsx = relstp; } jtrmcd = 2; if (rsx > steptl) { /* check iteration limit */ jtrmcd = 4; if (itncnt < itnlim) { /* check number of consecutive steps \ stepmx */ if (!mxtake) { *icscmx = 0; return 0; } else { ++(*icscmx); if (*icscmx < 5) return 0; jtrmcd = 5; } } } } return jtrmcd; } /* opt_stop */ static void optchk(int n, double *x, double *typsiz, double *sx, double *fscale, double gradtl, int *itnlim, int *ndigit, double epsm, double *dlt, int *method, int *iexp, int *iagflg, int *iahflg, double *stepmx, int *msg) { /* Check input for reasonableness. * Return *msg in {-1,-2,..,-7} if something is wrong * PARAMETERS : * n --> dimension of problem * x(n) --> on entry, estimate to root of fcn * typsiz(n) <--> typical size of each component of x * sx(n) <-- diagonal scaling matrix for x * fscale <--> estimate of scale of objective function fcn * gradtl --> tolerance at which gradient considered close * enough to zero to terminate algorithm * itnlim <--> maximum number of allowable iterations * ndigit <--> number of good digits in optimization function fcn * epsm --> machine epsilon * dlt <--> trust region radius * method <--> algorithm indicator * iexp <--> expense flag * iagflg <--> =1 if analytic gradient supplied * iahflg <--> =1 if analytic hessian supplied * stepmx <--> maximum step size * msg <--> message and error code * ipr --> device to which to send output */ int i; double stpsiz; /* check that parameters only take on acceptable values. if not, set them to default values. */ if (*method < 1 || *method > 3) *method = 1; if (*iagflg != 1) *iagflg = 0; if (*iahflg != 1) *iahflg = 0; if (*iexp != 0) *iexp = 1; if (*msg / 2 % 2 == 1 && *iagflg == 0) { *msg = -6; return;/* 830 write(ipr,906) msg,iagflg */ } if (*msg / 4 % 2 == 1 && *iahflg == 0) { *msg = -7; return;/* 835 write(ipr,907) msg,iahflg */ } /* check dimension of problem */ if (n <= 0) { *msg = -1; return;/* 805 write(ipr,901) n */ } if (n == 1 && *msg % 2 == 0) { *msg = -2; return;/* 810 write(ipr,902) */ } /* compute scale matrix */ for (i = 0; i < n; ++i) { if (typsiz[i] == 0.) typsiz[i] = 1.; else if (typsiz[i] < 0.) typsiz[i] = -typsiz[i]; sx[i] = 1. / typsiz[i]; } /* compute default maximum step size if not provided */ if (*stepmx <= 0.) { stpsiz = 0.; for (i = 0; i < n; ++i) stpsiz += x[i] * x[i] * sx[i] * sx[i]; *stepmx = 1000. * fmax2(sqrt(stpsiz), 1.); } /* check function scale */ if (*fscale == 0.) *fscale = 1.; else if (*fscale < 0.) *fscale = -(*fscale); /* check gradient tolerance */ if (gradtl < 0.) { *msg = -3; return;/* 815 write(ipr,903) gradtl */ } /* check iteration limit */ if (*itnlim <= 0) { *msg = -4; return;/* 820 write(ipr,904) itnlim */ } /* check number of digits of accuracy in function fcn */ if (*ndigit == 0) { *msg = -5; return;/* 825 write(ipr,905) ndigit */ } else if (*ndigit < 0) /* use default, = 15 for IEEE double */ *ndigit = (int) (-log10(epsm)); /* check trust region radius */ if (*dlt <= 0.) { *dlt = -1.; } else if (*dlt > *stepmx) { *dlt = *stepmx; } return; } /* optchk */ static void prt_result(int nr, int n, const double x[], double f, const double g[], const double *a, const double p[], int itncnt, int iflg) { /* * PURPOSE * * Print information on current iteration. * * PARAMETERS * * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> iterate x[k] * f --> function value at x[k] * g(n) --> gradient at x[k] * a(n,n) --> hessian at x[k] * p(n) --> step taken * itncnt --> iteration number k * iflg --> flag controlling info to print */ /* Print iteration number */ Rprintf("iteration = %d\n", itncnt); /* Print step */ if (iflg != 0) { Rprintf("Step:\n"); printRealVector((double *)p, n, 1); } /* Print current iterate */ Rprintf("Parameter:\n"); printRealVector((double *)x, n, 1); /* Print function value */ Rprintf("Function Value\n"); printRealVector((double *)&f, 1, 1); /* Print gradient */ Rprintf("Gradient:\n"); printRealVector((double *)g, n, 1); #ifdef NEVER /* Print Hessian */ /* We don't do this because the printRealMatrix code takes a SEXP rather than a double*. We could do something ugly like use fixed e format but that would be UGLY! */ if (iflg != 0) { } #endif Rprintf("\n"); } /* prt_result */ static void optdrv_end(int nr, int n, double *xpls, double *x, double *gpls, double *g, double *fpls, double f, double *a, double *p, int itncnt, int itrmcd, int *msg, void (*print_result)(int, int, const double *, double, const double *, const double *, const double *, int, int)) { int i; /* termination : reset xpls,fpls,gpls, if previous iterate solution */ if (itrmcd == 3) { *fpls = f; for (i = 0; i < n; ++i) { xpls[i] = x[i]; gpls[i] = g[i]; } } if (*msg / 8 % 2 == 0) (*print_result)(nr, n, xpls, *fpls, gpls, a, p, itncnt, 0); *msg = 0; } /* optdrv_end */ static void optdrv(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn, void *state, double *typsiz, double fscale, int method, int iexp, int *msg, int ndigit, int itnlim, int iagflg, int iahflg, double dlt, double gradtl, double stepmx, double steptl, double *xpls, double *fpls, double *gpls, int *itrmcd, double *a, double *udiag, double *g, double *p, double *sx, double *wrk0, double *wrk1, double *wrk2, double *wrk3, int *itncnt) { /* Driver for non-linear optimization problem -- called by optif0() & optif9() * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> on entry: estimate to a root of fcn * fcn --> name of subroutine to evaluate optimization function * must be declared external in calling routine * fcn: R^n --> R * d1fcn --> (optional) name of subroutine to evaluate gradient * of fcn. must be declared external in calling routine * d2fcn --> (optional) name of subroutine to evaluate * hessian of of fcn. must be declared external * in calling routine * state <--> information other than x and n that fcn, * d1fcn and d2fcn requires. * state is not modified in optdrv (but can be * modified by fcn, d1fcn or d2fcn). * typsiz(n) --> typical size for each component of x * fscale --> estimate of scale of objective function * method --> algorithm to use to solve minimization problem * =1 line search * =2 double dogleg * =3 more-hebdon * iexp --> =1 if optimization function fcn is expensive to * evaluate, =0 otherwise. if set then hessian will * be evaluated by secant update instead of * analytically or by finite differences * msg <--> on input: ( > 0) message to inhibit certain * automatic checks; see do_nlm() in optimize.c * on output: ( < 0) error code; =0 no error * ndigit --> number of good digits in optimization function fcn * itnlim --> maximum number of allowable iterations * iagflg --> =1 if analytic gradient supplied * iahflg --> =1 if analytic hessian supplied * dlt --> trust region radius * gradtl --> tolerance at which gradient considered close * enough to zero to terminate algorithm * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * xpls(n) <--> on exit: xpls is local minimum * fpls <--> on exit: function value at solution, xpls * gpls(n) <--> on exit: gradient at solution xpls * itrmcd <-- termination code * a(n,n) --> workspace for hessian (or estimate) * and its cholesky decomposition * udiag(n) --> workspace [for diagonal of hessian] * g(n) --> workspace (for gradient at current iterate) * p(n) --> workspace for step * sx(n) --> workspace (for diagonal scaling matrix) * wrk0(n) --> workspace * wrk1(n) --> workspace * wrk2(n) --> workspace * wrk3(n) --> workspace * itncnt current iteration, k {{was `internal'}} * internal variables * analtl tolerance for comparison of estimated and * analytical gradients and hessians * epsm machine epsilon * f function value: fcn(x) * rnf relative noise in optimization function fcn. * noise=10.**(-ndigit) */ Rboolean mxtake = FALSE, noupdt; int i, iretcd = 0, icscmx = 0; double dltp = 0., epsm, phip0 = 0., f, analtl; double dlpsav = 0., phisav = 0., dltsav = 0.;/* -Wall */ double amusav = 0., phpsav = 0.; /* -Wall */ double phi = 0., amu = 0., rnf, wrk; *itncnt = 0; epsm = DBL_EPSILON; optchk(n, x, typsiz, sx, &fscale, gradtl, &itnlim, &ndigit, epsm, &dlt, &method, &iexp, &iagflg, &iahflg, &stepmx, msg); if (*msg < 0) return; for (i = 0; i < n; ++i) p[i] = 0.; rnf = Rexp10(-(double)ndigit); rnf = fmax2(rnf, epsm); analtl = sqrt(rnf); analtl = fmax2(0.1, analtl); /* evaluate fcn(x) */ (*fcn)(n, x, &f, state); /* evaluate analytic or finite difference gradient and check analytic gradient, if requested. */ if (!iagflg) { fstofd(1, 1, n, x, (fcn_p)fcn, state, &f, g, sx, rnf, &wrk, 1); } else { /* analytic gradient */ (*d1fcn)(n, x, g, state); if (*msg / 2 % 2 == 0) { grdchk(n, x, (fcn_p)fcn, state, f, g, typsiz, sx, fscale, rnf, analtl, wrk1, msg); if (*msg < 0) return; } } iretcd = -1; *itrmcd = opt_stop(n, x, f, g, wrk1, *itncnt, &icscmx, gradtl, steptl, sx, fscale, itnlim, iretcd, /* mxtake = */FALSE, msg); if (*itrmcd != 0) { optdrv_end(nr, n, xpls, x, gpls, g, fpls, f, a, p, *itncnt, 3, msg, prt_result); return; } if (iexp) { /* if optimization function expensive to evaluate (iexp=1), then * hessian will be obtained by secant updates. get initial hessian.*/ hsnint(nr, n, a, sx, method); } else { /* evaluate analytic or finite difference hessian and check analytic * hessian if requested (only if user-supplied analytic hessian * routine d2fcn fills only lower triangular part and diagonal of a). */ if (!iahflg) { /* no analytic hessian */ if (iagflg) /* anal.gradient */ fstofd(nr, n, n, x, (fcn_p)d1fcn, state, g, a, sx, rnf, wrk1,3); else sndofd(nr, n, x, (fcn_p)fcn, state, f, a, sx, rnf, wrk1, wrk2); } else { /* analytic hessian */ if (*msg / 4 % 2 == 1) { (*d2fcn)(nr, n, x, a, state); } else { heschk(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn, (d2fcn_p)d2fcn, state, f, g, a, typsiz, sx, rnf, analtl, iagflg, udiag, wrk1, wrk2, msg); /* heschk evaluates d2fcn and checks it against the finite * difference hessian which it calculates by calling fstofd * (if iagflg == 1) or sndofd (otherwise). */ if (*msg < 0) return; } } } if (*msg / 8 % 2 == 0) prt_result(nr, n, x, f, g, a, p, *itncnt, 1); /* THE Iterations : */ while(1) { ++(*itncnt); /* find perturbed local model hessian and its LL+ decomposition * ( skip this step if line search or dogstep techniques being used * with secant updates (i.e. method == 1 or 2). * cholesky decomposition L already obtained from secfac.) */ if (iexp && method != 3) { goto L105; } L103: chlhsn(nr, n, a, epsm, sx, udiag); L105: /* solve for newton step: ap=-g */ for (i = 0; i < n; ++i) wrk1[i] = - g[i]; lltslv(nr, n, a, p, wrk1); /* decide whether to accept newton step xpls=x + p */ /* or to choose xpls by a global strategy. */ if (iagflg == 0 && method != 1) { dltsav = dlt; if (method != 2) {/* i.e. method = 3 */ amusav = amu; dlpsav = dltp; phisav = phi; phpsav = phip0; } } switch(method) { case 1: lnsrch(n, x, f, g, p, xpls, fpls, (fcn_p)fcn, state, &mxtake, &iretcd, stepmx, steptl, sx); break; case 2: dogdrv(nr, n, x, f, g, a, p, xpls, fpls, (fcn_p)fcn, state, sx, stepmx, steptl, &dlt, &iretcd, &mxtake, wrk0, wrk1, wrk2, wrk3, itncnt); break; case 3: hookdrv(nr, n, x, f, g, a, udiag, p, xpls, fpls, (fcn_p)fcn, state, sx, stepmx, steptl, &dlt, &iretcd, &mxtake, &amu, &dltp, &phi, &phip0, wrk0, wrk1 , wrk2, epsm, *itncnt); break; } /* if could not find satisfactory step and forward difference */ /* gradient was used, retry using central difference gradient. */ if (iretcd == 1 && iagflg == 0) { iagflg = -1; /* set iagflg for central differences */ fstocd(n, x, (fcn_p)fcn, state, sx, rnf, g); if (method == 1) goto L105; dlt = dltsav; if (method == 2) goto L105; /* else : method == 3 */ amu = amusav; dltp = dlpsav; phi = phisav; phip0 = phpsav; goto L103; } /* calculate step for output */ for (i = 0; i < n; ++i) p[i] = xpls[i] - x[i]; /* calculate gradient at xpls */ switch(iagflg) { case -1: /* central difference gradient */ fstocd(n, xpls, (fcn_p)fcn, state, sx, rnf, gpls); break; case 0: /* forward difference gradient */ fstofd(1, 1, n, xpls, (fcn_p)fcn, state, fpls, gpls, sx, rnf, &wrk, 1); break; default: /* analytic gradient */ (*d1fcn)(n, xpls, gpls, state); } /* check whether stopping criteria satisfied */ *itrmcd = opt_stop(n, xpls, *fpls, gpls, x, *itncnt, &icscmx, gradtl, steptl, sx, fscale, itnlim, iretcd, mxtake, msg); if(*itrmcd != 0) break; /* evaluate hessian at xpls */ if (iexp) { /* expensive obj.fun. */ if (method == 3) secunf(nr, n, x, g, a, udiag, xpls, gpls, epsm, *itncnt, rnf, iagflg, &noupdt, wrk1, wrk2, wrk3); else secfac(nr, n, x, g, a, xpls, gpls, epsm, *itncnt, rnf, iagflg, &noupdt, wrk0, wrk1, wrk2, wrk3); } else { /* iexp == 0 */ if (!iahflg) { if (iagflg) fstofd(nr, n, n, xpls, (fcn_p)d1fcn, state, gpls, a, sx, rnf, wrk1, 3); else /* (iagflg != 1) */ sndofd(nr, n, xpls, (fcn_p)fcn, state, *fpls, a, sx, rnf, wrk1, wrk2); } else /* analytic hessian */ (*d2fcn)(nr, n, xpls, a, state); } if (*msg / 16 % 2 == 1) prt_result(nr, n, xpls, *fpls, gpls, a, p, *itncnt, 1); /* x <-- xpls and g <-- gpls and f <-- fpls */ f = *fpls; for (i = 0; i < n; ++i) { x[i] = xpls[i]; g[i] = gpls[i]; } } /* END while(1) */ optdrv_end(nr, n, xpls, x, gpls, g, fpls, f, a, p, *itncnt, *itrmcd, msg, prt_result); } /* optdrv */ #if 0 static void dfault(int n, double *x, double *typsiz, double *fscale, int *method, int *iexp, int *msg, int *ndigit, int *itnlim, int *iagflg, int *iahflg, double *dlt, double *gradtl, double *stepmx, double *steptl) { /* Set default values for each input variable to minimization algorithm * for optif0() only. * PARAMETERS : * INPUT: * n dimension of problem * x(n) initial guess to solution (to compute max step size) * OUTPUT: * typsiz(n) typical size for each component of x * fscale estimate of scale of minimization function * method algorithm to use to solve minimization problem * iexp =0 if minimization function not expensive to evaluate * msg message to inhibit certain automatic checks + output * ndigit number of good digits in minimization function * itnlim maximum number of allowable iterations * iagflg =0 if analytic gradient not supplied * iahflg =0 if analytic hessian not supplied * dlt trust region radius * gradtl tolerance at which gradient considered close enough * to zero to terminate algorithm * stepmx value of zero to trip default maximum in optchk * steptl tolerance at which successive iterates considered * close enough to terminate algorithm */ double epsm; int i; /* set typical size of x and minimization function */ for (i = 0; i < n; ++i) typsiz[i] = 1.; *fscale = 1.; /* set tolerances */ epsm = DBL_EPSILON; /* for IEEE : = 2^-52 ~= 2.22 e-16 */ *gradtl = pow(epsm, 1./3.); /* for IEEE : = 2^(-52/3) ~= 6.055 e-6 */ *steptl = sqrt(epsm); /* for IEEE : = 2^-26 ~= 1.490 e-8 */ *stepmx = 0.;/* -> compute default in optchk() */ *dlt = -1.;/* (not needed for method 1) */ /* set flags */ *method = 1; *iexp = 1; *msg = 0; *ndigit = -1;/* -> compute default = floor(-log10(EPS)) in optchk() */ *itnlim = 150; *iagflg = 0;/* no gradient */ *iahflg = 0;/* no hessian */ } /* dfault() */ void optif0(int nr, int n, double *x, fcn_p fcn, void *state, double *xpls, double *fpls, double *gpls, int *itrmcd, double *a, double *wrk) { /* Provide simplest interface to minimization package. * User has no control over options. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> initial estimate of minimum * fcn --> name of routine to evaluate minimization function. * must be declared external in calling routine. * state <--> information other than x and n that fcn requires * state is not modified in optif0 (but can be * modified by fcn). * xpls(n) <-- local minimum * fpls <-- function value at local minimum xpls * gpls(n) <-- gradient at local minimum xpls * itrmcd <-- termination code * a(n,n) --> workspace * wrk(n,9) --> workspace */ int iexp, iagflg, iahflg; int ndigit, method, itnlim, itncnt; double fscale, gradtl, steptl, stepmx, dlt; int msg; /* Function Body */ dfault(n, x, &wrk[nr], &fscale, &method, &iexp, &msg, &ndigit, &itnlim, &iagflg, &iahflg, &dlt, &gradtl, &stepmx, &steptl); optdrv(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn_dum, (d2fcn_p)d2fcn_dum, state, &wrk[nr * 3], fscale, method, iexp, &msg, ndigit, itnlim, iagflg, iahflg, dlt, gradtl, stepmx, steptl, xpls, fpls, gpls, itrmcd, a, wrk, &wrk[nr], &wrk[nr * 2], &wrk[nr * 4], &wrk[nr * 5], &wrk[nr * 6], &wrk[nr * 7], &wrk[nr * 8], &itncnt); } /* optif0 */ #endif /* ---- this one is called from optimize.c : --------------- */ void optif9(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn, void *state, double *typsiz, double fscale, int method, int iexp, int *msg, int ndigit, int itnlim, int iagflg, int iahflg, double dlt, double gradtl, double stepmx, double steptl, double *xpls, double *fpls, double *gpls, int *itrmcd, double *a, double *wrk, int *itncnt) { /* provide complete interface to minimization package. * user has full control over options. * PARAMETERS : * nr --> row dimension of matrix * n --> dimension of problem * x(n) --> on entry: estimate to a root of fcn * fcn --> name of subroutine to evaluate optimization function * must be declared external in calling routine * fcn: r(n) --> r(1) * d1fcn --> (optional) name of subroutine to evaluate gradient * of fcn. must be declared external in calling routine * d2fcn --> (optional) name of subroutine to evaluate hessian of * of fcn. must be declared external in calling routine * state <--> information other than x and n that fcn, * d1fcn and d2fcn requires. * state is not modified in optif9 (but can be * modified by fcn, d1fcn or d2fcn). * typsiz(n) --> typical size for each component of x * fscale --> estimate of scale of objective function * method --> algorithm to use to solve minimization problem * =1 line search * =2 double dogleg * =3 more-hebdon * iexp --> =1 if optimization function fcn is expensive to * evaluate, =0 otherwise. if set then hessian will * be evaluated by secant update instead of * analytically or by finite differences * msg <--> on input: ( > 0) to inhibit certain automatic checks * on output: ( < 0) error code; =0 no error * ndigit --> number of good digits in optimization function fcn * itnlim --> maximum number of allowable iterations * iagflg --> =1 if analytic gradient supplied * iahflg --> =1 if analytic hessian supplied * dlt --> trust region radius * gradtl --> tolerance at which gradient considered close * enough to zero to terminate algorithm * stepmx --> maximum allowable step size * steptl --> relative step size at which successive iterates * considered close enough to terminate algorithm * xpls(n) <--> on exit: xpls is local minimum * fpls <--> on exit: function value at solution, xpls * gpls(n) <--> on exit: gradient at solution xpls * itrmcd <-- termination code (in 0..5 ; 0 is "perfect"); * see optcode() in optimize.c for meaning * a(n,n) --> workspace for hessian (or estimate) * and its cholesky decomposition * wrk(n,8) --> workspace * itncnt <--> iteration count */ optdrv(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn, (d2fcn_p)d2fcn, state, typsiz, fscale, method, iexp, msg, ndigit, itnlim, iagflg, iahflg, dlt, gradtl, stepmx, steptl, xpls, fpls, gpls, itrmcd, a, wrk, wrk + nr, wrk + nr * 2, wrk + nr * 3, wrk + nr * 4, wrk + nr * 5, wrk + nr * 6, wrk + nr * 7, itncnt); } /* optif9 */