dixmaane.cpp

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

//cF54                             DIXMAANE                              
//*                             
//*                    Initial point x0=[2.,2.,2...,2.].
//*
#include<iostream.h>
#include<math.h>
#define N 5
void CALF(int n,double x[N],double &f){
        double alpha,beta,gamma,delta;                      
        int  k1,k2,k3,k4;
    
        alpha=1.0;
        beta=0.0;
        gamma=0.1250;
        delta=0.1250;
                      
        k1=1;
        k2=0;
        k3=0;
        k4=1;
           
        int m = n/3;
        
        f=1.0;
        for(int i=0;i<n;i++)
          f=f+alpha*x[i]*x[i]*pow((float(i+1)/float(n)),k1);
   
    
        for(i=0;i<n-1;i++)
          f=f+beta*x[i]*x[i]*(x[i+1]+x[i+1]*x[i+1])*pow((float(i+1)/float(n)),k2);
   
        
        for(i=0;i<2*m;i++)
          f=f+gamma*x[i]*x[i]*(pow(x[i+m],4))*pow((float(i+1)/float(n)),k3);
    
        
       for(i=0;i<m;i++)
          f=f+delta*x[i]*x[i+2*m]*pow((float(i+1)/float(n)),k4);
}
//*-- Gradient    
void CALG(int n,double x[N],double g[N]){
        double alpha,beta,gamma,delta;                       
        int  k1,k2,k3,k4;
    
        alpha=1.0;
        beta=0.0;
        gamma=0.1250;
        delta=0.1250;
                      
        k1=1;
        k2=0;
        k3=0;
        k4=1;
           
        int m=n/3;

        g[0]=2.0*alpha*x[0]*pow((float(1)/float(n)),k1)+2.0*beta*x[0]*pow((x[0]+x[1]*x[1]),2)*pow((float(1)/float(n)),k2)+2.0*gamma*x[0]*pow(x[m],4)*pow((float(1)/float(n)),k3)+delta*x[2*m]*pow((float(1)/float(n)),k4);
       for(int i=1;i<m;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+delta*x[i+2*m]*pow((float(i+1)/float(n)),k4);
        
       for(i=m;i<2*m;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i+1)/float(n)),k2)+2.0*gamma*x[i]*pow(x[i+m],4)*pow((float(i+1)/float(n)),k3)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)*pow((float(i-m+1)/float(n)),k3);
        
       for(i=2*m;i<n-1;i++)
          g[i]=2.0*alpha*x[i]*pow((float(i+1)/float(n)),k1)+2.0*beta*pow(x[i-1],2)*(x[i]+x[i]*x[i])*(1.0+2.0*x[i])*pow((float(i)/float(n)),k2)+2.0*beta*x[i]*pow((x[i+1]+pow(x[i+1],2)),2)*pow((float(i)/float(n)),k2)+4.0*gamma*pow(x[i-m],2)*pow(x[i],3)* pow((float(i-m+1)/float(n)),k3)+delta*x[i-2*m]*pow((float(i-2*m+1)/float(n)),k4);
        
        
        g[n-1]=2.0*alpha*x[n-1]+2.0*beta*x[n-2]*x[n-2]*(x[n-1]+x[n-1]*x[n-1])*(1.0+2.0*x[n-1])*pow((float(n-1)/float(n)),k2)+4.0*gamma*x[2*m-1]*x[2*m-1]*x[n-1]*x[n-1]*x[n-1]*pow((float(2*m)/float(n)),k3)+delta*x[m-1]*pow((float(m)/float(n)),k4);
}
//************************************************************* 


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*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;
	int k=0;
	for(int i=0;i<N;i++)
		if(i%2==0)x[i]=2;
			else x[i]=2;
   CALG(N,x,G);
	for(i=0;i<N;i++)
		work[i]=-G[i];
	wolf1(x,work,G,alpha,c1,c2);
step2:
	if(wolf1(x,work,G,alpha,c1,c2)&&wolf2(x,work,G,alpha,c1,c2)){
	for(i=0;i<N;i++)
		   x[i]=x[i]+alpha*work[i];
	goto step5;}
	if(wolf1(x,work,G,alpha,c1,c2)==0){
		k++;
		goto step3;}
	if(wolf2(x,work,G,alpha,c1,c2)==0){
		k++;
		goto step4;}
step3:
	b=alpha;
	alpha=(a+alpha)/2;
	goto step2;
step4:
	a=alpha; 
	alpha=(2*alpha)>(alpha+b)/2?(alpha+b)/2:(2*alpha);
	goto step2;
step5:
	cout<<"迭带步数"<<k<<'\n';
	cout<<"alpha:"<<alpha<<'\n';

}