/* PGF got this code from Dave L Swofford. Thanks Dave! He got it
from a book, and re-coded the fortran into c. This is not GPL'd; I
assume it is copyrighted by DLS. PGF added the newBrent() and
freeBrent() functions. And a little more tweaking.
It is OK with Dave that this code be used as long as acknowlegement
is given to him.
-PGF
*/
#include "brent.h"
brent *newBrent(int dim)
{
brent *aBrent;
aBrent = malloc(sizeof(brent));
aBrent->vv = psdmatrix(dim);
if ((aBrent->d = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("praxis: unable to malloc aBrent->d\n");
exit(1);
}
if ((aBrent->q0 = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("praxis: unable to malloc aBrent->q0\n");
exit(1);
}
if ((aBrent->q1 = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("praxis: unable to malloc aBrent->q1\n");
exit(1);
}
if ((aBrent->xnew = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("praxis: unable to malloc aBrent->xnew\n");
exit(1);
}
if ((aBrent->y = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("praxis: unable to malloc aBrent->y\n");
exit(1);
}
if ((aBrent->z = (double *)malloc(dim * sizeof(double))) == NULL) {
printf("unable to malloc aBrent->z\n");
exit(1);
}
return aBrent;
}
void freeBrent(brent *aBrent)
{
if(aBrent->vv) free_psdmatrix(aBrent->vv);
if(aBrent->d) free(aBrent->d);
if(aBrent->q0) free(aBrent->q0);
if(aBrent->q1) free(aBrent->q1);
if(aBrent->xnew) free(aBrent->xnew);
if(aBrent->y) free(aBrent->y);
if(aBrent->z) free(aBrent->z);
free(aBrent);
}
double praxis(brent *aBrent, double tol, double h, int dim, double *xx, double (*theFunction)(double []))
{
int i, j, k, k2, illc, kl, kt, ktm;
double scbd, ldfac;
double sf, df, f1, lds, t2, sl, dn, s, sz;
aBrent->eps2 = DBL_EPSILON*DBL_EPSILON;
aBrent->m2 = sqrt(DBL_EPSILON);
aBrent->m4 = sqrt(aBrent->m2);
// machine-dependent initializations
aBrent->vsmall = aBrent->eps2*aBrent->eps2;
aBrent->large = 1.0/aBrent->eps2;
aBrent->vlarge = 1.0/aBrent->vsmall;
aBrent->n = dim;
aBrent->fxn = theFunction;
aBrent->gx = xx;
aBrent->toler = tol;
aBrent->htol = h;
// heuristic numbers:
// - if axes may be badly scaled (which should be avoided if possible),
// set scbd=10, otherwise 1
// - if the problem is known to be ill-conditioned, set illc=TRUE,
// otherwise FALSE
// - ktm+1 is the number of iterations without improvement before
// the algorithm terminates
// (see section 7.6 of Brent's book). ktm=4 is very cautious;
// usually ktm=1 is satisfactory.
scbd = 1.0;
illc = FALSE;
ktm = 1;
ldfac = illc ? 0.1 : 0.01;
kt = aBrent->nl = 0;
aBrent->qf1 = aBrent->gfx = (*aBrent->fxn)(aBrent->gx);
aBrent->toler = t2 = aBrent->eps2 + fabs(aBrent->toler);
aBrent->dmin = aBrent->eps2;
if (aBrent->htol < 100.0*aBrent->toler)
aBrent->htol = 100.0*aBrent->toler;
aBrent->ldt = aBrent->htol;
for(i = 0; i < aBrent->n; i++) {
for(j = 0; j < aBrent->n; j++) {
aBrent->vv[i][j] = 0.0;
if(i == j) aBrent->vv[i][j] = 1.0;
}
}
aBrent->d[0] = aBrent->qd0 = 0.0;
for (i = 0; i < aBrent->n; i++) {
aBrent->q0[i] = 0.0; // DLS: this wasn't included in Brent's code,
// but it's important
aBrent->q1[i] = aBrent->gx[i];
}
// ------ main loop ------ //
for(FOREVER)
{
//printf("a aBrent->nl = %i, aBrent->gfx = %f, aBrent->gx[0] = %f, aBrent->gx[1] = %f\n", aBrent->nl, aBrent->gfx, aBrent->gx[0], aBrent->gx[1]);
//printf("a aBrent->nl = %i, aBrent->gx[0] = %f\n", aBrent->nl, aBrent->gx[0]);
sf = aBrent->d[0];
aBrent->d[0] = s = 0.0;
// minimize along first direction
lineMin(aBrent, 0, 2, &aBrent->d[0], &s, aBrent->gfx, FALSE);
if (s < 0.0)
{
for (i = 0; i < aBrent->n; i++)
aBrent->vv[i][0] = -aBrent->vv[i][0];
}
if ((sf <= 0.9*aBrent->d[0]) || (0.9*sf >= aBrent->d[0]))
{
for (i = 1; i < aBrent->n; i++)
aBrent->d[i] = 0.0;
}
for (k = 1; k < aBrent->n; k++)
{
for (i = 0; i < aBrent->n; i++)
aBrent->y[i] = aBrent->gx[i];
sf = aBrent->gfx;
illc = illc || (kt > 0);
for (FOREVER)
{
kl = k;
df = 0.0;
if (illc)
{
// random step to get off resolution valley
for (i = 0; i < aBrent->n; i++)
{
s = aBrent->z[i] = (0.1*aBrent->ldt + t2*pow(10.0, kt))*(ranDoubleUpToOne() - 0.5);
for (j = 0; j < aBrent->n; j++)
aBrent->gx[j] += s*aBrent->vv[j][i];
// peter was here, Aug 24_00, these next few lines, try to avoid negs
//for (j = 0; j < aBrent->n; j++) {
// if(aBrent->gx[j] < 0.0) {
// aBrent->gx[j] = 1.0e-12 * ranDoubleUpToOne();
// printf("Brent deneggify 1 here.\n");
// }
//}
}
aBrent->gfx = (*aBrent->fxn)(aBrent->gx);
}
for (k2 = k; k2 < aBrent->n; k2++)
{
sl = aBrent->gfx;
s = 0.0;
// minimize along "non-conjugate" directions
lineMin(aBrent, k2, 2, &aBrent->d[k2], &s, aBrent->gfx, FALSE);
if (illc)
{
sz = s + aBrent->z[k2];
s = aBrent->d[k2]*sz*sz;
}
else
s = sl - aBrent->gfx;
if (df < s)
{
df = s;
