10.6 Conjugate Gradient Methods in Multidimensions 423 instead ofequation(10.6.5)."Wait,"you say."aren't they equal by the orthogonality conditions(10.6.3)?"They are equal for exact quadratic forms.In the real world, however,your function is not exactly a quadratic form.Arriving at the supposed minimum of the quadratic form,you may still need to proceed for another set of iterations.There is some evidence [2]that the Polak-Ribiere formula accomplishes the transition to further iterations more gracefully:When it runs out of steam.it tends to reset h to be down the local gradient,which is equivalent to beginning the conjugate-gradient procedure anew. The following routine implements the Polak-Ribiere variant,which we recom- mend;but changing one program line,as shown,will give you Fletcher-Reeves.The routine presumes the existence of a function func(p),where p[1..n]is a vector of length n,and also presumes the existence of a function dfunc(p,df)that sets 19881992 the vector gradient df [1..n]evaluated at the input point p. The routine calls linmin to do the line minimizations.As already discussed. -00 you may wish to use a modified version of linmin that uses dbrent instead of brent,i.e.,that uses the gradient in doing the line minimizations.See note below. 6 #include <math.h> #include "nrutil.h" (Nort serve #define ITMAX 200 #define EPS 1.0e-10 America computer, THE Here ITMAX is the maximum allowed number of iterations,while EPS is a small number to rectify the special case of converging to exactly zero function value. ART #define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n); Programs void frprmn(float p[],int n,float ftol,int *iter,float *fret, float (*func)(float [])void (*dfunc)(float ]float []) Given a starting point p[1..n],Fletcher-Reeves-Polak-Ribiere minimization is performed on a function func,using its gradient as calculated by a routine dfunc.The convergence tolerance on the function value is input as ftol.Returned quantities are p(the location of the minimum). iter (the number of iterations that were performed),and fret (the minimum value of the to dir function).The routine linmin is called to perform line minimizations. void linmin(float p[,float xi[],int n,float *fret, OF SCIENTIFIC COMPUTING(ISBN f1oat(*func)(f1oat☐))； int j,its; float gg,gam,fp,dggi float *g,*h,*xi; 10-621 g=vector(1,n); v@cambridge.org Numerical Recipes books or 1988-1992 by Numerical Recipes 43108 h=vector(1,n); r1=vector(1n】: fp=(*func)(p); Initializations (*dfunc)(p,xi); for(j=1:j<=n;j++){ g[1]=-x1[]; Software. xi[j]h[j]=g[j]; (outside North America) for (its=1;its<=ITMAX;its++){ Loop over iterations *itersits; linmin(p,xi,n,fret,func); Next statement is the normal return: if (2.0*fabs(*fret-fp)<=ftol*(fabs(*fret)+fabs(fp)+EPS)){ FREEALL return; fp=*fret; (*dfunc)(p,xi); dgg=gg=0.0; for(j=1;j<=n;j++)10.6 Conjugate Gradient Methods in Multidimensions 423 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). instead of equation (10.6.5). “Wait,” you say, “aren’t they equal by the orthogonality conditions (10.6.3)?” They are equal for exact quadratic forms. In the real world, however, your function is not exactly a quadratic form. Arriving at the supposed minimum of the quadratic form, you may still need to proceed for another set of iterations. There is some evidence [2] that the Polak-Ribiere formula accomplishes the transition to further iterations more gracefully: When it runs out of steam, it tends to reset h to be down the local gradient, which is equivalent to beginning the conjugate-gradient procedure anew. The following routine implements the Polak-Ribiere variant, which we recom￾mend; but changing one program line, as shown, will give you Fletcher-Reeves. The routine presumes the existence of a function func(p), where p[1..n] is a vector of length n, and also presumes the existence of a function dfunc(p,df) that sets the vector gradient df[1..n] evaluated at the input point p. The routine calls linmin to do the line minimizations. As already discussed, you may wish to use a modified version of linmin that uses dbrent instead of brent, i.e., that uses the gradient in doing the line minimizations. See note below. #include <math.h> #include "nrutil.h" #define ITMAX 200 #define EPS 1.0e-10 Here ITMAX is the maximum allowed number of iterations, while EPS is a small number to rectify the special case of converging to exactly zero function value. #define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n); void frprmn(float p[], int n, float ftol, int *iter, float *fret, float (*func)(float []), void (*dfunc)(float [], float [])) Given a starting point p[1..n], Fletcher-Reeves-Polak-Ribiere minimization is performed on a function func, using its gradient as calculated by a routine dfunc. The convergence tolerance on the function value is input as ftol. Returned quantities are p (the location of the minimum), iter (the number of iterations that were performed), and fret (the minimum value of the function). The routine linmin is called to perform line minimizations. { void linmin(float p[], float xi[], int n, float *fret, float (*func)(float [])); int j,its; float gg,gam,fp,dgg; float *g,*h,*xi; g=vector(1,n); h=vector(1,n); xi=vector(1,n); fp=(*func)(p); Initializations. (*dfunc)(p,xi); for (j=1;j<=n;j++) { g[j] = -xi[j]; xi[j]=h[j]=g[j]; } for (its=1;its<=ITMAX;its++) { Loop over iterations. *iter=its; linmin(p,xi,n,fret,func); Next statement is the normal return: if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) { FREEALL return; } fp= *fret; (*dfunc)(p,xi); dgg=gg=0.0; for (j=1;j<=n;j++) {