/* 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 #include #include #include #include #include #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; }