dlinmin.c

数值算法有原码,英文版

#include "nrutil.h"
#define TOL 2.0e-4 //Tolerance passed to dbrent.
int ncom; //Global variables communicate with df1dim.
float *pcom,*xicom,(*nrfunc)(float []);
void (*nrdfun)(float [], float []);
void dlinmin(float p[], float xi[], int n, float *fret, float (*func)(float []),
			 void (*dfunc)(float [], float []))
/*
Given an n-dimensional point p[1..n] and an n-dimensional direction xi[1..n], moves and
resets p to where the function func(p) takes on a minimum along the direction xi from p,
and replaces xi by the actual vector displacement that p was moved. Also returns as fret
the value of func at the returned location p. This is actually all accomplished by calling the
routines mnbrak and dbrent.
*/
{
	float dbrent(float ax, float bx, float cx,
		float (*f)(float), float (*df)(float), float tol, float *xmin);
	float f1dim(float x);
	float df1dim(float x);
	void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb,
		float *fc, float (*func)(float));
	int j;
	float xx,xmin,fx,fb,fa,bx,ax;
	ncom=n; //Define the global variables.
	pcom=vector(1,n);
	xicom=vector(1,n);
	nrfunc=func;
	nrdfun=dfunc;
	for (j=1;j<=n;j++)
	{
		pcom[j]=p[j];
		xicom[j]=xi[j];
	}
	ax=0.0; //Initial guess for brackets.
	xx=1.0;
	mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
	*fret=dbrent(ax,xx,bx,f1dim,df1dim,TOL,&xmin);
	for (j=1;j<=n;j++) 
	{ //Construct the vector results to return.
		xi[j] *= xmin;
		p[j] += xi[j];
	}
	free_vector(xicom,1,n);
	free_vector(pcom,1,n);
}


float df1dim(float x)
{
	int j;
	float df1=0.0;
	float *xt,*df;
	xt=vector(1,ncom);
	df=vector(1,ncom);
	for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
	(*nrdfun)(xt,df);
	for (j=1;j<=ncom;j++) df1 += df[j]*xicom[j];
	free_vector(df,1,ncom);
	free_vector(xt,1,ncom);
	return df1;
}