frprmn.cpp

实现共轭梯度法的实例

#include "stdafx.h"
#include <math.h>
#include <stdio.h>
#include "linmin.h"

#define ITMAX 1000
#define EPS 1.0e-10

void frprmn(double p[],int n,\
			double ftol,int *iter,double *fret,\
			double (*func)(double[]),\
			void (*dfunc)(double[],double[]))
{
	int j,its;
	double gg,gam,fp,dgg;
	double *g,*h,*xi;

	g=new double[n+1];
	h=new double[n+1];
	xi=new double[n+1];
	fp=(*func)(p);
	(*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++) {
		*iter=its;
		linmin(p,xi,n,fret,func);
		if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
			delete xi,g,h;
			return;
		}
		fp=(*func)(p);
		(*dfunc)(p,xi);
		dgg=gg=0.0;
		for (j=1;j<=n;j++) {
			gg += g[j]*g[j];
/*		  dgg += xi[j]*xi[j];	*/
			dgg += (xi[j]+g[j])*xi[j];
		}
		if (gg == 0.0) {
			delete xi,h,g;
			return;
		}
		gam=dgg/gg;
		for (j=1;j<=n;j++) {
			g[j] = -xi[j];
			xi[j]=h[j]=g[j]+gam*h[j];
		}
	}
	printf("Too many iterations in FRPRMN");
}

#undef ITMAX
#undef EPS