C语言实现：共轭梯度法与黄金分割法求解优化问题

#include<stdio.h>
#include <math.h>

#define MAXITER 1000000  //最大的迭代次数
#define alpha 0.8      //变化率
#define E2 0.00000001   //误差

double r=4,kk=0;

/******构造惩罚函数***********/
double p(double x[])
{
double zz;
//if(kk++<10)printf("迭代次数为%f  r=%lf\n \n",kk,r);
zz=x[1]*x[1]+2*x[2]*x[2]-4*x[1]-2*x[1]*x[2];   
return zz;
}
double fun(double f,double g)    //构造惩罚函数
{
double x[3],zz;
x[1]=f;x[2]=g;

zz=x[1]*x[1]+2*x[2]*x[2]-4*x[1]-2*x[1]*x[2];
return zz;
}

/******求最优歩长***********/
double buchang(double xx[],double ff,double ss)
{
double y[5],x[5];

double c=1,a=0,b=0,h=1;//定义c,h(为x初值，h为歩长）
double R=0.618,g=0.001;        //定义缩短率,精度
/*******外推法，求区间**********/
y[1]=fun((xx[1]-c*ff),(xx[2]-c*ss));x[1]=c;
x[2]=(c+h);//x[2]表示的是阿尔发，歩长大小
y[2]=fun((xx[1]-x[2]*ff),(xx[2]-x[2]*ss));
if(y[2]>y[1])
{
h=-h;x[3]=x[1];y[3]=y[1]; 
x[1]=x[2];y[1]=y[2];
x[2]=x[3];y[2]=y[3];   
}
x[3]=(x[2]+h);y[3]=fun((xx[1]-x[3]*ff),(xx[2]-x[3]*ss));
while(y[3]<y[2])
{
h=2*h;
x[1]=x[2]; y[1]=y[2];
x[2]=x[3]; y[2]=y[3];
x[3]=x[2]+h; y[3]=fun((xx[1]-x[3]*ff),(xx[2]-x[3]*ss));
}
if(x[3]<x[1]){a=x[3];b=x[1];}
else {a=x[1];b=x[3];}
//printf("区间为： [%f , %f]\n",a,b);     //搜索区间外推法

/*******黄金分割法**********/
x[1]=(b-R*(b-a)); y[1]=fun((xx[1]-x[1]*ff),(xx[2]-x[1]*ss));
x[2]=(a+R*(b-a)); y[2]=fun((xx[1]-x[2]*ff),(xx[2]-x[2]*ss));
do{
if(y[1]>=y[2])
{