/*
File name: conjg.c
Created by: Ljubomir Buturovic
Created: 11/27/2003
Purpose: conjugate gradient function minimizer.
*/
/*
Copyright 2004 Ljubomir J. Buturovic
Permission is hereby granted, free of charge, to any person
obtaining a copy of this software and associated documentation files
(the "Software"), to deal in the Software without restriction,
including without limitation the rights to use, copy, modify, merge,
publish, distribute, sublicense, and/or sell copies of the Software,
and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
/*
The conjugate-gradient optimization function cgfam() in this file is
a C version of the FORTRAN package CG+ (version 1.1), written by
Guanghui Liu, J. Nocedal and R. Waltz. The package can be downloaded
from http://www.ece.northwestern.edu/~rwaltz/software.html.
The CG+ software is in the public domain, however the authors
request that the original reference on which the package is based be
mentioned in all work using CG+. The reference is:
J. C. Gilbert and J. Nocedal. Global Convergence Properties of
Conjugate Gradient Methods for Optimization, (1992) SIAM J. on
Optimization, 2, 1.
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <errno.h>
#include "conjg.h"
/*
Forms the dot product of two vectors.
Uses unrolled loops for increments equal to one.
Jack Dongarra, Linpack, 3/11/78.
*/
static double ddot(int n, double *dx, int incx, double *dy, int incy)
{
int i,ix,iy,m,mp1;
int done;
double dot;
double dtemp;
done = 0;
dot = 0.0;
dtemp = 0.0;
if(n > 0)
{
if ((incx == 1) && (incy == 1))
{
/*
Code for both increments equal to 1.
Clean-up loop.
*/
m = n % 5;
if (m != 0)
{
for (i = 0; i < m; i++)
dtemp = dtemp + dx[i]*dy[i];
if (n < 5)
{
dot = dtemp;
done = 1;
}
}
if (!done)
{
mp1 = m;
for (i = mp1; i < n; i += 5)
dtemp = dtemp + dx[i]*dy[i]+dx[i+1]*dy[i+1]+
dx[i+2]*dy[i+2]+ dx[i+3]*dy[i+3]+dx[i+4]*dy[i+4];
dot = dtemp;
}
}
else
{
/*
Code for unequal increments or equal increments not equal
to 1.
*/
ix = 1;
iy = 1;
if (incx < 0)
ix = (-n+1)*incx+1;
if (incy < 0)
iy = (-n+1)*incy+1;
for (i = 0; i < n; i++)
{
dtemp = dtemp + dx[ix]*dy[iy];
ix = ix+incx;
iy = iy+incy;
}
dot = dtemp;
}
}
return dot;
}
static double min(double x, double y)
{
double z;
z = x;
if (y < x)
z = y;
return z;
}
static double max(double x, double y)
{
double z;
z = x;
if (y > x)
z = y;
return z;
}
static double three_max(double x, double y, double z)
{
double tmax;
tmax = x;
if (y > tmax)
tmax = y;
if (z > tmax)
tmax = z;
return tmax;
}
/*
The purpose of cstepm() is to compute a safeguarded step for a
linesearch and to update an interval of uncertainty for a minimizer
of the function.
The parameter stx contains the step with the least function
value. The parameter stp contains the current step. It is assumed
that the derivative at stx is negative in the direction of the
step. If brackt is set to 1 then a minimizer has been bracketed in
an interval of uncertainty with endpoints stx and sty.
stx, fx, and dx are variables which specify the step, the function,
and the derivative at the best step obtained so far. The derivative
must be negative in the direction of the step, that is, dx and
stp-stx must have opposite signs. On output these parameters are
updated appropriately.
sty, fy, and dy are variables which specify the step, the function,
and the derivative at the other endpoint of the interval of
uncertainty. On output these parameters are updated appropriately.
stp, fp, and dp are variables which specify the step, the function,
and the derivative at the current step. If brackt is set to 1 then
on input stp must be between stx and sty. On output stp is set to
the new step.
brackt is an input/output variable which specifies if a minimizer
has been bracketed. If the minimizer has not been bracketed then on
input brackt must be set to 0. If the minimizer is bracketed then on
output brackt is set to 1.
stpmin and stpmax specify lower and upper bounds for the step.
info is set as follows:
If *info = 1,2,3,4,5, then the step has been computed according to
one of the five cases below. Otherwise *info = 0, and this indicates
improper input parameters.
Argonne National Laboratory. Minpack project. June 1983
Jorge J. More', David J. Thuente
*/
static void cstepm(double *pstx, double *pfx, double *pdx, double *psty,
double *pfy, double *pdy, double *pstp, double fp, double dp,
int *pbrackt, double stpmin, double stpmax, int *info)
{
int bound;
double gamma,p,q,r,s,sgnd,stpc,stpf,stpq,theta;
double stx;
double fx;
double dx;
double sty;
double fy;
double dy;
double stp;
int brackt;
stx = *pstx;
fx = *pfx;
dx = *pdx;;
sty = *psty;
fy = *pfy;
dy = *pdy;
stp = *pstp;
brackt = *pbrackt;
*info = 0;
/*
Check the input parameters for errors.
*/
if ((brackt && (stp <= min(stx,sty) || stp >= max(stx,sty))) ||
(dx*(stp-stx) >= 0.0) || (stpmax < stpmin))
return;
/*
Determine if the derivatives have opposite sign.
*/
sgnd = dp*(dx/fabs(dx));
/*
First case. A higher function value. The minimum is bracketed. If
the cubic step is closer to stx than the quadratic step, the cubic
step is taken, else the average of the cubic and quadratic steps
is taken.
*/
if (fp > fx)
{
*info = 1;
bound = 1;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = three_max(fabs(theta), fabs(dx), fabs(dp));
gamma = s*sqrt((theta/s)*(theta/s) - (dx/s)*(dp/s));
if (stp < stx) gamma = -gamma;
p = (gamma - dx) + theta;
q = ((gamma - dx) + gamma) + dp;
r = p/q;
stpc = stx + r*(stp - stx);
stpq = stx + ((dx/((fx-fp)/(stp-stx)+dx))/2)*(stp - stx);
if (fabs(stpc-stx) < fabs(stpq-stx))
stpf = stpc;
else
stpf = stpc + (stpq - stpc)/2;
brackt = 1;
/*
Second case. A lower function value and derivatives of
opposite sign. The minimum is bracketed. If the cubic step is
closer to stx than the quadratic (secant) step, the cubic step
is taken, else the quadratic step is taken.
*/
}
else if (sgnd < 0.0)
{
*info = 2;
bound = 0;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = three_max(fabs(theta), fabs(dx), fabs(dp));
gamma = s*sqrt((theta/s)*(theta/s) - (dx/s)*(dp/s));
if (stp > stx) gamma = -gamma;
p = (gamma - dp) + theta;
q = ((gamma - dp) + gamma) + dx;
r = p/q;
stpc = stp + r*(stx - stp);
stpq = stp + (dp/(dp-dx))*(stx - stp);
if (fabs(stpc-stp) > fabs(stpq-stp))
stpf = stpc;
else
stpf = stpq;
brackt = 1;
}
/*
Third case. A lower function value, derivatives of the same sign,
and the magnitude of the derivative decreases. The cubic step is
only used if the cubic tends to infinity in the direction of the
step or if the minimum of the cubic is beyond stp. Otherwise the
cubic step is defined to be either stpmin or stpmax. The quadratic
(secant) step is also computed and if the minimum is bracketed
then the step closest to stx is taken, else the step farthest away
is taken.
*/
else if (fabs(dp) < fabs(dx))
{
*info = 3;
bound = 1;
theta = 3*(fx - fp)/(stp - stx) + dx + dp;
s = three_max(fabs(theta), fabs(dx), fabs(dp));
/*
The case gamma = 0 only arises if the cubic does not tend to
infinity in the direction of the step.
*/
gamma = s*sqrt(max(0.0,(theta/s)*(theta/s) - (dx/s)*(dp/s)));
if (stp > stx) gamma = -gamma;
p = (gamma - dp) + theta;
q = (gamma + (dx - dp)) + gamma;
r = p/q;
if ((r < 0.0) && (gamma != 0.0))
stpc = stp + r*(stx - stp);
else if (stp > stx)
stpc = stpmax;
else
stpc = stpmin;
stpq = stp + (dp/(dp-dx))*(stx - stp);
if (brackt)
{
if (fabs(stp-stpc) < fabs(stp-stpq))
stpf = stpc;
else
stpf = stpq;
}
else
{
if (fabs(stp-stpc) > fabs(stp-stpq))
stpf = stpc;
else
stpf = stpq;
}
}
else
{
/*
Fourth case. A lower function value, derivatives of the same
sign, and the magnitude of the derivative does not
decrease. If the minimum is not bracketed, the step is either
stpmin or stpmax, else the cubic step is taken.
*/
*info = 4;
bound = 0;
if (brackt)
{
theta = 3*(fp - fy)/(sty - stp) + dy + dp;
s = three_max(fabs(theta), fabs(dy), fabs(dp));
gamma = s*sqrt((theta/s)*(theta/s) - (dy/s)*(dp/s));
if (stp > sty) gamma = -gamma;
p = (gamma - dp) + theta;
q = ((gamma - dp) + gamma) + dy;
r = p/q;
stpc = stp + r*(sty - stp);
stpf = stpc;
}
else if (stp > stx)
stpf = stpmax;
else
stpf = stpmin;
}
/*
Update the interval of uncertainty. This update does not depend on
the new step or the case analysis above.
*/
if (fp > fx)
{
sty = stp;
fy = fp;
dy = dp;
}
else
{
if (sgnd < 0.0)
{
sty = stx;
fy = fx;
dy = dx;
}
stx = stp;
fx = fp;
dx = dp;
}
/*
Compute the new step and safeguard it.
*/
stpf = min(stpmax,stpf);
stpf = max(stpmin,stpf);
stp = stpf;
if (brackt && bound)
{
if (sty > stx)
stp = min(stx+0.66*(sty-stx),stp);
else
stp = max(stx+0.66*(sty-stx),stp);
}
*pstx = stx;
*pfx = fx;
*pdx = dx;;
*psty = sty;
*pfy = fy;
*pdy = dy;
*pstp = stp;
*pbrackt = brackt;
return;
}
/*
This is a modification of More's line search routine.
The purpose of cvsmod() is to find a step which satisfies a
sufficient decrease condition and a curvature condition. The user
must provide functions which calculate the function and the
gradient.
At each stage the subroutine updates an interval of uncertainty with
endpoints stx and sty. The interval of uncertainty is initially
chosen so that it contains a minimizer of the modified function
f(x+stp*s) - f(x) - ftol*stp*(gradf(x)'s).
If a step is obtained for which the modified function has a
nonpositive function value and nonnegative derivative, then the
interval of uncertainty is chosen so that it contains a minimizer of
f(x+stp*s).
The algorithm is designed to find a step which satisfies the
sufficient decrease condition
f(x+stp*s) <= f(x) + ftol*stp*(gradf(x)'s),
and the curvature condition
abs(gradf(x+stp*s)'s)) <= gtol*abs(gradf(x)'s).
If ftol is less than gtol and if, for example, the function is
bounded below, then there is always a step which satisfies both
conditions. If no step can be found which satisfies both conditions,
then the algorithm usually stops when rounding errors prevent
further progress. In this case stp only satisfies the sufficient
decrease condition.
'n' is the number of variables. 'x' is an array of length n. On
input it must contain the base point for the line search. On output
it contains x + stp*s. 'func' is the function to be optimized,
'gradient' the corresponding gradient. 'g' is an array of length
n. On output it contains the gradient of func at x +
stp*s. 'direction' is an input array of length n which specifies the
search direction. 'stp' is a nonnegative variable. On input stp
contains an initial estimate of a satisfactory step. On output it
contains the final estimate. 'ftol' and 'gtol' are nonnegative input
variables. Termination occurs when the sufficient decrease condition
and the directional derivative condition are satisfied. 'xtol' is a
nonnegative input variable. Termination occurs when the relative
width of the interval of uncertainty is at most xtol. 'stpmin' and
'stpmax' are nonnegative input variables which specify lower and
upper bounds for the step. 'maxfev' is the max. number of function
evaluations per iteration. 'wa' is a work array of length n. 'nfev'
is set to the number of function evaluations.
On output, cvsmod() returns the function value and sets *errc to 0
for success. On error, *errc is set as follows: *errc = EINVAL
improper input parameters, or error code set by
func()/gradient(). *info is set as follows: 1: the sufficient
decrease condition and the directional derivative condition hold; 2:
relative width of the interval of uncertainty is at most xtol; 3:
number of calls to fcn has reached maxfev; 4: the step is at the
lower bound stpmin; 5: the step is at the upper bound stpmax; 6:
rounding errors prevent further progress. There may not be a step
which satisfies the sufficient decrease and curvature
conditions. Tolerances may be too small.
The following two parameters are a modification to the code:
'dg' is the initial directional derivative (in the original code it
was computed in this routine); 'dgout' is the value of the
directional derivative when the Wolfe conditions hold, and an exit
is made to check descent.
Argonne National Laboratory. Minpack project, June 1983. Jorge
J. More', David J. Thuente
*/
static double cvsmod(double *x, double *direction, crit_dfunc *func,
double *g, int n, crit_gradient *gradient, double ftol,
double gtol, double xtol, double stpmin, double stpmax,
int maxfev, int *nfev, double dginit, double *dgout,
int iteration, void *parameters, double *alpha, double *wa,
int *info, int *errc)
{
int infoc;
int j;
int done;
int brackt;
int stage1;
double fval;
double dg,dgm,dgtest,dgx,dgxm,dgy,dgym;
double finit,ftest1,fm,fx,fxm,fy,fym,p5,p66,stx,sty,stp;
double stmin,stmax,width,width1,xtrapf,zero;
double *s;
fval = FLT_MAX;
done = 0;
if (errc)
*errc = 0;
p5 = 0.5;
p66 = 0.66;
xtrapf = 4.0;
zero = 0.0;
if (alpha)
stp = *alpha;
s = direction;
infoc = 1;
if(*info != 1)
{
/*
Check the input parameters for errors.
*/
if ((n <= 0) || (stp <= zero) || (ftol < zero) ||
(gtol < zero) || (xtol < zero) || (stpmin < zero) ||
(stpmax < stpmin) || (maxfev <= 0))
{
if (errc)
*errc = EINVAL;
}
else
{
/*
Compute the initial gradient in the search direction and
check that 'direction' is a descent direction.
*/
if (dginit >= zero)
{
if (alpha)
*alpha = stp;
return fval;
}
/*
Initialize local variables.
*/
brackt = 0;
stage1 = 1;
if (nfev)
*nfev = 0;
finit = func(x, n, iteration, parameters, errc);
if (errc && (*errc))
done = 1;
dgtest = ftol*dginit;
width = stpmax - stpmin;
width1 = width/p5;
for (j = 0; j < n; j++)
wa[j] = x[j];
/*
The variables stx, fx, dgx contain the values of the step,
function, and directional derivative at the best step.
The variables sty, fy, dgy contain the value of the step,
function, and derivative at the other endpoint of the
interval of uncertainty. The variables stp, f, dg contain
the values of the step, function, and derivative at the
current step.
*/
stx = zero;
fx = finit;
dgx = dginit;
sty = zero;
fy = finit;
dgy = dginit;
}
/*
Start of iteration.
*/
while (!done)
{
if (*info != 1)
{
/*
Set the minimum and maximum steps to correspond to the
present interval of uncertainty.
*/
if (brackt)
{
stmin = min(stx,sty);
stmax = max(stx,sty);
}
else
{
stmin = stx;
stmax = stp + xtrapf*(stp - stx);
}
/*
Force the step to be within the bounds stpmax and
stpmin.
*/
stp = max(stp,stpmin);
stp = min(stp,stpmax);
/*
If an unusual termination is to occur then let stp be
the lowest point obtained so far.
*/
if ((brackt && ((stp <= stmin) || (stp >= stmax))) ||
(nfev && (*nfev >= (maxfev-1))) || (infoc == 0) ||
(brackt && ((stmax-stmin) <= xtol*stmax)))
stp = stx;
/*
Evaluate the function and gradient at stp and compute
the directional derivative.
*/
for (j = 0; j < n; j++)
x[j] = wa[j]+stp*s[j];
*info = 0;
fval = func(x, n, iteration, parameters, errc);
if (errc && (*errc))
done = 1;
else
{
gradient(x, n, g, parameters, errc);
if (errc && (*errc))
done = 1;
else
{
if (nfev)
(*nfev)++;
dg = zero;
for (j = 0; j < n; j++)
dg = dg + g[j]*s[j];
ftest1 = finit + stp*dgtest;
/*
Test for convergence.
*/
if ((brackt && ((stp <= stmin) || (stp >= stmax)))
|| (infoc == 0))
*info = 6;
if ((stp == stpmax) && (fval <= ftest1) && (dg <= dgtest))
*info = 5;
if ((stp == stpmin) && ((fval > ftest1) || (dg >= dgtest)))
*info = 4;
if (nfev && (*nfev >= maxfev))
*info = 3;
if (brackt && ((stmax-stmin) <= xtol*stmax))
*info = 2;
/*
More's code has been modified so that at least
one new function value is computed during the
line search (enforcing at least one
interpolation is not easy, since the code may
override an interpolation).
*/
if ((fval <= ftest1) && (fabs(dg) <= gtol*(-dginit)) &&
(nfev && (*nfev > 1)))
*info = 1;
/*
Check for termination.
*/
if (*info != 0)
{
*dgout = dg;
if (alpha)
*alpha = stp;
return fval;
}
/*
In the first stage we seek a step for which
the modified function has a nonpositive value
and nonnegative derivative.
*/
if (stage1 && (fval <= ftest1) && (dg >= min(ftol, gtol)*dginit))
stage1 = 0;
/*
A modified function is used to predict the
step only if we have not obtained a step for
which the modified function has a nonpositive
function value and nonnegative derivative, and
if a lower function value has been obtained
but the decrease is not sufficient.
*/
if (stage1 && (fval <= fx) && (fval > ftest1))
{
/*
Define the modified function and
derivative values.
*/
fm = fval - stp*dgtest;
fxm = fx - stx*dgtest;
fym = fy - sty*dgtest;
dgm = dg - dgtest;
dgxm = dgx - dgtest;
dgym = dgy - dgtest;
/*
Call cstepm() to update the interval of
uncertainty and to compute the new step.
*/
cstepm(&stx, &fxm, &dgxm, &sty, &fym, &dgym, &stp, fm, dgm,
&brackt, stmin, stmax, &infoc);
/*
Reset the function and gradient values for
f.
*/
fx = fxm + stx*dgtest;
fy = fym + sty*dgtest;
dgx = dgxm + dgtest;
dgy = dgym + dgtest;
}
else
/*
Call cstepm() to update the interval of
uncertainty and to compute the new step.
*/
cstepm(&stx, &fx, &dgx, &sty, &fy, &dgy, &stp, fval, dg, &brackt,
stmin, stmax, &infoc);
/*
Force a sufficient decrease in the size of the
interval of uncertainty.
*/
if (brackt)
{
if (fabs(sty-stx) >= p66*width1)
stp = stx + p5*(sty - stx);
width1 = width;
width = fabs(sty-stx);
}
} /* grad() error */
} /* func() error */
/*
End of iteration.
*/
}
}
}
if (alpha)
*alpha = stp;
return fval;
}
/*
Conjugate gradient method for solving unconstrained nonlinear
optimization problems, as described in the paper:
Gilbert, J.C. and Nocedal, J. (1992). "Global convergence properties
of conjugate gradient methods", Siam Journal on Optimization,
vol. 2, pp. 21-42.
cgfam() returns minimum value of function 'func' of 'n' variables,
using the above method. On input, set 'irest' to 0 for no restarts,
1 to restart every n steps; 'eps' is the convergence constant;
'method' is the method code method (1 : Fletcher-Reeves; 2:
Polak-Ribiere, 3: positive Polak-Ribiere ( beta=max{beta,0} );
'itmax' is the maximum allowed number of iterations; 'parameters' is
a pointer to a structure used to pass additional arguments to the
criterion function.
The signatures of func() and grad() must be:
double func(double *x, int n, int iteration, void *parameters, int *errc);
void grad(double *x, int n, double *grad, void *parameters, int *errc);
Here, 'x' is the n-dimensional point at which the function is
evaluated. 'parameters' is a pointer to the structure used to pass
additional arguments to func() and grad(). 'iteration' is the
current iteration, set by cgfam(). In case of successful evaluation,
func() and grad() must set *errc to 0, and non-zero otherwise.
On success, cgfam() returns the minimum value, sets x to the minimum
point, iter to the number of iterations, and iflag and errc to 0; on
error, it sets iflag and/or errc to the appropriate error codes. The
iflag error codes are:
iflag = -1 : line search failure
iflag = -2 : descent was not obtained
The errc error codes are EINVAL for bad arguments, malloc() error
codes and function/gradient evaluation error codes as set by
func()/grad().
*/
double cgfam(int n, double *x, crit_dfunc *func, crit_gradient *grad,
int irest, double eps, int method, void *parameters, int itmax,
int *iter, int *iflag, int *errc)
{
int im,ndes;
int nfun,maxfev,info,i,nfev,nrst,ides;
int done;
int end_loop;
int finish;
int new;
int descent_direction;
double gtol,one,zero,gnorm,stp1,ftol,xtol,stpmin,
stpmax,stp,beta,betafr,betapr,dg0,gg,gg0,dgold,
dgout,dg,dg1;
double f;
double *g;
double *d;
double *gold;
double *wa;
double tlev;
if (n <= 0)
{
if (errc)
*errc = EINVAL;
}
else
{
if (errc)
*errc = 0;
done = 0;
one = 1.0;
zero = 0.0;
finish = 0;
*iflag = 0;
g = malloc(n*sizeof(double));
d = malloc(n*sizeof(double));
gold = malloc(n*sizeof(double));
wa = malloc(n*sizeof(double));
/*
im = number of times betapr was negative for method 2 or number of
times betapr was 0 for method 3
ndes = number of line search iterations after wolfe conditions
were satisfied
*/
*iter = 0;
f = func(x, n, *iter, parameters, errc);
if (errc && (*errc))
done = 1;
else
{
grad(x, n, g, parameters, errc);
if (errc && (*errc))
done = 1;
nfun = 1;
new = 1;
nrst = 0;
im = 0;
ndes = 0;
for (i = 0; i < n; i++)
d[i] = -g[i];
gnorm = sqrt(ddot(n,g,1,g,1));
stp1 = one/gnorm;
/*
parameters for line search routine
----------------------------------
ftol and gtol are nonnegative input variables. termination
occurs when the sufficient decrease condition and the
directional derivative condition are satisfied.
xtol is a nonnegative input variable. termination occurs
when the relative width of the interval of uncertainty
is at most xtol.
stpmin and stpmax are nonnegative input variables which
specify lower and upper bounds for the step.
maxfev is a positive integer input variable. termination
occurs when the number of calls to fcn is at least
maxfev by the end of an iteration.
*/
ftol = 1.0e-4;
gtol = 1.0e-1;
if(gtol <= 1.0e-04)
gtol = 1.0e-02;
xtol = 1.0e-17;
stpmin = 1.0e-20;
stpmax = 1.0e+20;
maxfev = 40;
}
/*
Main iteration loop.
*/
while (!done)
{
(*iter)++;
/*
When nrst > n and irest == 1 then restart.
*/
nrst = nrst+1;
info = 0;
/*
Call the line search routine of More and Thuente (modified
for the CG+ method).
JJ More and D Thuente, Linesearch Algorithms with
Guaranteed Sufficient Decrease. ACM Transactions on
Mathematical Software 20 (1994), pp 286-307.
*/
nfev = 0;
for (i = 0; i < n; i++)
gold[i] = g[i];
dg = ddot(n,d,1,g,1);
dgold = dg;
stp = one;
/*
Shanno-Phua formula for trial step.
*/
if(!new)
stp = dg0/dg;
if (*iter == 1)
stp = stp1;
ides = 0;
new = 0;
descent_direction = 0;
while (!descent_direction)
{
/*
Call to the line search subroutine.
*/
f = cvsmod(x, d, func, g, n, grad, ftol, gtol, xtol, stpmin,
stpmax, maxfev, &nfev, dg, &dgout, *iter, parameters, &stp,
wa, &info, errc);
#if (defined TESTING)
printf("stp: %f\n", stp);
#endif
/*
info is set in cvsmod() as follows:
info = 0 improper input parameters.
info = 1 the sufficient decrease condition and the
directional derivative condition hold.
info = 2 relative width of the interval of uncertainty
is at most xtol.
info = 3 number of calls to fcn has reached maxfev.
info = 4 the step is at the lower bound stpmin.
info = 5 the step is at the upper bound stpmax.
info = 6 rounding errors prevent further progress.
There may not be a step which satisfies the
sufficient decrease and curvature conditions.
Tolerances may be too small.
*/
if (info != 1)
{
*iflag = -1;
done = 1;
descent_direction = 1;
}
else
{
/*
Test if descent direction is obtained for methods 2 and 3.
*/
gg = ddot(n, g, 1, g, 1);
gg0 = ddot(n, g, 1, gold, 1);
betapr = (gg-gg0)/(gnorm*gnorm);
if ((irest == 1) && (nrst > n))
{
nrst = 0;
new = 1;
descent_direction = 1;
}
else if (method == 1)
descent_direction = 1;
else
{
dg1 = -gg+betapr*dgout;
if (dg1 < 0.0)
descent_direction = 1;
else
{
ides++;
if (ides > 5)
{
done = 1;
*iflag = -2;
descent_direction = 1;
}
}
}
}
} /* end while (descent_direction) */
if (!done)
{
/*
Determine correct beta value for method chosen.
im = number of times betapr was negative for method 2
or number of times betapr was 0 for method 3
ndes = number of line search iterations after Wolfe
conditions were satisfied
*/
nfun += nfev;
ndes = ndes + ides;
betafr = gg/(gnorm*gnorm);
if (nrst == 0)
beta = zero;
else
{
if (method == 1)
beta = betafr;
if (method == 2)
beta = betapr;
if (((method == 2) || (method == 3)) && (betapr < 0))
im++;
if (method == 3)
beta = max(zero,betapr);
}
/*
compute the new direction
--------------------------
*/
for (i = 0; i < n; i++)
d[i] = -g[i]+beta*d[i];
dg0 = dgold*stp;
gnorm = sqrt(ddot(n,g,1,g,1));
if (finish)
{
*iflag = 0;
done = 1;
}
else
{
/*
Termination test: if an element of the gradient vector is >
tlev, continue processing.
If all are <= tlev, do one more step and finish.
*/
tlev = eps*(one+fabs(f));
end_loop = 0;
for (i = 0; (i < n) && !end_loop; i++)
if (fabs(g[i]) > tlev)
end_loop = 1;
if (!end_loop)
finish = 1;
}
} /* end if (!done) */
if ((*iter) == itmax)
done = 1;
} /* main iteration loop */
free(g);
free(d);
free(gold);
free(wa);
} /* end if (n <= 0) */
return f;
}
syntax highlighted by Code2HTML, v. 0.9.1