dfpmin.cpp

实现共轭梯度法的实例

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

#define ITMAX 1000
#define EPS 1.0e-10

void dfpmin(double p[],int n,double ftol,int *iter,double *fret,\
			double (*func)(double p[]),\
			void (*dfunc)(double p[],double xi[]))
{
	int j,i,its;
	double fp,fae,fad,fac;
	double *xi,*g,*dg,*hdg;
	double **hessin;

	hessin=new double*[n+1];
	for(i=0;i<n+1;i++)
		hessin[i]=new double[n+1];
	xi=new double[n+1];
	g=new double[n+1];
	dg=new double[n+1];
	hdg=new double[n+1];
	fp=(*func)(p);
	(*dfunc)(p,g);
	for (i=1;i<=n;i++) {
		for (j=1;j<=n;j++) hessin[i][j]=0.0;
		hessin[i][i]=1.0;
		xi[i] = -g[i];
	}
	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 hdg,dg,g,xi;
			for(i=0;i<n+1;i++)
				delete hessin[i];
			delete hessin;
			return;
		}
		fp=(*fret);
		for (i=1;i<=n;i++) dg[i]=g[i];
		*fret=(*func)(p);
		(*dfunc)(p,g);
		for (i=1;i<=n;i++) dg[i]=g[i]-dg[i];
		for (i=1;i<=n;i++) {
			hdg[i]=0.0;
			for (j=1;j<=n;j++) hdg[i] += hessin[i][j]*dg[j];
		}
		fac=fae=0.0;
		for (i=1;i<=n;i++) {
			fac += dg[i]*xi[i];
			fae += dg[i]*hdg[i];
		}
		fac=1.0/fac;
		fad=1.0/fae;
		for (i=1;i<=n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
		for (i=1;i<=n;i++)
			for (j=1;j<=n;j++)
				hessin[i][j] += fac*xi[i]*xi[j]
					-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
		for (i=1;i<=n;i++) {
			xi[i]=0.0;
			for (j=1;j<=n;j++) xi[i] -= hessin[i][j]*g[j];
		}
	}
	printf("Too many iterations in DFPMIN");
}

#undef ITMAX
#undef EPS