/*
 *  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>
// as in <Defn.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 ``cholesky + tolerance''
 * CC	 ------
 * CC	 it should make use of BLAS routines as [linkpack's dpofa!] */


void fdhess(int n, double *x, double fval, fcn_p fun, void *state,
	    double *h, int nfd, double *step, double *f,
	    int ndigit, double *typx)
{
/*	calculates a numerical approximation to the upper triangular
 *	portion of the second derivative matrix (the hessian).

 *	Algorithm A5.6.2 from Dennis and Schnabel (1983), numerical methods
 *	for unconstrained optimization and nonlinear equations,
 *	prentice-hall, 321-322.

 *	programmed by richard h. jones, january 11, 1989

 * INPUT to subroutine

 *	n      the number of parameters
 *	x      vector of parameter values
 *	fval   double precision value of function at x
 *	fun    a function provided by the user which must be declared as
 *	       external in the calling program.	 its call must
 *	       be of the call fun(n,x,state,fval) where fval is the
 *	       computed value of the function
 *	state  information other than x and n that fun requires.
 *	       state is not modified in fdhess (but can be modified by fun).
 *	nfd    first dimension of h in the calling program

 * OUTPUT from subroutine

 *	h     an n by n matrix of the approximate hessian

 * Work space :

 *	step	a real array of length n
 *	f	a double precision array of length n
 */

    int i, j;
    double tempi, tempj, fii, eta, fij;

    eta = Rexp10(-ndigit/3.0);
    for (i = 0; i < n; ++i) {
	step[i] = eta * fmax2(x[i], typx[i]);
	if (typx[i] < 0.)
	    step[i] = -step[i];

	tempi = x[i];
	x[i] += step[i];
	step[i] = x[i] - tempi;
	(*fun)(n, x, &f[i], state);
	x[i] = tempi;
    }
    for (i = 0; i < n; ++i) {
	tempi = x[i];
	x[i] += step[i] * 2.;
	(*fun)(n, x, &fii, state);
	h[i + i * nfd] = (fval - f[i] + (fii - f[i]))/(step[i] * step[i]);
	x[i] = tempi + step[i];
	for (j = i + 1; j < n; ++j) {
	    tempj = x[j];
	    x[j] += step[j];
	    (*fun)(n, x, &fij, state);
	    h[i + j * nfd] = (fval - f[i] + (fij - f[j]))/(step[i] * step[j]);
	    x[j] = tempj;
	}
	x[i] = tempi;
    }
} /* fdhess */

#if 0
static void d1fcn_dum(int n, double *x, double *g, void *state)
{
/*	dummy routine to prevent unsatisfied external diagnostic
 *	when specific analytic gradient function not supplied. */
}

static void d2fcn_dum(int nr, int n, double *x, double *h, void *state)
{
/*	dummy routine to prevent unsatisfied external diagnostic
 *	when specific analytic hessian function not supplied. */
}
#endif

static void mvmltl(int nr, int n, double *a, double *x, double *y)
{
/*	compute y = l x
 *	where l is a lower triangular matrix stored in a

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	     --> lower triangular (n*n) matrix
 *	x(n)	     --> operand vector
 *	y(n)	    <--	 result vector

 *	note

 *	x and y cannot share storage */

  int i, j;
  double sum;

  for (i = 0; i < n; ++i) {
    sum = 0.;
    for (j = 0; j <= i; ++j)
      sum += a[i + j * nr] * x[j];
    y[i] = sum;
  }
} /* mvmltl */

static void mvmltu(int nr, int n, double *a, double *x, double *y)
{
/*	compute y = (L+) x
 *	where L is a lower triangular matrix stored in a
 *	(L-transpose (L+) is taken implicitly)

 * ARGUMENTS :

 *	nr	 --> row dimension of matrix
 *	n	 --> dimension of problem
 *	a(nr,1)	 --> lower triangular (n*n) matrix
 *	x(n)	 --> operand vector

 *	y(n)	<-- result vector

 * NOTE :	x and y cannot share storage */

  int i, length, one = 1;

  for (i = 0, length = n; i < n; --length, ++i)
    y[i] = F77_CALL(ddot)(&length, &a[i + i * nr], &one, &x[i], &one);
} /* mvmltu */


static void mvmlts(int nr, int n, double *a, double *x, double *y)
{
/*	compute y=ax
 *	where "a" is a symmetric (n*n) matrix stored in its lower
 *	triangular part and x,y are n-vectors

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	     --> symmetric (n*n) matrix stored in
 *			 lower triangular part and diagonal
 *	x(n)	     --> operand vector
 *	y(n)	    <--	 result vector

 * NOTE:	x and y cannot share storage.
 */

    int i, j;
    double sum;

    for (i = 0; i < n; ++i) {
	sum = 0.;
	for (j = 0; j <= i; ++j)
	    sum += a[i + j * nr] * x[j];
	for (j = i+1; j < n; ++j)
	    sum += a[j + i * nr] * x[j];
	y[i] = sum;
    }
} /* mvmlts */

static void lltslv(int nr, int n, double *a, double *x, double *b)
{
/*	solve ax=b where a has the form l(l-transpose)
 *	but only the lower triangular part, l, is stored.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	     --> matrix of form l(l-transpose).
 *			 on return a is unchanged.
 *	x(n)	    <--	 solution vector
 *	b(n)	     --> right-hand side vector

 *	note

 *	if b is not required by calling program, then
 *	b and x may share the same storage. */

    int job = 0, info;

    if (x != b) Memcpy(x, b, n);
    F77_CALL(dtrsl)(a, &nr, &n, x, &job, &info);
    job = 10;
    F77_CALL(dtrsl)(a, &nr, &n, x, &job, &info);
} /* lltslv */

static void
choldc(int nr, int n, double *a, double diagmx, double tol, double *addmax)
{
/* Find the perturbed L(L-transpose) [written LL+] decomposition
 * of a+d, where d is a non-negative diagonal matrix added to a if
 * necessary to allow the cholesky decomposition to continue.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	    <--> on entry: matrix for which to find perturbed
 *			      cholesky decomposition
 *			 on exit:  contains L of LL+ decomposition
 *			 in lower triangular part and diagonal of "a"
 *	diagmx	     --> maximum diagonal element of "a"
 *	tol	     --> tolerance
 *	addmax	    <--	 maximum amount implicitly added to diagonal of "a"
 *			 in forming the cholesky decomposition of a+d
 *	internal variables

 *	aminl	 smallest element allowed on diagonal of L
 *	amnlsq	 =aminl**2
 *	offmax	 maximum off-diagonal element in column of a


 *	description

 *	the normal cholesky decomposition is performed.	 however, if at any
 *	point the algorithm would attempt to set L(i,i)=sqrt(temp)
 *	with temp < tol*diagmx, then L(i,i) is set to sqrt(tol*diagmx)
 *	instead.  this is equivalent to adding tol*diagmx-temp to a(i,i)
 */

    double tmp1, tmp2;
    int i, j, k;
    double aminl, offmax, amnlsq;
    double sum;

    *addmax = 0.0;
    aminl = sqrt(diagmx * tol);
    amnlsq = aminl * aminl;

    /*	form row i of L */

    for (i = 0; i < n; ++i) {

	// A[i,j] := * || find i,j element of lower triangular matrix L
	for (j = 0; j < i; ++j) {
	    sum = 0.;
	    for (k = 0; k < j; ++k)
		sum += a[i + k * nr] * a[j + k * nr];
	    a[i + j * nr] = (a[i + j * nr] - sum) / a[j + j * nr];
	}

	// A[i,i] := * || find diagonal elements of L
	sum = 0.;
	for (k = 0; k < i; ++k)
	    sum += a[i + k * nr] * a[i + k * nr];

	tmp1 = a[i + i * nr] - sum;
	if (tmp1 >= amnlsq) { // normal Cholesky
	    a[i + i * nr] = sqrt(tmp1);
	}
	else { // augment diagonal of L
	    /*	find maximum off-diagonal element in row */
	    offmax = 0.;
	    for (j = 0; j < i; ++j) {
		if(offmax < (tmp2 = fabs(a[i + j * nr])))
		    offmax = tmp2;
	    }
	    if (offmax <= amnlsq) offmax = amnlsq;

	    /* add to diagonal element to
	     * allow cholesky decomposition to continue */
	    a[i + i * nr] = sqrt(offmax);
	    if(*addmax < (tmp2 = offmax - tmp1)) *addmax = tmp2;
	}
    }
} /* choldc */

static void qraux1(int nr, int n, double *r, int i)
{
/* Interchange rows i,i+1 of the upper hessenberg matrix r, columns i to n .

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of matrix
 *	r[n*n]	    <--> upper hessenberg matrix
 *	i	     --> index of row to interchange (i < n-1)
*/
  double tmp;
  double *r1, *r2;
  /* pointer arithmetic : */
  r1 = r + i + i * nr;
  r2 = r1 + 1;

  while(n-- > i) {
    tmp = *r1; *r1 = *r2; *r2 = tmp;
    r1 += nr;
    r2 += nr;
  }
} /* qraux1 */

static void qraux2(int nr, int n, double *r, int i, double a, double b)
{
/* Pre-multiply r by the jacobi rotation j(i,i+1,a,b) .

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of matrix
 *	r(n,n)	    <--> upper hessenberg matrix
 *	i	     --> index of row
 *	a	     --> scalar
 *	b	     --> scalar */

  double c, s;
  double y, z, den;
  double *r1, *r2;

  den = hypot(a,b);
  c = a / den;
  s = b / den;

  /* pointer arithmetic : */
  r1 = r + i + i*nr;
  r2 = r1 + 1;

  while(n-- > i) {
    y = *r1;
    z = *r2;
    *r1 = c * y - s * z;
    *r2 = s * y + c * z;
    r1 += nr;
    r2 += nr;
  }
} /* qraux2 */

static void
qrupdt(int nr, int n, double *a, double *u, double *v)
{
/* Find an orthogonal (n*n) matrix (q*) and an upper triangular (n*n)
 * matrix (r*) such that (q*)(r*)=r+u(v+)

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	    <--> on input:  contains r
 *			 on output: contains (r*)
 *	u(n)	     --> vector
 *	v(n)	     --> vector */

    int i, j, k;
    double t1, t2;
    int ii;

    /*	determine last non-zero in u(.) */

    for(k = n-1; k > 0 && u[k] == 0.0; k--)
	;

    /*	(k-1) jacobi rotations transform
     *	    r + u(v+) --> (r*) + (u(1)*e1)(v+)
     *	which is upper hessenberg */

    if (k > 0) {
	ii = k;
	while(ii > 0) {
	    i = ii - 1;
	    if (u[i] == 0.0) {
		qraux1(nr, n, a, i);
		u[i] = u[ii];
	    } else {
		qraux2(nr, n, a, i, u[i], -u[ii]);
		u[i] = hypot(u[i], u[ii]);
	    }
	    ii = i;
	}
    }

    /*	r <-- r + (u(1)*e1)(v+) */

    for (j = 0; j < n; ++j)
	a[j * nr] += u[0] * v[j];

    /*	(k-1) jacobi rotations transform upper hessenberg r
     *	to upper triangular (r*) */

    for (i = 0; i < k; ++i) {
	if (a[i + i * nr] == 0.)
	    qraux1(nr, n, a, i);
	else {
	    t1 = a[i + i * nr];
	    t2 = -a[i + 1 + i * nr];
	    qraux2(nr, n, a, i, t1, t2);
	}
    }
} /* qrupdt */

static void
tregup(int nr, int n, double *x, double f, double *g, double *a, fcn_p fcn,
       void *state, double *sc, double *sx, Rboolean nwtake,
       double stepmx, double steptl, double *dlt, int *iretcd,
       double *xplsp, double *fplsp, double *xpls, double *fpls,
       Rboolean *mxtake,
       int method, double *udiag)
{
/* TRust REGion UPdating
 * ==	  ==	==
 * Decide whether to accept xpls = x+sc as the next iterate and
 * update the trust region radius dlt.
 * Used iff method == 2 or 3

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate x[k-1]
 *	f	     --> function value at old iterate, f(x)
 *	g(n)	     --> gradient at old iterate, g(x), or approximate
 *	a(n,n)	     --> cholesky decomposition of hessian in
 *			 lower triangular part and diagonal.
 *			 hessian or approx in upper triangular part
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in tregup (but can be
 *			 modified by fcn).
 *	sc(n)	     --> current step
 *	sx(n)	     --> diagonal scaling matrix for x
 *	nwtake	     --> boolean, = TRUE if newton step taken
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	dlt	    <--> trust region radius
 *	iretcd	    <--> return code
 *			   =0 xpls accepted as next iterate;
 *			      dlt trust region for next iteration.
 *			   =1 xpls unsatisfactory but accepted as next iterate
 *			      because xpls-x < smallest allowable step length.
 *			   =2 f(xpls) too large.  continue current iteration
 *			      with new reduced dlt.
 *			   =3 f(xpls) sufficiently small, but quadratic model
 *			      predicts f(xpls) sufficiently well to continue
 *			      current iteration with new doubled dlt.
 *	xplsp(n)    <--> workspace [value needs to be retained between
 *			 succesive calls of k-th global step]
 *	fplsp	    <--> [retain value between successive calls]
 *	xpls(n)	    <--	 new iterate x[k]
 *	fpls	    <--	 function value at new iterate, f(xpls)
 *	mxtake	    <--	 boolean flag indicating step of maximum length used
 *	ipr	     --> device to which to send output
 *	method	     --> algorithm to use to solve minimization problem
 *			   =1 line search
 *			   =2 double dogleg
 *			   =3 more-hebdon
 *	udiag(n)     --> diagonal of hessian in a(.,.) */

    double dltf;
    double temp1;
    int i, j, one = 1;
    double dltfp, dltmp;
    double rln, slp;

    *mxtake = FALSE;
    for (i = 0; i < n; ++i)
	xpls[i] = x[i] + sc[i];

    (*fcn)(n, xpls, fpls, state);
    dltf = *fpls - f;
    slp = F77_CALL(ddot)(&n, g, &one, sc, &one);

    /*	next statement added for case of compilers which do not optimize
	evaluation of next "if" statement (in which case fplsp could be
	undefined).
	if (*iretcd == 4) {
	*fplsp = 0.;
	}
    */
    if (*iretcd == 3 && (*fpls >= *fplsp || dltf > slp * 1e-4)) {
	/*	reset xpls to xplsp and terminate global step */

	*iretcd = 0;
	for (i = 0; i < n; ++i)
	    xpls[i] = xplsp[i];

	*fpls = *fplsp;
	*dlt *= .5;
    }
    else {

	/*	fpls too large */

	if (dltf > slp * 1e-4) {
	    rln = 0.;
	    for (i = 0; i < n; ++i) {
		temp1 = fabs(sc[i])/fmax2(fabs(xpls[i]), 1./sx[i]);
		if(rln < temp1) rln = temp1;
	    }
	    if (rln < steptl) {
		/* cannot find satisfactory xpls sufficiently distinct from x */
		*iretcd = 1;
	    }
	    else {
		/*  reduce trust region and continue global step */
		*iretcd = 2;
		dltmp = -slp * *dlt / ((dltf - slp) * 2.);
		if (dltmp < *dlt * .1)
		    *dlt *= .1;
		else
		    *dlt = dltmp;
	    }
	}
	else {
	    /* fpls sufficiently small */

	    dltfp = 0.;
	    if (method == 2) {
		for (i = 0; i < n; ++i) {
		    temp1 = 0.;
		    for (j = i; j < n; ++j)
			temp1 += a[j + i * nr] * sc[j];
		    dltfp += temp1 * temp1;
		}
	    }
	    else { /* method != 2 */
		for (i = 0; i < n; ++i) {
		    dltfp += udiag[i] * sc[i] * sc[i];
		    temp1 = 0.;
		    for (j = i+1; j < n; ++j)
			temp1 += a[i + j * nr] * sc[i] * sc[j];
		    dltfp += temp1 * 2.;
		}
	    }
	    dltfp = slp + dltfp / 2.;
	    if (*iretcd != 2 && fabs(dltfp - dltf) <= fabs(dltf) * 0.1
		&& nwtake && *dlt <= stepmx * .99) {

		/* double trust region and continue global step */

		*iretcd = 3;
		for (i = 0; i < n; ++i)
		    xplsp[i] = xpls[i];
		*fplsp = *fpls;
		temp1 = *dlt * 2.0;
		*dlt = fmin2(temp1, stepmx);
	    }
	    else {
		/* accept xpls as next iterate.	 choose new trust region. */

		*iretcd = 0;
		if (*dlt > stepmx * .99)
		    *mxtake = TRUE;
		if (dltf >= dltfp * .1) {
		    /* decrease trust region for next iteration */

		    *dlt *= .5;
		}
		else {
		    /* check whether to increase trust region for next iteration */
		    if (dltf <= dltfp * .75) {
			temp1 = *dlt * 2.0;
			*dlt = fmin2(temp1, stepmx);
		    }
		}
	    }
	}
    }
} /* tregup */

static void
lnsrch(int n, double *x, double f, double *g, double *p, double *xpls,
       double *fpls, fcn_p fcn, void *state, Rboolean *mxtake, int *iretcd,
       double stepmx, double steptl, double *sx)
{
/* Find a next newton iterate by line search.  (iff  method == 1)

 * PARAMETERS :

 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate:	x[k-1]
 *	f	     --> function value at old iterate, f(x)
 *	g(n)	     --> gradient at old iterate, g(x), or approximate
 *	p(n)	     --> non-zero newton step
 *	xpls(n)	    <--	 new iterate x[k]
 *	fpls	    <--	 function value at new iterate, f(xpls)
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in lnsrch (but can be
 *			 modified by fcn).
 *	iretcd	    <--	 return code
 *	mxtake	    <--  boolean flag indicating step of maximum length used
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	sx(n)	     --> diagonal scaling matrix for x

 *	internal variables

 *	sln		 newton length
 *	rln		 relative length of newton step
*/

    int i, one = 1;
    Rboolean firstback = TRUE;
    double disc;
    double a3, b;
    double t1, t2, t3, lambda, tlmbda, rmnlmb;
    double scl, rln, sln, slp;
    double temp1;
    double pfpls = 0., plmbda = 0.; /* -Wall */

    temp1 = 0.;
    for (i = 0; i < n; ++i)
	temp1 += sx[i] * sx[i] * p[i] * p[i];
    sln = sqrt(temp1);
    if (sln > stepmx) {
	/*	newton step longer than maximum allowed */
	scl = stepmx / sln;
	F77_CALL(dscal)(&n, &scl, p, &one);
	sln = stepmx;
    }
    slp = F77_CALL(ddot)(&n, g, &one, p, &one);
    rln = 0.;
    for (i = 0; i < n; ++i) {
	temp1 = fabs(p[i])/ fmax2(fabs(x[i]), 1./sx[i]);
	if(rln < temp1) rln = temp1;
    }
    rmnlmb = steptl / rln;
    lambda = 1.0;

    /*	check if new iterate satisfactory.  generate new lambda if necessary. */

    *mxtake = FALSE;
    *iretcd = 2;
    do {
	for (i = 0; i < n; ++i)
	    xpls[i] = x[i] + lambda * p[i];
	(*fcn)(n, xpls, fpls, state);
	if (*fpls <= f + slp * 1e-4 * lambda) { /* solution found */

	    *iretcd = 0;
	    if (lambda == 1. && sln > stepmx * .99) *mxtake = TRUE;
	    return;
	}
	/* else : solution not (yet) found */

	/* First find a point with a finite value */

	if (lambda < rmnlmb) {
	    /* no satisfactory xpls found sufficiently distinct from x */

	    *iretcd = 1;
	    return;
	}
	else { /*	calculate new lambda */

	    /* modifications by BDR 2000/01/05 to cover non-finite values
	     * ">=" instead of "==" :  MM 2001/07/24 */
	    if (*fpls >= DBL_MAX) {
		lambda *= 0.1;
		firstback = TRUE;
	    }
	    else {
		if (firstback) { /*	first backtrack: quadratic fit */
		    tlmbda = -lambda * slp / ((*fpls - f - slp) * 2.);
		    firstback = FALSE;
		}
		else { /*	all subsequent backtracks: cubic fit */
		    t1 = *fpls - f - lambda * slp;
		    t2 = pfpls - f - plmbda * slp;
		    t3 = 1. / (lambda - plmbda);
		    a3 = 3. * t3 * (t1 / (lambda * lambda)
				   - t2 / (plmbda * plmbda));
		    b = t3 * (t2 * lambda / (plmbda * plmbda)
			      - t1 * plmbda / (lambda * lambda));
		    disc = b * b - a3 * slp;
		    if (disc > b * b)
			/* only one positive critical point, must be minimum */
			tlmbda = (-b + ((a3 < 0)? -sqrt(disc): sqrt(disc))) /a3;
		    else
			/* both critical points positive, first is minimum */
			tlmbda = (-b + ((a3 < 0)? sqrt(disc): -sqrt(disc))) /a3;

		    if (tlmbda > lambda * .5)
			tlmbda = lambda * .5;
		}
		plmbda = lambda;
		pfpls = *fpls;
		if (tlmbda < lambda * .1)
		    lambda *= .1;
		else
		    lambda = tlmbda;
	    }
	}
    } while(*iretcd > 1);
} /* lnsrch */

static void
dog_1step(int nr, int n, double *g, double *a, double *p, double *sx,
       double rnwtln, double *dlt, Rboolean *nwtake, Rboolean *fstdog,
       double *ssd, double *v, double *cln, double *eta, double *sc,
       double stepmx)
{
/* Find new step by double dogleg algorithm  (iff method == 2);
 * repeatedly called by dogdrv() only.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	g(n)	     --> gradient at current iterate, g(x)
 *	a(n,n)	     --> cholesky decomposition of hessian in
 *			 lower part and diagonal
 *	p(n)	     --> newton step
 *	sx(n)	     --> diagonal scaling matrix for x
 *	rnwtln	     --> newton step length
 *	dlt	    <--> trust region radius
 *	nwtake	    <--> boolean, =.true. if newton step taken
 *	fstdog	    <--> boolean, =.true. if on first leg of dogleg
 *	ssd(n)	    <--> workspace [cauchy step to the minimum of the
 *			 quadratic model in the scaled steepest descent
 *			 direction] [retain value between successive calls]
 *	v(n)	    <--> workspace  [retain value between successive calls]
 *	cln	    <--> cauchy length
 *			 [retain value between successive calls]
 *	eta		 [retain value between successive calls]
 *	sc(n)	    <--	 current step
 *	ipr	     --> device to which to send output
 *	stepmx	     --> maximum allowable step size

 *	internal variables

 *	cln		 length of cauchy step */

    int i, j, one = 1;
    double alam, bet, alpha, tmp, dot1, dot2;

    /*	can we take newton step */
    *nwtake = (rnwtln <= *dlt);

    if (*nwtake) {
	for (i = 0; i < n; ++i)
	    sc[i] = p[i];
	*dlt = rnwtln;
	return;
    }

    /* else *nwtake = FALSE :
     * newton step too long -- cauchy step is on double dogleg curve */

    if (*fstdog) {
	/*	  calculate double dogleg curve (ssd) */
	*fstdog = FALSE;
	alpha = 0.;
	for (i = 0; i < n; ++i)
	    alpha += g[i] * g[i] / (sx[i] * sx[i]);
	bet = 0.;
	for (i = 0; i < n; ++i) {
	    tmp = 0.;
	    for (j = i; j < n; ++j)
		tmp += a[j + i * nr] * g[j] / (sx[j] * sx[j]);
	    bet += tmp * tmp;
	}
	for (i = 0; i < n; ++i)
	    ssd[i] = -(alpha / bet) * g[i] / sx[i];
	*cln = alpha * sqrt(alpha) / bet;
	*eta = (.8 * alpha * alpha /
		(-bet * F77_CALL(ddot)(&n, g, &one, p, &one))) + .2;
	for (i = 0; i < n; ++i)
	    v[i] = *eta * sx[i] * p[i] - ssd[i];
	if (*dlt == -1.) *dlt = fmin2(*cln, stepmx);
    }

    if (*eta * rnwtln <= *dlt) {
	/*	  take partial step in newton direction */
	for (i = 0; i < n; ++i)
	    sc[i] = *dlt / rnwtln * p[i];
    }
    else if (*cln >= *dlt) {
	/*	    take step in steepest descent direction */
	for (i = 0; i < n; ++i)
	    sc[i] = *dlt / *cln * ssd[i] / sx[i];
    }
    else {
	/* calculate convex combination of ssd and eta*p
	   which has scaled length dlt */
	dot1 = F77_CALL(ddot)(&n, v, &one, ssd, &one);
	dot2 = F77_CALL(ddot)(&n, v, &one, v, &one);
	alam = (-dot1 + sqrt(dot1 * dot1 - dot2 * (*cln * *cln - *dlt * *dlt)))
	    / dot2;
	for (i = 0; i < n; ++i)
	    sc[i] = (ssd[i] + alam * v[i]) / sx[i];
    }
} /* dog_1step */

static void
dogdrv(int nr, int n, double *x, double f, double *g, double *a, double *p,
       double *xpls, double *fpls, fcn_p fcn, void *state, double *sx,
       double stepmx, double steptl, double *dlt, int *iretcd, Rboolean *mxtake,
       double *sc, double *wrk1, double *wrk2, double *wrk3, int *itncnt)
{
/* Find a next newton iterate (xpls) by the double dogleg method
 * (iff	 method == 2 ).

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate x[k-1]
 *	f	     --> function value at old iterate, f(x)
 *	g(n)	     --> gradient  at old iterate, g(x), or approximate
 *	a(n,n)	     --> cholesky decomposition of hessian
 *			 in lower triangular part and diagonal
 *	p(n)	     --> newton step
 *	xpls(n)	    <--	 new iterate x[k]
 *	fpls	    <--	 function value at new iterate, f(xpls)
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in dogdrv (but can be
 *			 modified by fcn).
 *	sx(n)	     --> diagonal scaling matrix for x
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	dlt	    <--> trust region radius
 *			 [retain value between successive calls]
 *	iretcd	    <--	 return code
 *			   =0 satisfactory xpls found
 *			   =1 failed to find satisfactory xpls sufficiently
 *			      distinct from x
 *	mxtake	    <--	 boolean flag indicating step of maximum length used
 *	sc(n)	     --> workspace [current step]
 *	wrk1(n)	     --> workspace (and place holding argument to tregup)
 *	wrk2(n)	     --> workspace
 *	wrk3(n)	     --> workspace
 *	ipr	     --> device to which to send output */

    Rboolean fstdog, nwtake;
    int i;
    double fplsp = 0.0, rnwtln, eta = 0.0, cln = 0.0, tmp; /* -Wall */

    tmp = 0.;
    for (i = 0; i < n; ++i)
	tmp += sx[i] * sx[i] * p[i] * p[i];
    rnwtln = sqrt(tmp);

    *iretcd = 4;
    fstdog = TRUE;
    do {
	/*	find new step by double dogleg algorithm */

	dog_1step(nr, n, g, a, p, sx, rnwtln, dlt, &nwtake,
		  &fstdog, wrk1, wrk2, &cln, &eta, sc, stepmx);

	/*	check new point and update trust region */

	tregup(nr, n, x, f, g, a, (fcn_p)fcn, state, sc, sx, nwtake, stepmx,
	       steptl, dlt, iretcd, wrk3, &fplsp, xpls, fpls, mxtake,
	       2, wrk1);
    } while(*iretcd > 1);
} /* dogdrv */


static void
hook_1step(int nr, int n, double *g, double *a, double *udiag, double *p,
       double *sx, double rnwtln, double *dlt, double *amu, double dltp,
       double *phi, double *phip0, Rboolean *fstime, double *sc,
       Rboolean *nwtake, double *wrk0, double epsm)
{
/* Find new step by more-hebdon algorithm  (iff	 method == 3);
 * repeatedly called by hookdrv() only.

 * PARAMETERS :

 *	 nr	      --> row dimension of matrix
 *	 n	      --> dimension of problem
 *	 g(n)	      --> gradient at current iterate, g(x)
 *	 a(n,n)	      --> cholesky decomposition of hessian in
 *			  lower triangular part and diagonal.
 *			  hessian or approx in upper triangular part
 *	 udiag(n)     --> diagonal of hessian in a(.,.)
 *	 p(n)	      --> newton step
 *	 sx(n)	      --> diagonal scaling matrix for n
 *	 rnwtln	      --> newton step length
 *	 dlt	     <--> trust region radius
 *	 amu	     <--> [retain value between successive calls]
 *	 dltp	      --> trust region radius at last exit from this routine
 *	 phi	     <--> [retain value between successive calls]
 *	 phip0	     <--> [retain value between successive calls]
 *	 fstime	     <--> boolean. =.true. if first entry to this routine
 *			  during k-th iteration
 *	 sc(n)	     <--  current step
 *	 nwtake	     <--  boolean, =.true. if newton step taken
 *	 wrk0(n)      --> workspace
 *	 epsm	      --> machine epsilon
 */
    int one = 1, job = 0, info;
    int i, j;
    double phip;
    double amulo, amuup;
    double addmax, stepln;
    double temp1;
    const double hi = 1.5, alo = 0.75;
    /*	 hi and alo are constants used in this routine. */
    /*	 change here if other values are to be substituted. */

    /*	shall we take newton step ? */
    *nwtake = (rnwtln <= hi * *dlt);

    if (*nwtake) { /*	take newton step */
	for (i = 0; i < n; ++i)
	    sc[i] = p[i];
	*dlt = fmin2(*dlt, rnwtln);
	*amu = 0.;
	return;
    }

    /* else *nwtake = FALSE :	newton step not taken */
    if (*amu > 0.)
	*amu -= (*phi + dltp) * (dltp - *dlt + *phi) / (*dlt * *phip0);

    *phi = rnwtln - *dlt;
    if (*fstime) {
	for (i = 0; i < n; ++i)
	    wrk0[i] = sx[i] * sx[i] * p[i];
	/*	  solve l*y = (sx**2)*p */
	F77_CALL(dtrsl)(a, &nr, &n, wrk0, &job, &info);
	/* Computing 2nd power */
	temp1 = F77_CALL(dnrm2)(&n, wrk0, &one);
	*phip0 = -(temp1 * temp1) / rnwtln;
	*fstime = FALSE;
    }
    phip = *phip0;
    amulo = -(*phi) / phip;
    amuup = 0.;
    for (i = 0; i < n; ++i)
	amuup += g[i] * g[i] / (sx[i] * sx[i]);
    amuup = sqrt(amuup) / *dlt;

    while (1) {
	/*	test value of amu; generate next amu if necessary */
	if (*amu < amulo || *amu > amuup) {
	    *amu = fmax2(sqrt(amulo * amuup), amuup * .001);
	}

	/*	copy (h,udiag) to l */
	/*	where h <-- h+amu*(sx**2) [do not actually change (h,udiag)] */
/* The original code was
   http://people.sc.fsu.edu/~jburkardt/f77_src/uncmin/uncmin.f
	do 100 j=1,n
	  a(j,j)=udiag(j) + amu*sx(j)*sx(j)
	  if(j.eq.n) go to 100
	  jp1=j+1
	  do i=jp1,n
	    a(i,j)=a(j,i)
	  end do
  100   continue
*/
	for (i = 0; i < n; ++i) {
	    a[i + i * nr] = udiag[i] + *amu * sx[i] * sx[i];
	    for (j = 0; j < i; ++j) // changed from ++i 2012-05-31, PR#
		a[i + j * nr] = a[j + i * nr];
	}

	/*	factor h=l(l+) */

	temp1 = sqrt(epsm);
	choldc(nr, n, a, 0.0, temp1, &addmax);

	/*	solve h*p = l(l+)*sc = -g */
	for (i = 0; i < n; ++i)
	    wrk0[i] = -g[i];
	lltslv(nr, n, a, sc, wrk0);

	/* reset h.  note since udiag has not been destroyed we need do */
	/* nothing here.	h is in the upper part and in udiag, still intact */

	stepln = 0.;
	for (i = 0; i < n; ++i)
	    stepln += sx[i] * sx[i] * sc[i] * sc[i];
	stepln = sqrt(stepln);
	*phi = stepln - *dlt;
	for (i = 0; i < n; ++i)
	    wrk0[i] = sx[i] * sx[i] * sc[i];
	F77_CALL(dtrsl)(a, &nr, &n, wrk0, &job, &info);
	temp1 = F77_CALL(dnrm2)(&n, wrk0, &one);
	phip = -(temp1 * temp1) / stepln;
	if ((alo * *dlt <= stepln && stepln <= hi * *dlt)
	    || (amuup - amulo > 0.)) {
	    /*	  sc is acceptable hookstep */
	    break;
	}
	else { /*	  sc not acceptable hookstep.  select new amu */
	    temp1 = (*amu - *phi)/phip;
	    amulo = fmax2(amulo, temp1);
	    if (*phi < 0.)
		amuup = fmin2(amuup,*amu);
	    *amu -= stepln * *phi / (*dlt * phip);
	}
    }
} /* hook_1step */

static void
hookdrv(int nr, int n, double *x, double f, double *g, double *a,
	double *udiag, double *p, double *xpls, double *fpls, fcn_p fcn,
	void *state, double *sx, double stepmx, double steptl, double *dlt,
	int *iretcd, Rboolean *mxtake, double *amu, double *dltp,
	double *phi, double *phip0, double *sc, double *xplsp,
	double *wrk0, double epsm, int itncnt)
{
/* Find a next newton iterate (xpls) by the more-hebdon method.
 * (iff	 method == 3)

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate x[k-1]
 *	f	     --> function value at old iterate, f(x)
 *	g(n)	     --> gradient at old iterate, g(x), or approximate
 *	a(n,n)	     --> cholesky decomposition of hessian in lower
 *			 triangular part and diagonal.
 *			 hessian in upper triangular part and udiag.
 *	udiag(n)     --> diagonal of hessian in a(.,.)
 *	p(n)	     --> newton step
 *	xpls(n)	    <--	 new iterate x[k]
 *	fpls	    <--	 function value at new iterate, f(xpls)
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in hookdrv (but can be
 *			 modified by fcn).
 *	sx(n)	     --> diagonal scaling matrix for x
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	dlt	    <--> trust region radius
 *	iretcd	    <--	 return code
 *			   =0 satisfactory xpls found
 *			   =1 failed to find satisfactory xpls sufficiently
 *			      distinct from x
 *	mxtake	    <--	 boolean flag indicating step of maximum length used
 *	amu	    <--> [retain value between successive calls]
 *	dltp	    <--> [retain value between successive calls]
 *	phi	    <--> [retain value between successive calls]
 *	phip0	    <--> [retain value between successive calls]
 *	sc(n)	     --> workspace
 *	xplsp(n)     --> workspace
 *	wrk0(n)	     --> workspace
 *	epsm	     --> machine epsilon
 *	itncnt	     --> iteration count
 *	ipr	     --> device to which to send output
 */

    Rboolean fstime, nwtake;
    int i, j;
    double bet, alpha, fplsp = 0.0 /* -Wall */, rnwtln, tmp;

    tmp = 0.;
    for (i = 0; i < n; ++i)
	tmp += sx[i] * sx[i] * p[i] * p[i];
    rnwtln = sqrt(tmp);

    if (itncnt == 1) {
	*amu = 0.;

	/*	if first iteration and trust region not provided by user,
		compute initial trust region. */

	if (*dlt == -1.) {

	    alpha = 0.;
	    for (i = 0; i < n; ++i)
		alpha += g[i] * g[i] / (sx[i] * sx[i]);
	    bet = 0.;
	    for (i = 0; i < n; ++i) {
		tmp = 0.;
		for (j = i; j < n; ++j)
		    tmp += a[j + i * nr] * g[j] / (sx[j] * sx[j]);

		bet += tmp * tmp;
	    }
	    *dlt = alpha * sqrt(alpha) / bet;
	    if(*dlt > stepmx) *dlt = stepmx;
	}
    }
    *iretcd = 4;
    fstime = TRUE;
    do {
	/*	find new step by more-hebdon algorithm */

	hook_1step(nr, n, g, a, udiag, p, sx, rnwtln, dlt, amu,
		   *dltp, phi, phip0, &fstime, sc, &nwtake, wrk0, epsm);
	*dltp = *dlt;

	/*	check new point and update trust region */

	tregup(nr, n, x, f, g, a, (fcn_p)fcn, state, sc, sx, nwtake, stepmx,
	       steptl, dlt, iretcd, xplsp, &fplsp, xpls, fpls, mxtake,
	       3, udiag);
    } while(*iretcd > 1);
} /* hookdrv */

static void
secunf(int nr, int n, double *x, double *g, double *a, double *udiag,
       double *xpls, double *gpls, double epsm, int itncnt, double rnf,
       int iagflg, Rboolean *noupdt, double *s, double *y, double *t)
{
/* Update hessian by the bfgs unfactored method	 (only when  method == 3)

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate, x[k-1]
 *	g(n)	     --> gradient or approximate at old iterate
 *	a(n,n)	    <--> on entry: approximate hessian at old iterate
 *			   in upper triangular part (and udiag)
 *			 on exit:  updated approx hessian at new iterate
 *			   in lower triangular part and diagonal
 *			 [lower triangular part of symmetric matrix]
 *	udiag	     --> on entry: diagonal of hessian
 *	xpls(n)	     --> new iterate, x[k]
 *	gpls(n)	     --> gradient or approximate at new iterate
 *	epsm	     --> machine epsilon
 *	itncnt	     --> iteration count
 *	rnf	     --> relative noise in optimization function fcn
 *	iagflg	     --> =1 if analytic gradient supplied, =0 otherwise
 *	noupdt	    <--> boolean: no update yet
 *			 [retain value between successive calls]
 *	s(n)	     --> workspace
 *	y(n)	     --> workspace
 *	t(n)	     --> workspace
 */

    double ynrm2, snorm2;
    int i, j, one = 1;
    Rboolean skpupd;
    double gam, tol, den1, den2;


    /*	copy hessian in upper triangular part and udiag to
	lower triangular part and diagonal */
    for (i = 0; i < n; ++i) {
	a[i + i * nr] = udiag[i];
	for (j = 0; j < i; ++j)
	    a[i + j * nr] = a[j + i * nr];
    }

    *noupdt = (itncnt == 1);

    for (i = 0; i < n; ++i) {
	s[i] = xpls[i] - x[i];
	y[i] = gpls[i] - g[i];
    }
    den1 = F77_CALL(ddot)(&n, s, &one, y, &one);
    snorm2 = F77_CALL(dnrm2)(&n, s, &one);
    ynrm2  = F77_CALL(dnrm2)(&n, y, &one);
    if (den1 < sqrt(epsm) * snorm2 * ynrm2) return;

    mvmlts(nr, n, a, s, t);
    den2 = F77_CALL(ddot)(&n, s, &one, t, &one);
    if (*noupdt) {
	/*	  h <-- [(s+)y/(s+)hs]h */

	gam = den1 / den2;
	den2 *= gam;
	for (j = 0; j < n; ++j) {
	    t[j] *= gam;
	    for (i = j; i < n; ++i)
		a[i + j * nr] *= gam;
	}
	*noupdt = FALSE;
    }
    skpupd = TRUE;

    /*	check update condition on row i */

    for (i = 0; i < n; ++i) {
	tol = rnf * fmax2(fabs(g[i]), fabs(gpls[i]));
	if (iagflg == 0)
	    tol /= sqrt(rnf);
	if (fabs(y[i] - t[i]) >= tol) {
	    skpupd = FALSE;
	    break;
	}
    }
    if (skpupd) return;

    /*	  bfgs update */

    for (j = 0; j < n; ++j) {
	for (i = j; i < n; ++i)
	    a[i + j * nr] += y[i] * y[j] / den1 - t[i] * t[j] / den2;
    }
} /* secunf */

static void
secfac(int nr, int n, double *x, double *g, double *a, double *xpls,
       double *gpls, double epsm, int itncnt, double rnf, int iagflg,
       Rboolean *noupdt, double *s, double *y, double *u, double *w)
{
/* Update hessian by the bfgs factored method  (only when  method == 1 or 2)

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> old iterate, x[k-1]
 *	g(n)	     --> gradient or approximate at old iterate
 *	a(n,n)	    <--> on entry: cholesky decomposition of hessian in
 *			   lower part and diagonal.
 *			 on exit:  updated cholesky decomposition of hessian
 *			   in lower triangular part and diagonal
 *	xpls(n)	     --> new iterate, x[k]
 *	gpls(n)	     --> gradient or approximate at new iterate
 *	epsm	     --> machine epsilon
 *	itncnt	     --> iteration count
 *	rnf	     --> relative noise in optimization function fcn
 *	iagflg	     --> =1 if analytic gradient supplied, =0 itherwise
 *	noupdt	    <--> boolean: no update yet
 *			 [retain value between successive calls]
 *	s(n)	     --> workspace
 *	y(n)	     --> workspace
 *	u(n)	     --> workspace
 *	w(n)	     --> workspace */

    double ynrm2;
    int i, j, one = 1;
    Rboolean skpupd;
    double snorm2, reltol;
    double alp, den1, den2;

    *noupdt = (itncnt == 1);

    for (i = 0; i < n; ++i) {
	s[i] = xpls[i] - x[i];
	y[i] = gpls[i] - g[i];
    }
    den1 = F77_CALL(ddot)(&n, s, &one, y, &one);
    snorm2 = F77_CALL(dnrm2)(&n, s, &one);
    ynrm2  = F77_CALL(dnrm2)(&n, y, &one);

    if (den1 < sqrt(epsm) * snorm2 * ynrm2)
	return;

    mvmltu(nr, n, a, s, u);
    den2 = F77_CALL(ddot)(&n, u, &one, u, &one);

    /*	l <-- sqrt(den1/den2)*l */

    alp = sqrt(den1 / den2);
    if (*noupdt) {
	for (j = 0; j < n; ++j) {
	    u[j] = alp * u[j];
	    for (i = j; i < n; ++i) {
		a[i + j * nr] *= alp;
	    }
	}
	*noupdt = FALSE;
	den2 = den1;
	alp = 1.;
    }
    /*	w = l(l+)s = hs */

    mvmltl(nr, n, a, u, w);
    if (iagflg == 0)
	reltol = sqrt(rnf);
    else
	reltol = rnf;

    skpupd = TRUE;
    for (i = 0; i < n; ++i) {
	skpupd = (fabs(y[i] - w[i]) <
		  reltol * fmax2(fabs(g[i]), fabs(gpls[i])));
	if(!skpupd)
	    break;
    }

    if(skpupd)
	return;

    /*	  w = y-alp*l(l+)s */
    for (i = 0; i < n; ++i)
	w[i] = y[i] - alp * w[i];


    /*	  alp=1/sqrt(den1*den2) */
    alp /= den1;

    /*	  u=(l+)/sqrt(den1*den2) = (l+)s/sqrt((y+)s * (s+)l(l+)s) */

    for (i = 0; i < n; ++i)
	u[i] *= alp;

    /*	  copy l into upper triangular part.  zero l. */

    for (i = 1; i < n; ++i) {
	for (j = 0; j < i; ++j) {
	    a[j + i * nr] = a[i + j * nr];
	    a[i + j * nr] = 0.;
	}
    }

    /*	  find q, (l+) such that  q(l+) = (l+) + u(w+) */

    qrupdt(nr, n, a, u, w);

    /* upper triangular part and diagonal of a[] now contain updated
     * cholesky decomposition of hessian.
     * copy back to lower triangular part.
     */
    for (i = 1; i < n; ++i)
	for (j = 0; j < i; ++j)
	    a[i + j * nr] = a[j + i * nr];
} /* secfac */

static void
chlhsn(int nr, int n, double *a, double epsm, double *sx, double *udiag)
{
/*	find the l(l-transpose) [written LL+] decomposition of the perturbed
 *	model hessian matrix a+mu*i(where mu\0 and i is the identity matrix)
 *	which is safely positive definite.  if a is safely positive definite
 *	upon entry, then mu=0.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	    <--> on entry; "a" is model hessian (only lower
 *			 triangular part and diagonal stored)
 *			 on exit:  a contains l of LL+ decomposition of
 *			 perturbed model hessian in lower triangular
 *			 part and diagonal and contains hessian in upper
 *			 triangular part and udiag
 *	epsm	     --> machine epsilon
 *	sx(n)	     --> diagonal scaling matrix for x
 *	udiag(n)    <--	 on exit: contains diagonal of hessian

 * INTERNAL VARIABLES

 *	tol		 tolerance
 *	diagmn		 minimum element on diagonal of a
 *	diagmx		 maximum element on diagonal of a
 *	offmax		 maximum off-diagonal element of a
 *	offrow		 sum of off-diagonal elements in a row of a
 *	evmin		 minimum eigenvalue of a
 *	evmax		 maximum eigenvalue of a

 * DESCRIPTION

 *	1. if "a" has any negative diagonal elements, then choose mu>0
 *	such that the diagonal of a:=a+mu*i is all positive
 *	with the ratio of its smallest to largest element on the
 *	order of sqrt(epsm).

 *	2. "a" undergoes a perturbed cholesky decomposition which
 *	results in an LL+ decomposition of a+d, where d is a
 *	non-negative diagonal matrix which is implicitly added to
 *	"a" during the decomposition if "a" is not positive definite.
 *	"a" is retained and not changed during this process by
 *	copying l into the upper triangular part of "a" and the
 *	diagonal into udiag.  then the cholesky decomposition routine
 *	is called.  on return, addmax contains maximum element of d.

 *	3. if addmax=0, "a" was positive definite going into step 2
 *	and return is made to calling program.	otherwise,
 *	the minimum number sdd which must be added to the
 *	diagonal of a to make it safely strictly diagonally dominant
 *	is calculated.	since a+addmax*i and a+sdd*i are safely
 *	positive definite, choose mu=fmin2(addmax,sdd) and decompose
 *	a+mu*i to obtain l.
 */

    int i, j;
    double evmin, evmax;
    double addmax, diagmn, diagmx, offmax, offrow, posmax;
    double sdd, amu, tol, tmp;


    /*	scale hessian */
    /*	pre- and post- multiply "a" by inv(sx) */

    for (j = 0; j < n; ++j)
	for (i = j; i < n; ++i)
	    a[i + j * nr] /= sx[i] * sx[j];

    /*	step1
     *	-----
     *	note:  if a different tolerance is desired throughout this
     *	algorithm, change tolerance here: */

    tol = sqrt(epsm);

    diagmx = a[0];
    diagmn = a[0];
    if (n > 1) {
	for (i = 1; i < n; ++i) {
	    tmp = a[i + i * nr];
	    if(diagmn > tmp)	diagmn = tmp;
	    if(diagmx < tmp)	diagmx = tmp;
	}
    }
    posmax = fmax2(diagmx, 0.0);

    if (diagmn <= posmax * tol) {
	amu = tol * (posmax - diagmn) - diagmn;
	if (amu == 0.) {
	    /*	find largest off-diagonal element of a */
	    offmax = 0.;
	    for (i = 1; i < n; ++i) {
		for (j = 0; j < i; ++j)
		    if (offmax < (tmp = fabs(a[i + j * nr]))) offmax = tmp;
	    }

	    if (offmax == 0.)
		amu = 1.;
	    else
		amu = offmax * (tol + 1.);
	}
	/*	a=a + mu*i */
	for (i = 0; i < n; ++i)
	    a[i + i * nr] += amu;

	diagmx += amu;
    }

    /*	copy lower triangular part of "a" to upper triangular part */
    /*	and diagonal of "a" to udiag */
    for (i = 0; i < n; ++i) {
	udiag[i] = a[i + i * nr];
	for (j = 0; j < i; ++j)
	    a[j + i * nr] = a[i + j * nr];
    }
    choldc(nr, n, a, diagmx, tol, &addmax);


    /*	step3

     *	if addmax=0, "a" was positive definite going into step 2,
     *	the LL+ decomposition has been done, and we return.
     *	otherwise, addmax>0.  perturb "a" so that it is safely
     *	diagonally dominant and find LL+ decomposition */

    if (addmax > 0.0) {

	/*	restore original "a" (lower triangular part and diagonal) */

	for (i = 0; i < n; ++i) {
	    a[i + i * nr] = udiag[i];
	    for (j = 0; j < i; ++j)
		a[i + j * nr] = a[j + i * nr];
	}

	/*	find sdd such that a+sdd*i is safely positive definite */
	/*	note:  evmin<0 since a is not positive definite; */

	evmin = 0.;
	evmax = a[0];
	for (i = 0; i < n; ++i) {
	    offrow = 0.;
	    for (j = 0; j < i; ++j)
		offrow += fabs(a[i + j * nr]);
	    for (j = i+1; j < n; ++j)
		offrow += fabs(a[j + i * nr]);
	    tmp = a[i + i * nr] - offrow;
	    if(evmin > tmp) evmin = tmp;
	    tmp = a[i + i * nr] + offrow;
	    if(evmax < tmp) evmax = tmp;
	}
	sdd = tol * (evmax - evmin) - evmin;

	/*	perturb "a" and decompose again */

	amu = fmin2(sdd, addmax);
	for (i = 0; i < n; ++i) {
	    a[i + i * nr] += amu;
	    udiag[i] = a[i + i * nr];
	}

	/*	 "a" now guaranteed safely positive definite */

	choldc(nr, n, a, 0.0, tol, &addmax);
    }
    /*	unscale hessian and cholesky decomposition matrix */

    for (j = 0; j < n; ++j) {
	for (i = j; i < n; ++i)
	    a[i + j * nr] *= sx[i];
	for (i = 0; i < j; ++i)
	    a[i + j * nr] *= sx[i] * sx[j];
	udiag[j] *= sx[j] * sx[j];
    }
} /* chlhsn */

static void
hsnint(int nr, int n, double *a, double *sx, int method)
{
/* Provide initial hessian when using secant updates .

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	a(n,n)	    <--	 initial hessian (lower triangular matrix)
 *	sx(n)	     --> diagonal scaling matrix for x
 *	method	     --> algorithm to use to solve minimization problem
 *			   =1,2 factored secant method used
 *			   =3	unfactored secant method used */
    int i, j;

    for (i = 0; i < n; ++i) {
	if (method == 3)
	    a[i + i * nr] = sx[i] * sx[i];
	else
	    a[i + i * nr] = sx[i];

	for (j = 0; j < i; ++j)
	    a[i + j * nr] = 0.;
    }
} /* hsnint */


static void
fstofd(int nr, int m, int n, double *xpls, fcn_p fcn, void *state,
       const double *fpls, double *a, double *sx, double rnoise,
       double *fhat, int icase)
{
/*	find first order forward finite difference approximation "a" to the
 *	first derivative of the function defined by the subprogram "fname"
 *	evaluated at the new iterate "xpls".

 *	for optimization use this routine to estimate:
 *	1) the first derivative (gradient) of the optimization function "fcn
 *	   analytic user routine has been supplied;
 *	2) the second derivative (hessian) of the optimization function
 *	   if no analytic user routine has been supplied for the hessian but
 *	   one has been supplied for the gradient ("fcn") and if the
 *	   optimization function is inexpensive to evaluate

 *	note

 *	_m=1 (optimization) algorithm estimates the gradient of the function
 *	     (fcn).   fcn(x) # f: r(n)-->r(1)
 *	_m=n (systems) algorithm estimates the jacobian of the function
 *	     fcn(x) # f: r(n)-->r(n).
 *	_m=n (optimization) algorithm estimates the hessian of the optimizatio
 *	     function, where the hessian is the first derivative of "fcn"

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	m	     --> number of rows in a
 *	n	     --> number of columns in a; dimension of problem
 *	xpls(n)	     --> new iterate:  x[k]
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in fstofd (but can be
 *			 modified by fcn).
 *	fpls(m)	     --> _m=1 (optimization) function value at new iterate:
 *			      fcn(xpls)
 *			 _m=n (optimization) value of first derivative
 *			      (gradient) given by user function fcn
 *			 _m=n (systems)	 function value of associated
 *			      minimization function
 *	a(nr,n)	    <--	 finite difference approximation (see note).  only
 *			 lower triangular matrix and diagonal are returned
 *	sx(n)	     --> diagonal scaling matrix for x
 *	rnoise	     --> relative noise in fcn [f(x)]
 *	fhat(m)	     --> workspace
 *	icase	     --> =1 optimization (gradient)
 *			 =2 systems
 *			 =3 optimization (hessian)

 *	internal variables

 *	stepsz - stepsize in the j-th variable direction
 */

    int i, j;
    double xtmpj, stepsz, temp1, temp2;

    /*	find j-th column of a
	each column is derivative of f(fcn) with respect to xpls(j) */

    for (j = 0; j < n; ++j) {
	temp1 = fabs(xpls[j]);
	temp2  = 1.0/sx[j];
	stepsz = sqrt(rnoise) * fmax2(temp1, temp2);
	xtmpj = xpls[j];
	xpls[j] = xtmpj + stepsz;
	(*fcn)(n, xpls, fhat, state);
	xpls[j] = xtmpj;
	for (i = 0; i < m; ++i)
	    a[i + j * nr] = (fhat[i] - fpls[i]) / stepsz;
    }
    if (icase == 3 && n > 1) {
	/*	if computing hessian, a must be symmetric */
	for (i = 1; i < m; ++i)
	    for (j = 0; j < i; ++j)
		a[i + j * nr] = (a[i + j * nr] + a[j + i * nr]) / 2.0;
    }
} /* fstofd */

static void fstocd(int n, double *x, fcn_p fcn, void *state,
		   double *sx, double rnoise, double *g)
{
/* Find central difference approximation g to the first derivative
 * (gradient) of the function defined by fcn at the point x.

 * PARAMETERS :

 *	n	     --> dimension of problem
 *	x	     --> point at which gradient is to be approximated.
 *	fcn	     --> name of subroutine to evaluate function.
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in fstocd (but can be
 *			 modified by fcn).
 *	sx	     --> diagonal scaling matrix for x.
 *	rnoise	     --> relative noise in fcn [f(x)].
 *	g	    <--	 central difference approximation to gradient.
 */
    int i;
    double stepi, fplus, fminus, xtempi, temp1, temp2;

    /*	find i th  stepsize, evaluate two neighbors in direction of i th */
    /*	unit vector, and evaluate i th	component of gradient. */

    for (i = 0; i < n; ++i) {
	xtempi = x[i];
	temp1 = fabs(xtempi);
	temp2 = 1.0/sx[i];
	stepi = pow(rnoise, 1.0/3.0) * fmax2(temp1, temp2);
	x[i] = xtempi + stepi;
	(*fcn)(n, x, &fplus, state);
	x[i] = xtempi - stepi;
	(*fcn)(n, x, &fminus, state);
	x[i] = xtempi;
	g[i] = (fplus - fminus) / (stepi * 2.);
    }
} /* fstocd */

static void sndofd(int nr, int n, double *xpls, fcn_p fcn, void *state,
		   double fpls, double *a, double *sx, double rnoise,
		   double *stepsz, double *anbr)
{
/* Find second order forward finite difference approximation "a"
 * to the second derivative (hessian) of the function defined by the subp
 * "fcn" evaluated at the new iterate "xpls" .

 * For optimization use this routine to estimate
 *	1) the second derivative (hessian) of the optimization function
 *	   if no analytical user function has been supplied for either
 *	   the gradient or the hessian and if the optimization function
 *	   "fcn" is inexpensive to evaluate.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	xpls(n)	     --> new iterate:	x[k]
 *	fcn	     --> name of subroutine to evaluate function
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in sndofd (but can be
 *			 modified by fcn).
 *	fpls	     --> function value at new iterate, f(xpls)
 *	a(n,n)	    <--	 finite difference approximation to hessian
 *			 only lower triangular matrix and diagonal
 *			 are returned
 *	sx(n)	     --> diagonal scaling matrix for x
 *	rnoise	     --> relative noise in fname [f(x)]
 *	stepsz(n)    --> workspace (stepsize in i-th component direction)
 *	anbr(n)	     --> workspace (neighbor in i-th direction)
 */
    double fhat;
    int i, j;
    double xtmpi, xtmpj;

    /*	find i-th stepsize and evaluate neighbor in direction
	of i-th unit vector. */

    for (i = 0; i < n; ++i) {
	xtmpi = xpls[i];
	stepsz[i] = pow(rnoise, 1.0/3.0) * fmax2(fabs(xtmpi), 1./sx[i]);
	xpls[i] = xtmpi + stepsz[i];
	(*fcn)(n, xpls, &anbr[i], state);
	xpls[i] = xtmpi;
    }

    /*	calculate row i of a */

    for (i = 0; i < n; ++i) {
	xtmpi = xpls[i];
	xpls[i] = xtmpi + stepsz[i] * 2.;
	(*fcn)(n, xpls, &fhat, state);
	a[i + i * nr] = ((fpls - anbr[i]) + (fhat - anbr[i])) /
	    (stepsz[i] * stepsz[i]);

	/*	calculate sub-diagonal elements of column */
	if(i == 0) {
	    xpls[i] = xtmpi;
	    continue;
	}
	xpls[i] = xtmpi + stepsz[i];
	for (j = 0; j < i; ++j) {
	    xtmpj = xpls[j];
	    xpls[j] = xtmpj + stepsz[j];
	    (*fcn)(n, xpls, &fhat, state);
	    a[i + j*nr] = ((fpls - anbr[i]) + (fhat - anbr[j])) /
		(stepsz[i]*stepsz[j]);
	    xpls[j] = xtmpj;
	}
	xpls[i] = xtmpi;
    }
} /* sndofd */

static void
grdchk(int n, double *x, fcn_p fcn, void *state, double f, double *g,
       double *typsiz, double *sx, double fscale, double rnf,
       double analtl, double *wrk1, int *msg)
{
/* Check analytic gradient against estimated gradient

 * PARAMETERS :

 *	n	     --> dimension of problem
 *	x(n)	     --> estimate to a root of fcn
 *	fcn	     --> name of subroutine to evaluate optimization function
 *			 must be declared external in calling routine
 *			      fcn:  r(n) --> r(1)
 *	state	    <--> information other than x and n that fcn requires.
 *			 state is not modified in grdchk (but can be
 *			 modified by fcn).
 *	f	     --> function value:  fcn(x)
 *	g(n)	     --> gradient:  g(x)
 *	typsiz(n)    --> typical size for each component of x
 *	sx(n)	     --> diagonal scaling matrix:  sx(i)=1./typsiz(i)
 *	fscale	     --> estimate of scale of objective function fcn
 *	rnf	     --> relative noise in optimization function fcn
 *	analtl	     --> tolerance for comparison of estimated and
 *			 analytical gradients
 *	wrk1(n)	     --> workspace
 *	msg	    <--	 message or error code
 *			   on output: =-21, probable coding error of gradient
 */
    int i;
    double gs, wrk;

    /*	compute first order finite difference gradient
	and compare to analytic gradient. */

    fstofd(1, 1, n, x, (fcn_p)fcn, state, &f, wrk1, sx, rnf, &wrk, 1);
    for (i = 0; i < n; ++i) {
	gs = fmax2(fabs(f), fscale) / fmax2(fabs(x[i]), typsiz[i]);
	if (fabs(g[i] - wrk1[i]) > fmax2(fabs(g[i]), gs) * analtl) {
	    *msg = -21; return;
	}
    }
} /* grdchk */

static void
heschk(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn,
       void *state, double f, double *g, double *a, double *typsiz,
       double *sx, double rnf, double analtl, int iagflg, double *udiag,
       double *wrk1, double *wrk2, int *msg)
{
/* Check analytic hessian against estimated hessian
 *	 (this may be done only if the user supplied analytic hessian
 *	  d2fcn fills only the lower triangular part and diagonal of a)

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> estimate to a root of fcn
 *	fcn	     --> name of subroutine to evaluate optimization function
 *			 must be declared external in calling routine
 *			      fcn:  r(n) --> r(1)
 *	d1fcn	     --> name of subroutine to evaluate gradient of fcn.
 *			 must be declared external in calling routine
 *	d2fcn	     --> name of subroutine to evaluate hessian of fcn.
 *			 must be declared external in calling routine
 *	state	    <--> information other than x and n that fcn,
 *			 d1fcn and d2fcn requires.
 *			 state is not modified in heschk (but can be
 *			 modified by fcn, d1fcn or d2fcn).
 *	f	     --> function value:  fcn(x)
 *	g(n)	    <--	 gradient:  g(x)
 *	a(n,n)	    <--	 on exit:  hessian in lower triangular part and diag
 *	typsiz(n)    --> typical size for each component of x
 *	sx(n)	     --> diagonal scaling matrix:  sx(i)=1./typsiz(i)
 *	rnf	     --> relative noise in optimization function fcn
 *	analtl	     --> tolerance for comparison of estimated and
 *			 analytical gradients
 *	iagflg	     --> =1 if analytic gradient supplied
 *	udiag(n)     --> workspace
 *	wrk1(n)	     --> workspace
 *	wrk2(n)	     --> workspace
 *	msg	    <--> message or error code
 *			   on input : if = 1xx do not compare anal + est hess
 *			   on output: = -22, probable coding error of hessian
 */
    int i, j;
    double hs, temp1, temp2;

    /* compute finite difference approximation a to the hessian. */

    if (iagflg)
	fstofd(nr, n, n, x, (fcn_p)d1fcn, state, g, a, sx, rnf, wrk1, 3);
    else
	sndofd(nr, n, x, (fcn_p)fcn, state, f, a, sx, rnf, wrk1, wrk2);

    /*	copy lower triangular part of "a" to upper triangular part
	and diagonal of "a" to udiag */
    for (j = 0; j < n; ++j) {
	udiag[j] = a[j + j * nr];
	for (i = j+1; i < n; ++i)
	    a[j + i * nr] = a[i + j * nr];
    }

    /* compute analytic hessian and compare to finite difference approximation. */
    (*d2fcn)(nr, n, x, a, state);
    for (j = 0; j < n; ++j) {
	hs = fmax2(fabs(g[j]), 1.0) / fmax2(fabs(x[j]), typsiz[j]);
	if (fabs(a[j + j * nr] - udiag[j]) >
	    fmax2(fabs(udiag[j]), hs) * analtl) {
	    *msg = -22; return;
	}
	for (i = j+1; i < n; ++i) {
	    temp1 = a[i + j * nr];
	    temp2 = fabs(temp1 - a[j + i * nr]);
	    temp1 = fabs(temp1);
	    if (temp2 > fmax2(temp1, hs) * analtl) {
		*msg = -22; return;
	    }
	}
    }
} /* heschk */

static
int opt_stop(int n, double *xpls, double fpls, double *gpls, double *x,
	     int itncnt, int *icscmx, double gradtl, double steptl,
	     double *sx, double fscale, int itnlim,
	     int iretcd, Rboolean mxtake, int *msg)
{
/* Unconstrained minimization stopping criteria :

 * Find whether the algorithm should terminate, due to any
 * of the following:
 *	1) problem solved within user tolerance
 *	2) convergence within user tolerance
 *	3) iteration limit reached
 *	4) divergence or too restrictive maximum step (stepmx) suspected

 * ARGUMENTS :

 *	n	     --> dimension of problem
 *	xpls(n)	     --> new iterate x[k]
 *	fpls	     --> function value at new iterate f(xpls)
 *	gpls(n)	     --> gradient at new iterate, g(xpls), or approximate
 *	x(n)	     --> old iterate x[k-1]
 *	itncnt	     --> current iteration k
 *	icscmx	    <--> number consecutive steps >= stepmx
 *			 [retain value between successive calls]
 *	gradtl	     --> tolerance at which relative gradient considered close
 *			 enough to zero to terminate algorithm
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	sx(n)	     --> diagonal scaling matrix for x
 *	fscale	     --> estimate of scale of objective function
 *	itnlim	     --> maximum number of allowable iterations
 *	iretcd	     --> return code
 *	mxtake	     --> boolean flag indicating step of maximum length used
 *	msg	     --> if msg includes a term 8, suppress output
 *
 * VALUE :
 *	`itrmcd' : termination code
 */

    int i, jtrmcd;
    double d, relgrd, relstp, rgx, rsx;

    /*	last global step failed to locate a point lower than x */
    if (iretcd == 1)
	return 3;

    /* else : */

    /* find direction in which relative gradient maximum. */

    /* check whether within tolerance */
    d = fmax2(fabs(fpls), fscale);
    rgx = 0.;
    for (i = 0; i < n; ++i) {
	relgrd = fabs(gpls[i]) * fmax2(fabs(xpls[i]), 1./sx[i]) / d;
	if(rgx < relgrd) rgx = relgrd;
    }
    jtrmcd = 1;
    if (rgx > gradtl) {

	if (itncnt == 0)
	    return 0;

	/* find direction in which relative stepsize maximum */

	/* check whether within tolerance. */
	rsx = 0.;
	for (i = 0; i < n; ++i) {
	    relstp = fabs(xpls[i] - x[i]) / fmax2(fabs(xpls[i]), 1./sx[i]);
	    if(rsx < relstp) rsx = relstp;
	}
	jtrmcd = 2;
	if (rsx > steptl) { /*	check iteration limit */
	    jtrmcd = 4;
	    if (itncnt < itnlim) {
		/*	check number of consecutive steps \ stepmx */
		if (!mxtake) {
		    *icscmx = 0; return 0;
		} else {
		    ++(*icscmx);
		    if (*icscmx < 5) return 0;
		    jtrmcd = 5;
		}
	    }
	}
    }
    return jtrmcd;
} /* opt_stop */

static void
optchk(int n, double *x, double *typsiz, double *sx, double *fscale,
       double gradtl, int *itnlim, int *ndigit, double epsm, double *dlt,
       int *method, int *iexp, int *iagflg, int *iahflg, double *stepmx,
       int *msg)
{
/* Check input for reasonableness.
 * Return *msg in {-1,-2,..,-7}	 if something is wrong

 * PARAMETERS :

 *	n	     --> dimension of problem
 *	x(n)	     --> on entry, estimate to root of fcn
 *	typsiz(n)   <--> typical size of each component of x
 *	sx(n)	    <--	 diagonal scaling matrix for x
 *	fscale	    <--> estimate of scale of objective function fcn
 *	gradtl	     --> tolerance at which gradient considered close
 *			 enough to zero to terminate algorithm
 *	itnlim	    <--> maximum number of allowable iterations
 *	ndigit	    <--> number of good digits in optimization function fcn
 *	epsm	     --> machine epsilon
 *	dlt	    <--> trust region radius
 *	method	    <--> algorithm indicator
 *	iexp	    <--> expense flag
 *	iagflg	    <--> =1 if analytic gradient supplied
 *	iahflg	    <--> =1 if analytic hessian supplied
 *	stepmx	    <--> maximum step size
 *	msg	    <--> message and error code
 *	ipr	     --> device to which to send output
 */
    int i;
    double stpsiz;

    /*	check that parameters only take on acceptable values.
	if not, set them to default values. */

    if (*method < 1 || *method > 3) *method = 1;
    if (*iagflg != 1)		*iagflg = 0;
    if (*iahflg != 1)		*iahflg = 0;
    if (*iexp != 0)		*iexp = 1;
    if (*msg / 2 % 2 == 1 && *iagflg == 0) {
	*msg = -6;  return;/* 830 write(ipr,906) msg,iagflg */
    }
    if (*msg / 4 % 2 == 1 && *iahflg == 0) {
	*msg = -7;  return;/* 835 write(ipr,907) msg,iahflg */
    }

    /*	check dimension of problem */

    if (n <= 0) {
	*msg = -1;  return;/* 805 write(ipr,901) n */
    }
    if (n == 1 && *msg % 2 == 0) {
	*msg = -2;  return;/* 810 write(ipr,902) */
    }

    /*	compute scale matrix */
    for (i = 0; i < n; ++i) {
	if (typsiz[i] == 0.)
	    typsiz[i] = 1.;
	else if (typsiz[i] < 0.)
	    typsiz[i] = -typsiz[i];
	sx[i] = 1. / typsiz[i];
    }

    /*	compute default maximum step size if not provided */
    if (*stepmx <= 0.) {
	stpsiz = 0.;
	for (i = 0; i < n; ++i)
	    stpsiz += x[i] * x[i] * sx[i] * sx[i];
	*stepmx = 1000. * fmax2(sqrt(stpsiz), 1.);
    }

    /*	check function scale */
    if (*fscale == 0.)
	*fscale = 1.;
    else if (*fscale < 0.)
	*fscale = -(*fscale);

    /*	check gradient tolerance */
    if (gradtl < 0.) {
	*msg = -3;  return;/* 815 write(ipr,903) gradtl */
    }

    /*	check iteration limit */
    if (*itnlim <= 0) {
	*msg = -4;  return;/* 820 write(ipr,904) itnlim */
    }

    /*	check number of digits of accuracy in function fcn */
    if (*ndigit == 0) {
	*msg = -5;  return;/* 825 write(ipr,905) ndigit */
    }
    else if (*ndigit < 0) /* use default, = 15 for IEEE double */
	*ndigit = (int) (-log10(epsm));

    /*	check trust region radius */
    if (*dlt <= 0.) {
	*dlt = -1.;
    } else if (*dlt > *stepmx) {
	*dlt = *stepmx;
    }
    return;
} /* optchk */

static void
prt_result(int nr, int n, const double x[], double f, const double g[],
	   const double *a, const double p[], int itncnt, int iflg)
{
/*
 *  PURPOSE
 *
 *  Print information on current iteration.
 *
 *  PARAMETERS
 *
 *  nr	   --> row dimension of matrix
 *  n	   --> dimension of problem
 *  x(n)   --> iterate x[k]
 *  f	   --> function value at x[k]
 *  g(n)   --> gradient at x[k]
 *  a(n,n) --> hessian at x[k]
 *  p(n)   --> step taken
 *  itncnt --> iteration number k
 *  iflg   --> flag controlling info to print
 */

    /* Print iteration number */

    Rprintf("iteration = %d\n", itncnt);

    /* Print step */

    if (iflg != 0) {
	Rprintf("Step:\n");
	printRealVector((double *)p, n, 1);
    }

    /* Print current iterate */

    Rprintf("Parameter:\n");
    printRealVector((double *)x, n, 1);

    /* Print function value */

    Rprintf("Function Value\n");
    printRealVector((double *)&f, 1, 1);

    /* Print gradient */

    Rprintf("Gradient:\n");
    printRealVector((double *)g, n, 1);

#ifdef NEVER
    /*	Print Hessian */

    /*	We don't do this because the printRealMatrix
	code takes a SEXP rather than a double*.
	We could do something ugly like use fixed e format
	but that would be UGLY!
     */
    if (iflg != 0) {
    }
#endif

    Rprintf("\n");
} /* prt_result */

static void
optdrv_end(int nr, int n, double *xpls, double *x, double *gpls,
	   double *g, double *fpls, double f, double *a, double *p,
	   int itncnt, int itrmcd, int *msg,
	   void (*print_result)(int, int, const double *, double,
				const double *, const double *,
				const double *, int, int))
{
    int i;

    /*	termination :
	reset xpls,fpls,gpls,  if previous iterate solution */
    if (itrmcd == 3) {
	*fpls = f;
	for (i = 0; i < n; ++i) {
	    xpls[i] = x[i];
	    gpls[i] = g[i];
	}
    }
    if (*msg / 8 % 2 == 0)
	(*print_result)(nr, n, xpls, *fpls, gpls, a, p, itncnt, 0);

    *msg = 0;
} /* optdrv_end */

static void
optdrv(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn,
       void *state, double *typsiz, double fscale, int method,
       int iexp, int *msg, int ndigit, int itnlim, int iagflg, int iahflg,
       double dlt, double gradtl, double stepmx, double steptl,
       double *xpls, double *fpls, double *gpls, int *itrmcd,
       double *a, double *udiag, double *g, double *p, double *sx,
       double *wrk0, double *wrk1, double *wrk2, double *wrk3, int *itncnt)
{
/* Driver for non-linear optimization problem  -- called by optif0() & optif9()

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> on entry: estimate to a root of fcn
 *	fcn	     --> name of subroutine to evaluate optimization function
 *			 must be declared external in calling routine
 *				   fcn: R^n --> R
 *	d1fcn	     --> (optional) name of subroutine to evaluate gradient
 *			 of fcn.  must be declared external in calling routine
 *	d2fcn	     --> (optional) name of subroutine to evaluate
 *			 hessian of of fcn.  must be declared external
 *			 in calling routine
 *	state	    <--> information other than x and n that fcn,
 *			 d1fcn and d2fcn requires.
 *			 state is not modified in optdrv (but can be
 *			 modified by fcn, d1fcn or d2fcn).
 *	typsiz(n)    --> typical size for each component of x
 *	fscale	     --> estimate of scale of objective function
 *	method	     --> algorithm to use to solve minimization problem
 *			   =1 line search
 *			   =2 double dogleg
 *			   =3 more-hebdon
 *	iexp	     --> =1 if optimization function fcn is expensive to
 *			 evaluate, =0 otherwise.  if set then hessian will
 *			 be evaluated by secant update instead of
 *			 analytically or by finite differences
 *	msg	    <--> on input:  ( > 0) message to inhibit certain
 *			   automatic checks; see do_nlm() in optimize.c
 *			 on output: ( < 0) error code; =0 no error
 *	ndigit	     --> number of good digits in optimization function fcn
 *	itnlim	     --> maximum number of allowable iterations
 *	iagflg	     --> =1 if analytic gradient supplied
 *	iahflg	     --> =1 if analytic hessian supplied
 *	dlt	     --> trust region radius
 *	gradtl	     --> tolerance at which gradient considered close
 *			 enough to zero to terminate algorithm
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	xpls(n)	    <--> on exit:  xpls is local minimum
 *	fpls	    <--> on exit:  function value at solution, xpls
 *	gpls(n)	    <--> on exit:  gradient at solution xpls
 *	itrmcd	    <--	 termination code
 *	a(n,n)	     --> workspace for hessian (or estimate)
 *			 and its cholesky decomposition
 *	udiag(n)     --> workspace [for diagonal of hessian]
 *	g(n)	     --> workspace (for gradient at current iterate)
 *	p(n)	     --> workspace for step
 *	sx(n)	     --> workspace (for diagonal scaling matrix)
 *	wrk0(n)	     --> workspace
 *	wrk1(n)	     --> workspace
 *	wrk2(n)	     --> workspace
 *	wrk3(n)	     --> workspace
 *	itncnt	     current iteration, k  {{was `internal'}}


 *	internal variables

 *	analtl		 tolerance for comparison of estimated and
 *			 analytical gradients and hessians
 *	epsm		 machine epsilon
 *	f		 function value: fcn(x)
 *	rnf		 relative noise in optimization function fcn.
 *			      noise=10.**(-ndigit)
 */
    Rboolean mxtake = FALSE, noupdt;
    int i, iretcd = 0, icscmx = 0;
    double dltp = 0., epsm, phip0 = 0., f, analtl;
    double dlpsav = 0., phisav = 0., dltsav = 0.;/* -Wall */
    double amusav = 0., phpsav = 0.;		/* -Wall */
    double phi = 0., amu = 0., rnf, wrk;

    *itncnt = 0;
    epsm = DBL_EPSILON;
    optchk(n, x, typsiz, sx, &fscale, gradtl, &itnlim, &ndigit, epsm,
	   &dlt, &method, &iexp, &iagflg, &iahflg, &stepmx, msg);
    if (*msg < 0) return;

    for (i = 0; i < n; ++i)
	p[i] = 0.;

    rnf = Rexp10(-(double)ndigit);
    rnf = fmax2(rnf, epsm);
    analtl = sqrt(rnf);
    analtl = fmax2(0.1, analtl);

    /*	evaluate fcn(x) */

    (*fcn)(n, x, &f, state);

    /*	evaluate analytic or finite difference gradient and
	check analytic gradient, if requested. */

    if (!iagflg) {
	fstofd(1, 1, n, x, (fcn_p)fcn, state, &f, g, sx, rnf, &wrk, 1);
    }
    else { /* analytic gradient */
	(*d1fcn)(n, x, g, state);
	if (*msg / 2 % 2 == 0) {
	    grdchk(n, x, (fcn_p)fcn, state, f, g, typsiz, sx, fscale, rnf,
		   analtl, wrk1, msg);
	    if (*msg < 0) return;
	}
    }
    iretcd = -1;
    *itrmcd = opt_stop(n, x, f, g, wrk1, *itncnt, &icscmx,
		       gradtl, steptl, sx, fscale, itnlim, iretcd,
		       /* mxtake = */FALSE, msg);
    if (*itrmcd != 0) {
	optdrv_end(nr, n, xpls, x, gpls, g, fpls, f, a, p, *itncnt,
		   3, msg, prt_result);
	return;
    }

    if (iexp) {
	/* if optimization function expensive to evaluate (iexp=1), then
	 * hessian will be obtained by secant updates.	get initial hessian.*/
	hsnint(nr, n, a, sx, method);
    }
    else {
	/* evaluate analytic or finite difference hessian and check analytic
	 * hessian if requested (only if user-supplied analytic hessian
	 * routine d2fcn fills only lower triangular part and diagonal of a).
	 */
	if (!iahflg) { /* no analytic hessian */
	    if (iagflg) /* anal.gradient */
		fstofd(nr, n, n, x, (fcn_p)d1fcn, state, g, a, sx, rnf, wrk1,3);
	    else
		sndofd(nr, n, x, (fcn_p)fcn, state, f, a, sx, rnf, wrk1, wrk2);
	}
	else { /* analytic hessian */
	    if (*msg / 4 % 2 == 1) {
		(*d2fcn)(nr, n, x, a, state);
	    }
	    else {
		heschk(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn, (d2fcn_p)d2fcn,
		       state, f, g, a, typsiz, sx, rnf, analtl, iagflg,
		       udiag, wrk1, wrk2, msg);
		/* heschk evaluates d2fcn and checks it against the finite
		 * difference hessian which it calculates by calling fstofd
		 * (if iagflg == 1) or sndofd (otherwise).
		 */
		if (*msg < 0) return;
	    }
	}
    }
    if (*msg / 8 % 2 == 0)
	prt_result(nr, n, x, f, g, a, p, *itncnt, 1);


    /* THE Iterations : */

    while(1) {

	++(*itncnt);

	/* find perturbed local model hessian and its LL+ decomposition
	 * ( skip this step if line search or dogstep techniques being used
	 *   with secant updates (i.e. method == 1 or 2).
	 *   cholesky decomposition L already obtained from secfac.)
	 */
	if (iexp && method != 3) {
	    goto L105;
	}
    L103:
	chlhsn(nr, n, a, epsm, sx, udiag);

    L105:
	/* solve for newton step: ap=-g */
	for (i = 0; i < n; ++i) wrk1[i] = - g[i];
	lltslv(nr, n, a, p, wrk1);

	/*	decide whether to accept newton step  xpls=x + p */
	/*	or to choose xpls by a global strategy. */

	if (iagflg == 0 && method != 1) {
	    dltsav = dlt;
	    if (method != 2) {/* i.e. method = 3 */
		amusav = amu;
		dlpsav = dltp;
		phisav = phi;
		phpsav = phip0;
	    }
	}
	switch(method) {
	case 1:
	    lnsrch(n, x, f, g, p, xpls, fpls, (fcn_p)fcn, state, &mxtake,
		   &iretcd, stepmx, steptl, sx);
	    break;
	case 2:
	    dogdrv(nr, n, x, f, g, a, p, xpls, fpls, (fcn_p)fcn, state,
		   sx, stepmx, steptl, &dlt, &iretcd, &mxtake, wrk0, wrk1,
		   wrk2, wrk3, itncnt);
	    break;
	case 3:
	    hookdrv(nr, n, x, f, g, a, udiag, p, xpls, fpls, (fcn_p)fcn,
		    state, sx, stepmx, steptl, &dlt, &iretcd, &mxtake, &amu,
		    &dltp, &phi, &phip0, wrk0, wrk1 , wrk2, epsm, *itncnt);
	    break;
	}

	/*	if could not find satisfactory step and forward difference */
	/*	gradient was used, retry using central difference gradient. */

	if (iretcd == 1 && iagflg == 0) {
	    iagflg = -1; /*	 set iagflg for central differences */
	    fstocd(n, x, (fcn_p)fcn, state, sx, rnf, g);
	    if (method == 1)	goto L105;
	    dlt = dltsav;
	    if (method == 2)	goto L105;
	    /* else : method == 3 */
	    amu = amusav;
	    dltp = dlpsav;
	    phi = phisav;
	    phip0 = phpsav;
				goto L103;
	}
	/*	calculate step for output */
	for (i = 0; i < n; ++i)
	    p[i] = xpls[i] - x[i];

	/*	calculate gradient at xpls */
	switch(iagflg) {
	case -1:
	    /* central difference gradient */
	    fstocd(n, xpls, (fcn_p)fcn, state, sx, rnf, gpls);
	    break;
	case 0:
	    /* forward difference gradient */
	    fstofd(1, 1, n, xpls, (fcn_p)fcn, state, fpls, gpls, sx, rnf,
		   &wrk, 1);
	    break;
	default:
	    /* analytic gradient */
	    (*d1fcn)(n, xpls, gpls, state);
	}

	/*	check whether stopping criteria satisfied */
	*itrmcd = opt_stop(n, xpls, *fpls, gpls, x, *itncnt, &icscmx,
			   gradtl, steptl, sx, fscale, itnlim, iretcd, mxtake,
			   msg);
	if(*itrmcd != 0) break;

	/*	evaluate hessian at xpls */
	if (iexp) { /* expensive obj.fun. */
	    if (method == 3)
		secunf(nr, n, x, g, a, udiag, xpls, gpls, epsm, *itncnt, rnf,
		       iagflg, &noupdt, wrk1, wrk2, wrk3);
	    else
		secfac(nr, n, x, g, a, xpls, gpls, epsm, *itncnt, rnf, iagflg,
		       &noupdt, wrk0, wrk1, wrk2, wrk3);
	}
	else { /* iexp == 0 */
	    if (!iahflg) {
		if (iagflg)
		    fstofd(nr, n, n, xpls, (fcn_p)d1fcn, state, gpls, a, sx,
			   rnf, wrk1, 3);
		else /* (iagflg != 1) */
		    sndofd(nr, n, xpls, (fcn_p)fcn, state, *fpls, a, sx, rnf,
			   wrk1, wrk2);
	    }
	    else /* analytic hessian */
		(*d2fcn)(nr, n, xpls, a, state);
	}
	if (*msg / 16 % 2 == 1)
	    prt_result(nr, n, xpls, *fpls, gpls, a, p, *itncnt, 1);

	/*	x <-- xpls  and	 g <-- gpls  and  f <-- fpls */
	f = *fpls;
	for (i = 0; i < n; ++i) {
	    x[i] = xpls[i];
	    g[i] = gpls[i];
	}
    } /* END while(1) */

    optdrv_end(nr, n, xpls, x, gpls, g, fpls, f, a, p, *itncnt,
	       *itrmcd, msg, prt_result);
} /* optdrv */

#if 0
static void
dfault(int n, double *x,
       double *typsiz, double *fscale,
       int *method, int *iexp, int *msg,
       int *ndigit, int *itnlim, int *iagflg, int *iahflg,
       double *dlt, double *gradtl, double *stepmx, double *steptl)
{
/* Set default values for each input variable to minimization algorithm
 * for optif0() only.

 * PARAMETERS :

 * INPUT:
 *	n	  dimension of problem
 *	x(n)	  initial guess to solution (to compute max step size)
 * OUTPUT:
 *	typsiz(n) typical size for each component of x
 *	fscale	  estimate of scale of minimization function
 *	method	  algorithm to use to solve minimization problem
 *	iexp	  =0 if minimization function not expensive to evaluate
 *	msg	  message to inhibit certain automatic checks + output
 *	ndigit	  number of good digits in minimization function
 *	itnlim	  maximum number of allowable iterations
 *	iagflg	  =0 if analytic gradient not supplied
 *	iahflg	  =0 if analytic hessian not supplied
 *	dlt	  trust region radius
 *	gradtl	  tolerance at which gradient considered close enough
 *		  to zero to terminate algorithm
 *	stepmx	  value of zero to trip default maximum in optchk
 *	steptl	  tolerance at which successive iterates considered
 *		  close enough to terminate algorithm
 */
  double epsm;
  int i;

  /* set typical size of x and minimization function */
  for (i = 0; i < n; ++i)
    typsiz[i] = 1.;

  *fscale = 1.;

  /* set tolerances */

  epsm = DBL_EPSILON;		/* for IEEE : = 2^-52	  ~= 2.22  e-16 */
  *gradtl = pow(epsm, 1./3.);	/* for IEEE : = 2^(-52/3) ~= 6.055 e-6 */
  *steptl = sqrt(epsm);		/* for IEEE : = 2^-26	  ~= 1.490 e-8 */
  *stepmx = 0.;/* -> compute default in optchk() */
  *dlt = -1.;/* (not needed for method 1) */

  /* set flags */

  *method = 1;
  *iexp	  = 1;
  *msg	  = 0;
  *ndigit = -1;/* -> compute default = floor(-log10(EPS)) in optchk() */
  *itnlim = 150;
  *iagflg = 0;/* no gradient */
  *iahflg = 0;/* no hessian */
} /* dfault() */

void
optif0(int nr, int n, double *x, fcn_p fcn, void *state,
       double *xpls, double *fpls, double *gpls, int *itrmcd,
       double *a, double *wrk)
{
/* Provide simplest interface to minimization package.
 * User has no control over options.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> initial estimate of minimum
 *	fcn	     --> name of routine to evaluate minimization function.
 *			 must be declared external in calling routine.
 *	state	    <--> information other than x and n that fcn requires
 *			 state is not modified in optif0 (but can be
 *			 modified by fcn).
 *	xpls(n)	    <--	 local minimum
 *	fpls	    <--	 function value at local minimum xpls
 *	gpls(n)	    <--	 gradient at local minimum xpls
 *	itrmcd	    <--	 termination code
 *	a(n,n)	     --> workspace
 *	wrk(n,9)     --> workspace
 */
  int iexp, iagflg, iahflg;
  int ndigit, method, itnlim, itncnt;
  double fscale, gradtl, steptl, stepmx, dlt;
  int msg;

  /* Function Body */
  dfault(n, x, &wrk[nr], &fscale, &method, &iexp, &msg, &ndigit,
	 &itnlim, &iagflg, &iahflg, &dlt, &gradtl, &stepmx, &steptl);
  optdrv(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn_dum, (d2fcn_p)d2fcn_dum,
	 state, &wrk[nr * 3], fscale, method, iexp, &msg, ndigit,
	 itnlim, iagflg, iahflg, dlt, gradtl, stepmx, steptl,
	 xpls, fpls, gpls, itrmcd, a, wrk, &wrk[nr], &wrk[nr * 2],
	 &wrk[nr * 4], &wrk[nr * 5], &wrk[nr * 6], &wrk[nr * 7],
	 &wrk[nr * 8], &itncnt);
} /* optif0 */
#endif

/* ---- this one is called from optimize.c : --------------- */
void
optif9(int nr, int n, double *x, fcn_p fcn, fcn_p d1fcn, d2fcn_p d2fcn,
       void *state, double *typsiz, double fscale, int method,
       int iexp, int *msg, int ndigit, int itnlim, int iagflg, int iahflg,
       double dlt, double gradtl, double stepmx, double steptl,
       double *xpls, double *fpls, double *gpls, int *itrmcd, double *a,
       double *wrk, int *itncnt)
{
/*	provide complete interface to minimization package.
 *	user has full control over options.

 * PARAMETERS :

 *	nr	     --> row dimension of matrix
 *	n	     --> dimension of problem
 *	x(n)	     --> on entry: estimate to a root of fcn
 *	fcn	     --> name of subroutine to evaluate optimization function
 *			 must be declared external in calling routine
 *				   fcn: r(n) --> r(1)
 *	d1fcn	     --> (optional) name of subroutine to evaluate gradient
 *			 of fcn.  must be declared external in calling routine
 *	d2fcn	     --> (optional) name of subroutine to evaluate hessian of
 *			 of fcn.  must be declared external in calling routine
 *	state	    <--> information other than x and n that fcn,
 *			 d1fcn and d2fcn requires.
 *			 state is not modified in optif9 (but can be
 *			 modified by fcn, d1fcn or d2fcn).
 *	typsiz(n)    --> typical size for each component of x
 *	fscale	     --> estimate of scale of objective function
 *	method	     --> algorithm to use to solve minimization problem
 *			   =1 line search
 *			   =2 double dogleg
 *			   =3 more-hebdon
 *	iexp	     --> =1 if optimization function fcn is expensive to
 *			 evaluate, =0 otherwise.  if set then hessian will
 *			 be evaluated by secant update instead of
 *			 analytically or by finite differences
 *	msg	    <--> on input:  ( > 0) to inhibit certain automatic checks
 *			 on output: ( < 0) error code; =0 no error
 *	ndigit	     --> number of good digits in optimization function fcn
 *	itnlim	     --> maximum number of allowable iterations
 *	iagflg	     --> =1 if analytic gradient supplied
 *	iahflg	     --> =1 if analytic hessian supplied
 *	dlt	     --> trust region radius
 *	gradtl	     --> tolerance at which gradient considered close
 *			 enough to zero to terminate algorithm
 *	stepmx	     --> maximum allowable step size
 *	steptl	     --> relative step size at which successive iterates
 *			 considered close enough to terminate algorithm
 *	xpls(n)	    <--> on exit:  xpls is local minimum
 *	fpls	    <--> on exit:  function value at solution, xpls
 *	gpls(n)	    <--> on exit:  gradient at solution xpls
 *	itrmcd	    <--	 termination code (in 0..5 ; 0 is "perfect");
 *			see optcode() in optimize.c for meaning
 *	a(n,n)	     --> workspace for hessian (or estimate)
 *			 and its cholesky decomposition
 *	wrk(n,8)     --> workspace
 *	itncnt	    <--> iteration count
 */
  optdrv(nr, n, x, (fcn_p)fcn, (fcn_p)d1fcn, (d2fcn_p)d2fcn, state,
	 typsiz, fscale, method, iexp, msg, ndigit, itnlim, iagflg,
	 iahflg, dlt, gradtl, stepmx, steptl, xpls, fpls, gpls,
	 itrmcd, a, wrk, wrk + nr,
	 wrk + nr * 2, wrk + nr * 3, wrk + nr * 4,
	 wrk + nr * 5, wrk + nr * 6, wrk + nr * 7, itncnt);
} /* optif9 */