frprmn.c

数值算法有原码,英文版

#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++) 
		{
			gg += g[j]*g[j];
			/* dgg += xi[j]*xi[j]; */ //This statement for Fletcher-Reeves.
			dgg += (xi[j]+g[j])*xi[j]; //This statement for Polak-Ribiere.
		}
		if (gg == 0.0) 
		{// Unlikely. If gradient is exactly zero then we are already done. 
			FREEALL
			return;
		}
		gam=dgg/gg;
		for (j=1;j<=n;j++) 
		{
			g[j] = -xi[j];
			xi[j]=h[j]=g[j]+gam*h[j];
		}
	}
	nrerror("Too many iterations in frprmn");
}