bdqrtic.cpp

共轭梯度的几种方法对标准函数的测试结果比较 附比较结果

#include<iostream.h>
#include<math.h>
#define N 6
//cF43                              BDQRTIC (CUTE)   
//*                     Initial point x0=[1.,1.,...,1.].
//*
//*
void CALF(int n,double x[N],double &f){
      double t[11000];
      int n4=n-4;
      f=0.0;
      for(int i=0;i<n4;i++)
        t[i]=x[i]*x[i]+2.0*x[i+1]*x[i+1]+3.0*x[i+2]*x[i+2]+4.0*x[i+3]*x[i+3]+5.0*x[n-1]*x[n-1];
      for(i=0;i<n4;i++)
        f=f+(-4.0*x[i]+3.0)*(-4.0*x[i]+3.0)+t[i]*t[i];
}
//*-- Gradient
void CALG(int n,double x[N],double g[N]){
      double t[11000],tsum;
      int n4=n-4;
      for(int i=0;i<n4;i++)
        t[i]=x[i]*x[i]+2.0*x[i+1]*x[i+1]+3.0*x[i+2]*x[i+2]+4.0*x[i+3]*x[i+3]+5.0*x[n-1]*x[n-1];    
      g[0]=-8.0*(-4.0*x[0]+3.0)+(4.0*t[0])*x[0];        
      g[1]=-8.0*(-4.0*x[1]+3.0)+(8.0*t[1]+4.0*t[1])*x[1];
      g[2]=-8.0*(-4.0*x[2]+3.0)+(12.0*t[0]+8.0*t[1]+4.0*t[2])*x[2];
      g[3]=-8.0*(-4.0*x[3]+3.0)+(16.0*t[0]+12.0*t[1]+8.0*t[2]+4.0*t[3])*x[3];
      for(i=4;i<n4;i++)
        g[i]=-8.0*(-4.0*x[i]+3.0)+(16.0*t[i-3]+12.0*t[i-2]+8.0*t[i-1]+4.0*t[i])*x[i];                                                         
      g[n4]=(16.0*t[n4-2]+12.0*t[n4-2]+8.0*t[n4-1])*x[n4]; 
      g[n4+1]=(16.0*t[n4-2]+12.0*t[n4-1])*x[n4+1];
      g[n4+2]=(16.0*t[n4-1])*x[n4+2]; 
      tsum=0.0;
      for(i=0;i<n4;i++)
        tsum=tsum+t[i];
      g[n-1]=20.0*tsum*x[n-1];
}
//************************************************************   

double DOIT(double m1[N],double m2[N]){
	double sum=0;
	for(int i=0;i<N;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
int wolf1(double x[],double work[],double g[],double alpha,double c1,double c2){
	double f1,f2;
	double temp[N];
	CALF(N,x,f1);
	for(int i=0;i<N;i++)
		temp[i]=x[i]+alpha*work[i];
	CALF(N,temp,f2);
	double temp1=DOIT(g,work);
	double temp2=-c1*alpha*temp1;
	if((f1-f2)>=temp2)return 1;
	else return 0;
}//wolf1
int wolf2(double x[],double work[],double g[],double alpha,double c1,double c2){
	double ng[N];
	double temp[N];
	for(int i=0;i<N;i++)
		temp[i]=x[i]+alpha*work[i];
	CALG(N,temp,ng);
	double temp1=DOIT(ng,work);
	double temp2=c2*fabs(DOIT(g,work));
	if(temp1<=temp2)return 1;
	else return 0;
}//wolf2
void main()
{	
	double x[N],G[N],work[N];
	double c1=0.1,c2=0.5;
	double alpha=1;
	double a=0,b=1000;
	double temp,temp1,temp2,temp3,temp_x[N],temp_g[N];
	int k=0;
	for(int i=0;i<N;i++)
		if(i%2==0)x[i]=0;
			else x[i]=0;
    CALG(N,x,G);
	for(i=0;i<N;i++)
		work[i]=-G[i];
	CALF(N,x,temp1);
	temp2=DOIT(G,work);
step2:
    for(i=0;i<N;i++)
		   temp_x[i]=x[i]+alpha*work[i];
	CALF(N,temp_x,temp);
	if(wolf1(x,work,G,alpha,c1,c2))
		goto step3;
	else{
		b=alpha;
		alpha=a+(alpha-a)/(2*(1+(temp1-temp)/((alpha-a)*temp2)));
		k++;
		goto step2;
	}
step3:
	for(i=0;i<N;i++)
		   temp_x[i]=x[i]+alpha*work[i];
	CALG(N,temp_x,temp_g);
	temp3=DOIT(temp_g,work);
	if(wolf2(x,work,G,alpha,c1,c2))
		goto step4;
	else{
		 a=alpha;
		 alpha=alpha+((alpha-a)*temp3)/(temp2-temp3);
		 temp1=temp;
		 temp2=temp3;
		 k++;
		 goto step2;
		
	}
step4:
	cout<<"迭带步数"<<k<<'\n';
	cout<<"alpha:"<<alpha<<'\n';

}