Google Git
Sign in
third-party-mirror / R / 116db572fde53bb7b4e7596511f5804c5746bc81 / . / src / appl / uncmin.c
blob: f911c58fcfb2f3f61f54c86b08a569273bd2ebe0 [file] [log] [blame]
/*
* R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1999-2017 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 <math.h>
#include <float.h> /* DBL_MAX */
#include <R_ext/Applic.h>
#include <R_ext/Boolean.h>
#include <R_ext/Print.h> /* Rprintf */
#include <R_ext/PrtUtil.h> /* printRealVector */
#include <R_ext/Linpack.h> /* ddot, dnrm2, dtrsl, dscal */
#include <Rmath.h>
#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 ``choleski + 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 */