cg.cpp

共梯度算法程序代码.数学上

#include "math.h"
#include"iostream.h"
void grad(int n,double *a,double *b,double eps,double *x)
{ 
	int i,j,k;
	double *p,*r,*s,*q,alpha,beta,d,e;
	p=new double[n];
	r=new double[n];
	s=new double[n];
	q=new double[n];
	for (i=0; i<=n-1; i++)
	{
		x[i]=0.0; p[i]=b[i]; r[i]=b[i]; 
	}
	i=0;
	while (i<=n-1)
	{ 
		for (k=0; k<=n-1; k++)
		{ 
			s[k]=0.0;
			for (j=0; j<=n-1; j++)
				s[k]=s[k]+a[k*n+j]*p[j];
		}
		d=0.0; e=0.0;
		for (k=0; k<=n-1; k++)
		{ 
			d=d+p[k]*b[k]; e=e+p[k]*s[k]; 
		}
		alpha=d/e;
		for (k=0; k<=n-1; k++)
			x[k]=x[k]+alpha*p[k];
		for (k=0; k<=n-1; k++)
		{ 
			q[k]=0.0;
			for (j=0; j<=n-1; j++)
				q[k]=q[k]+a[k*n+j]*x[j];
		}
		d=0.0;
		for (k=0; k<=n-1; k++)
		{
			r[k]=b[k]-q[k]; d=d+r[k]*s[k]; 
		}
		beta=d/e; d=0.0;
		for (k=0; k<=n-1; k++) d=d+r[k]*r[k];
		d=sqrt(d);
		if (d<eps) 
		{ 
			delete []p;
			delete []r;
			delete []s;
			delete []q;
			return;
		}
		for (k=0; k<=n-1; k++)
			p[k]=r[k]-beta*p[k];
		i=i+1;
	}
	delete []p;
	delete []r;
	delete []s;
	delete []q;
	return;
}

void main()
{ 
	int i;
	double eps,x[4];
	static double a[4][4]={
		{5.0,7.0,6.0,5.0},
		{7.0,10.0,8.0,7.0},
		{6.0,8.0,10.0,9.0},
		{5.0,7.0,9.0,10.0}};
		static double b[4]={23.0,32.0,33.0,31.0};
		eps=0.000001;
		grad(4,*a,b,eps,x); 
		for (i=0; i<=3; i++)
			cout<<"x["<<i<<"]="<<x[i]<<endl;		
}