The R Project SVN
R
(root)/branches/R-grid/src/appl/uncmin.c – Rev 26266
Rev 24068 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
/** R : A Computer Language for Statistical Data Analysis* Copyright (C) 1997-2001 Saikat DebRoy and the* R Development Core Team** 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, write to the Free Software* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307* 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 <math.h>#include <R_ext/Applic.h>#include <R_ext/Boolean.h>#include <R_ext/PrtUtil.h> /* printRealVector */#include <R_ext/Linpack.h> /* ddot, dnrm2, dtrsl, dscal */#define MATHLIB_PRIVATE#include <Rmath.h> /* for dimach */#undef MATHLIB_PRIVATE/* CC subroutines mvmlt[lsu] should be REPLACED by BLAS ones!* CC* CC--- choldc(nr,n,a,diagmx,tol,addmax) is ``choleski + tolerance''* CC ------* CC it should make use of BLAS routines as [linkpack's dpofa!] */#ifndef HAVE_HYPOT# define hypot pythag#endifvoid 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 = pow(10.0, -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 */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. */}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 voidcholdc(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) {/* 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) {a[i + i * nr] = sqrt(tmp1);}else {/* 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;}/* find i,j element of lower triangular matrix */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];}}} /* 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 voidqrupdt(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 voidtregup(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 optimizeevaluation of next "if" statement (in which case fplsp could beundefined).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 voidlnsrch(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 */*mxtake = FALSE;*iretcd = 2;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. */while(*iretcd > 1) {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;elselambda = tlmbda;}}}} /* lnsrch */static voiddog_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*pwhich 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 voiddogdrv(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, rnwtln, eta, cln, tmp;*iretcd = 4;fstdog = TRUE;tmp = 0.;for (i = 0; i < n; ++i)tmp += sx[i] * sx[i] * p[i] * p[i];rnwtln = sqrt(tmp);while(*iretcd > 1) {/* 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);}} /* dogdrv */static voidhook_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)] */for (i = 0; i < n; ++i) {a[i + i * nr] = udiag[i] + *amu * sx[i] * sx[i];for (j = 0; j < i; ++i)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 voidhookdrv(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, rnwtln, tmp;*iretcd = 4;fstime = TRUE;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;}}while(*iretcd > 1) {/* 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);}} /* hookdrv */static voidsecunf(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 tolower 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 voidsecfac(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);elsereltol = 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 voidchlhsn(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.;elseamu = 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 voidhsnint(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];elsea[i + i * nr] = sx[i];for (j = 0; j < i; ++j)a[i + j * nr] = 0.;}} /* hsnint */static voidfstofd(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 aeach 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 directionof 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 voidgrdchk(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 gradientand 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 voidheschk(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);elsesndofd(nr, n, x, (fcn_p)fcn, state, f, a, sx, rnf, wrk1, wrk2);/* copy lower triangular part of "a" to upper triangular partand 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 */staticint 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 voidoptchk(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 voidprt_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 printRealMatrixcode takes a SEXP rather than a double*.We could do something ugly like use fixed e formatbut that would be UGLY!*/if (iflg != 0) {}#endifRprintf("\n");} /* prt_result */static voidoptdrv_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 voidoptdrv(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 ../main/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, noupdt;int i, iretcd, icscmx;double dltp, epsm, phip0, f, analtl;double dlpsav = 0., phisav = 0., dltsav = 0.;/* -Wall */double amusav = 0., phpsav = 0.; /* -Wall */double phi, amu, rnf, wrk;for (i = 0; i < n; ++i)p[i] = 0.;*itncnt = 0;iretcd = -1;epsm = d1mach(4);optchk(n, x, typsiz, sx, &fscale, gradtl, &itnlim, &ndigit, epsm,&dlt, &method, &iexp, &iagflg, &iahflg, &stepmx, msg);if (*msg < 0) return;rnf = pow(10., -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 andcheck 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;}}*itrmcd = opt_stop(n, x, f, g, wrk1, *itncnt, &icscmx,gradtl, steptl, sx, fscale, itnlim, iretcd, mxtake,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);elsesndofd(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);elsesecfac(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 */static voiddfault(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 = d1mach(4); /* 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() */voidoptif0(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 *//* ---- this one is called from ../main/optimize.c : --------------- */voidoptif9(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 ../main/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 */