dfp.cpp

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

#include<iostream.h>
#include<math.h>
#include<time.h>
#define N 50
//cF41                              FLETCBV3  
//*
//*                Initial Point: [1/(n+1),2/(n+1).,...,n/(n+1)].
//*
      void CALF(int n,double x[],double &f)
	  { 
		  double kappa, objscale,p,h, c1,c2;
		  kappa = 1.0;
		  objscale = 100000000.0;
		  h = 1.0/(float(n)+1.0);
		  p = 1.0/objscale;                                 
		  
		  c1 = p*(h*h+2.0)/(h*h);
		  c2 = kappa*p/(h*h);                
                       
		  f=0.50*p*(x[0]*x[0]+x[n-1]*x[n-1]);
      
		  for(int i=0;i<n-1;i++)
			  f = f+0.50*p*pow((x[i]-x[i+1]),2);
      
		  for(i=0;i<n;i++)
			  f = f-c1*x[i]-c2*cos(x[i]);
      }

//*-- Gradient
      void CALG(int n,double x[],double g[])
	  { 
     
		  double kappa, objscale,p,h, c1,c2;
      
		  kappa = 1.0;
		  objscale = 100000000.0;
		  h = 1.0/(float(n)+1.0);
		  p = 1.0/objscale;                                
		  
		  c1 = p*(h*h+2.0)/(h*h);
		  c2 = kappa*p/(h*h);
		  
      g[0] = 2.0*p*x[0] - p*x[1] - c1 + c2*sin(x[0]);
      
     for(int i=1;i<n-1;i++)
        g[i] = -p*(x[i-1]-x[i]) + p*(x[i]-x[i+1]) - c1 + c2*sin(x[i]);

      
      g[n-1] = 2.0*p*x[n-1] - p*x[n-2] - c1 + c2*sin(x[n-1]);            
                                       
	  }
//*
void DOIT1(double H[N][N],double g[N],double p[N]){
	//DOIT1
	int i,j;
	for(i=0;i<N;i++)
		p[i]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			p[i]+=H[i][j]*g[j];

}//DOIT1
double DDOIT(int n,double m1[N],double m2[N]){
	//DOIT
	double sum=0;
	for(int i=0;i<N;i++)
		sum+=m1[i]*m2[i];
	return sum;
}//DOIT
void DOIT2(double y1[N],double y2[N],double temp[N][N]){
	//DOIT2
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp[i][j]=y1[i]*y2[j];
}//DOIT2
void DOIT3(double temp1[N][N],double temp2[N][N],double temp3[N][N]){
	//DOIT3
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i][j]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			for(int k=0;k<N;k++)
			temp3[i][j]+=temp1[i][k]*temp2[k][j];
}//DOIT3
void DOIT4(double temp1[N],double temp2[N][N],double temp3[N]){
	//DOIT4
	int i,j;
	for(i=0;i<N;i++)
		temp3[i]=0;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i]+=temp1[j]*temp2[i][j];
}
int if_g(double g[N]){
	int i=0;
	double e=0.001;
	for(i=0;i<N;i++)
		if(fabs(g[i])>e)return 0;
		return 1;
}//if_g
void INEXLS(int n,double X[N],double D[N],double &ALPHA,int &count1,int &count2){
	//INEXLS
      int  A21;
	  double C1=0.1,C2=0.5;
	  double DD,XX;
	  double F,FO,F1,F2;
	  double TRY;
	  double DG,DGO,DG1,DG2;
	  double T1,T2,T3,T4;
	  double A1,A2;
	  double RATIO1,RATIO2;
	  double G[N],WORK[N];
	  int INFO;
	  int NCALF=0,NCALG=0,NCALFN;
	  int i;
	  
//
//   N      -  INTEGER, NUMBER OF VARIABLES
//   X      -  ARRAY, STORE THE CURRENT ITERATE POINT
//   F      -  FUNCTION VALUE
//   G      -  GRADIENT
//   D      -  ARRARY, THE SEARCH DIRECTION
//   WORK      -  ARRARY, working space
//   C1,C2  -  PARAMETERS FOR LINE SEARCH: 0<C1<C2<1 AND C1<=0.5
//   ALPHA  -  STEPSIZE
//   NCALF  -  NUMBER OF FUNCTION EVALUATION
//   NCALG  -  NUMBER OF GRADIENT EVALUATION
//   IPRINT -  INTEGER, INFORMATION OF CALCULATIONS WILL BE PRINTED IF IPRINT=0
//   INFO   -  INTEGER, ON RETURN:
//                INFO = 0 :   LINE SEARCH SUCCEEDED;
//                INFO = 1 :   AN UPPER-HILL SEARCH DIRECTION IS GIVEN;
//                INFO = 2 :   line search fails after NLSMAX TRIAL STEPSIZES;
//                INFO = 3 :   SEARCH DIRECTION IS TOO SHORT.
//

/*                                    */
/*       step 1  set initial value    */
/*                                    */
      int NLSMAX=50;
      INFO=0;
	  CALF(N,X,F);
	  CALG(N,X,G);
      DGO=DDOIT(N,D,G);

      if(DGO>=0)goto l1100;
      FO=F;
      for(i=0;i<N;i++)
        WORK[i]=X[i];
      A1=0.0;
      F1=F;
      DG1=DGO;
      A2=-1.0;
      A21=0;
      RATIO1=1.0;
      RATIO2=0.1;
      NCALFN=0;
      ALPHA=1.0;
 
//
//       step 2  COMPUTE THE NEW TRIAL POINT
//
l40:      for(i=0;i<N;i++)
        X[i]=WORK[i]+ALPHA*D[i];
      CALF(N,X,F);
      NCALFN=NCALFN+1;
      if(NCALFN>NLSMAX)goto l1300;
      if(F<=FO+C1*ALPHA*DGO)goto l200;
//
//        Step 3   If FIRST INEQUALITY DOES NOT HOLD, 
//           REDUCE STEPSIZE BY QUADRATIC FIT
//
      RATIO1=RATIO1*0.10;
      if(DGO/(1.0+fabs(FO))<-1.0e-16)goto l80;
      DD=DDOIT(N,D,D);
      XX=DDOIT(N,X,X);
      if(ALPHA*sqrt(DD)/(1.00+sqrt(XX))<=1.0e14)goto l1400;
l80:   A2=ALPHA;
      F2=F;
      A21=1;
l90:   TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
      TRY=-DG1/TRY/2.0;
      goto l240; 
//
//       Step 4 IF FIRST INEQUALITY HOLDS, TEST THE SECOND ONE
//
l200:  CALG(N,X,G);
      NCALG=NCALG+1;
      DG=DDOIT(N,D,G);
      if(fabs(DG)<=-C2*DGO)goto l1500;
      if(DG>=0.0)goto l210;
      if(A2<0.0)goto l260;
      RATIO2=RATIO2*0.10;
      A1=ALPHA;
      F1=F;
      DG1=DG;
      RATIO1=0.10;
      if(!A21)goto l90;
      goto l220;
//
//          Step 5 USING CUBIC INTERPOLATION TO FIND 
//                  A NEW POINT IN [A1,A2]
//
l210:  RATIO1=0.10;
      A2=ALPHA;
      F2=F;
      DG2=DG;
      A21=1;
l220:  T1=(F2-F1)/(A2-A1)-DG1;
      T2=DG2-DG1;
      T3=3.0*T1-T2;
      T4=T2-T1-T1;
      TRY=T3*T3-3.00*T4*DG1;
      if(TRY<0.0)goto l230;
      if(sqrt(TRY)<=T3)goto l230;
      TRY=-DG1*(A2-A1)/(sqrt(TRY)+T3);
      goto l240;
l230:  TRY=((F2-F1)/(A2-A1)-DG1)/(A2-A1);
      TRY=-DG1/TRY/2.0;
//
//            TRUNCATE THE STEP IF NECESSARY
//
l240:  if(TRY<RATIO1*(A2-A1))TRY=RATIO1*(A2-A1);
      if(TRY>(1.0-RATIO2)*(A2-A1))TRY=(1.0-RATIO2)*(A2-A1);
      ALPHA=TRY+A1;
      goto l40;
//
//          Step 6   ALPHA IS THE LARGEST POINT TRIED, USE 
//            CUBIC INTERPOLATION TO FIND ANOTHER TRIAL STEP > ALPHA
//
l260:  T1=(F-F1)/(ALPHA-A1)-DG;
      T2=DG1-DG;
      T3=3.0*T1-T2;
      T4=T2-T1-T1;
      TRY=T3*T3-3.0*T4*DG;
      if(TRY<0.0)goto l275;
      if(sqrt(TRY)<=T3)goto l275;
      TRY=-DG*(ALPHA-A1)/(sqrt(TRY)-T3);
      if(TRY<=10.0*(ALPHA-A1))goto l280;
l275:  TRY=10.0*(ALPHA-A1);
l280:  if(A2<0.0)goto l290;
      if(TRY<0.1*(A2-A1))TRY=0.10*(A2-A1);
      if(TRY>0.90*(A2-A1))TRY=0.90*(A2-A1);
l290:  A1=ALPHA;
      F1=F;
      DG1=DG;
      ALPHA=A1+TRY;
      goto l40;


//
//  UPHILL SEARCH DIRECTION IS USED OR PARAMETERS ARE INCORRECT
//
l1100:  //WRITE(6,2050)
      INFO=1;
      goto l1500;
l1300:  //WRITE(6,2070)
      INFO=2;
      goto l1500;
l1400: CALG(N,X,G);
      NCALG=NCALG+1;
     // WRITE(6,2080)
      INFO=3;
l1500:  NCALF=NCALF+NCALFN;
       count1=NCALF+NCALG;
		count2=NCALG;
    
}
void DFP_design(double H[N][N],double s[N],double y[N]){
	//DFP_design
	double t1,t2;
	double temp1[N],temp2[N];
	double temp3[N][N],temp4[N][N],temp5[N][N];
	int i,j;
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			temp3[i][j]=H[i][j];
		DOIT1(H,y,temp1);
		DOIT2(temp1,y,temp4);
		DOIT3(temp4,H,temp5);
		DOIT4(y,temp3,temp2);
		t1=DDOIT(N,temp2,y);
		DOIT2(s,s,temp3);
		t2=DDOIT(N,y,s);
	for(i=0;i<N;i++)
		for(j=0;j<N;j++){
			temp5[i][j]/=t1;
			temp3[i][j]/=t2;
		}
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			H[i][j]=H[i][j]-temp5[i][j]+temp3[i][j];
}
void main()
{
	int i,j;
	int k=1,k1=0;
	int count1=0,count2=0,count3=0;
	clock_t start,finish;
	double f;
	double alpha;
	double p[N];//search 
	double g[N];
	double s[N],y[N],tempg[N],tempx[N];
//step1:
	double x[N];
	double H[N][N];
for(i=0;i<N;i++)
	x[i]=(i+1)/(N+1);
	
	CALG(N,x,g);
	for(i=0;i<N;i++)
		for(j=0;j<N;j++)
			if(i==j)H[i][j]=1;
			else H[i][j]=0;
	start=clock();   //时间起点
step2:
	DOIT1(H,g,p);
	for(i=0;i<N;i++)
			p[i]=-p[i];
//step3:
	for(i=0;i<N;i++)
		tempx[i]=x[i];
	INEXLS(N,tempx,p,alpha,count1,count2);
	k1+=count1;
	count3+=count2;
	//	for(i=0;i<N;i++)
//	cout<<x[i]<<'\t';
	cout<<"alpha:"<<alpha<<'\n';
//step4 
	for(i=0;i<N;i++)
		x[i]=x[i]+alpha*p[i];
//step5:
	if(if_g(g))goto end;
	CALG(N,x,tempg);
	count3++;
	for(i=0;i<N;i++){
		s[i]=alpha*p[i];
		y[i]=tempg[i]-g[i];
	}
	for(i=0;i<N;i++)
		g[i]=tempg[i];
//step6:
	DFP_design(H,s,y);
	k++;
	goto step2;
end:
  finish=clock();  //时间终点
  CALF(N,x,f);
  cout<<"the opt solution:\n";
  for(i=0;i<N;i++)
	cout<<x[i]<<'\t';
  cout<<endl;
  cout<<"The opt value:"<<f<<endl;
  cout<<"number of iter:"<<k<<endl;
  k1=k1+k;
  cout<<"number of elavation:"<<k1<<endl;
  cout<<"NUMBER OF GRADIENT EVALUATION:"<<count3<<endl;
  cout<<"The duration is :"<<double(finish-start)/CLOCKS_PER_SEC<<'m';
}