不好意思，忘了。程序如下：
#define N 5
#define v 21
#define fjh1 1e-5
#define fjh2 1e-5
#define fjh3 1e-4
#include <stdio.h>
#include <math.h>
/*函数声明*/
/*不精确的直线搜索*/
double line(double *X0,double *p,double *g0,double *f,double *xx);
/*计算函数的梯度*/
void G(double x[N],double f[v],double g[N]);
/*计算向量的模*/
double Mod(double m[N]);
/*计算行向量与列向量相乘*/
double Multiply1(double m[N],double b[N]);
void Multiply5(double a,double b[N],double MUL5[N],double xx[N]);
/*计算矩阵加法*/
void Add(double m[N][N],double b[N][N],double ADD[N][N]);
/*初始函数以及函数值*/
double Fun(double m[N],double f[v]);
void GETD(FILE *fp);
void main(void)
{
FILE *fp;
if((fp=fopen("共轭梯度算法结果.doc","a+"))==NULL)
{
printf("\n\ncan not open file");
}
printf("******************************************************************\n");
printf("********本算法采用不精确直线搜索法确定函数的最优迭代步长；********\n");
printf("********每次迭代方向向量由共轭梯度算法产生。程序采用C语言编写，***\n");
printf("********在visual c++6.0环境下运行测试通过。
********\n");
printf("******************************************************************\n");
printf("********程序编写：机械电子工程 付金海 0800511 ************************\n");
printf("******************************************************************\n");
fprintf(fp,"******************************************************************\n");
fprintf(fp,"********本算法采用不精确的直线搜索法确定函数的最优迭代步长；********\n");
fprintf(fp,"********每次迭代方向向量由共轭梯度算法产生。程序采用C语言编写，*****\n");
fprintf(fp,"********在visual c++6.0环境下运行测试通过。
********\n");
fprintf(fp,"******************************************************************\n");
fprintf(fp,"********程序编写：机械电子工程 付金海 0800511*************************\n");
fprintf(fp,"******************************************************************\n");
GETD(fp);
}
void GETD(FILE *fp)
{
double x0[N];
double f0=0.0;
double f1=0.0;
double t=0.0;
double val1=0.0;
double val2=0.0;
double val3=0.0;
double x[N];
double xx[N];
double g0[N];
double f[v];
double g[N];
double p[N];
double p1[N];
double p2[N];
double g1[N];
double z2;
double z3;
int k=0;
int i=0;
int j=0;
int step=0;
printf("请输入函数的系数，要保证为正定的21个系数。格式为a*x1*x1+b*x2*x2+c*x3*x3+d*x4*x4+e*x5*x5+g*x1*x2+h*x1*x3+i*x1*x4+j*x1*x5:\n");
for(i=0;i<v;i++)
scanf("%lf",&f[i]);
for(i=0;i<v;i++)
{printf("输入函数的系数为f[%d]=%13lf\n",i,f[i]);
fprintf(fp,"输入函数的系数为f[%d]=%13lf\n",i,f[i]);
}
fprintf(fp,"***************************************************\n");
printf("请输入要搜索的初始点x0：\n");
for(i=0;i<N;i++)
scanf("%lf",&x0[i]);
for(i=0;i<N;i++)
{printf("输入搜索的初始点为x0[%d]=%13lf\n",i,x0[i]);
fprintf(fp,"输入搜索的初始点为x0[%d]=%13lf\n",i,x0[i]);
}
fprintf(fp,"***************************************************\n");
f0=Fun(x0,f);
for(i=0;i<N;i++)
x[i]=x0[i];
G(x,f,g);
for(i=0;i<N;i++)
g0[i]=g[i];
for(i=0;i<N;i++)
{p[i]=(-1)*g0[i];
p1[i]=p[i];
}
k=0;
LOOP: ++step;
val1=Multiply1(g0,g0);
t=line(x0,p1,g0,f,xx);/*这里的t是直线搜素后得到的步长*/
z2=Fun(xx,f); /*这里的z2是直线搜素后得到的函数值*/
G(xx,f,g0); /*这里是直线搜素后得 到的梯度值*/
val2=Multiply1(g0,g0);
printf("===============================================\n");
fprintf(fp,"===============================================\n");
printf("现在函数迭代了第%d步\n",step);
fprintf(fp,"现在函数迭代了第%d步\n",step);
printf("当前函数迭代的最优步长为：%20.16e\n",t);
fprintf(fp,"当前函数迭代的最优步长为：%20.16e\n",t);
printf("当前函数值为：%20.16e\n",z2);
fprintf(fp,"当前函数值为：%20.16e\n",z2);
for(i=0;i<N;i++)
{
fprintf(fp,"xx值的每个分量xx[%d]的值为：%20.16e\n",(i+1),xx[i]);
printf("x值的每个分量xx[%d]的值为：%20.16e\n",(i+1),xx[i]);
}
if(
(Mod(g0)<fjh3) && (fabs(z2-f0)/(fabs(f0)+1)<fjh2) && (fabs(Mod(xx)-Mod(x0)))/(fabs(Mod(x0)+1))<fjh1) /*H准则判别*/
{
fprintf(fp,"\n\n****************************************************************\n");
printf("\n\n****************************************************************\n");
fprintf(fp,"最终的迭代结果为：\n\n");
printf("最终的迭代结果为：\n\n");
fprintf(fp,"函数一共迭代了%d步 \n\n",step);
printf("函数一共迭代了%d步 \n\n",step);
fprintf(fp,"函数的精度要求为：%3.1e\n\n",fjh1);
printf("函数的精度要求为：%3.1e\n\n",fjh1);
fprintf(fp,"\n****************************************************************\n");
for(i=0;i<N;i++)
{
fprintf(fp,"xx值的最终迭代结果中每个分量xx[%d]的值为：%20.16e\n",(i+1),xx[i]);
printf("x向量的最终迭代结果中每个分量xx[%d]的值为：%20.16e\n",(i+1),xx[i]);
}
fprintf(fp,"函数Fun的最优值为:%20.16e\n",z2);
printf("\n函数Fun的最优值为:%20.16e\n",z2);
fprintf(fp,"****************************************************************\n\n");
printf("****************************************************************\n\n");
}
else
{
if(k==N)
{
for(i=0;i<N;i++)
{
x0[i]=xx[i];
f0=z2;
g1[i]=g0[i];
g0[i]=g1[i];
p1[i]=(-1)*g0[i];
}
k=0;
goto LOOP;
}
else
{ z3=(val2)/(val1) ;
for(i=0;i<N;i++)
p1[i]=(-1)*g0[i]+z3*p1[i];
val3=fabs(Multiply1(p1,g0));
if(val3<=fjh1)
{
for(i=0;i<N;i++)
{
x0[i]=xx[i];
g1[i]=g0[i];
g0[i]=g1[i];
p1[i]=(-1)*g0[i];
}
k=0;
f0=z2;
goto LOOP;
}
else if(fabs(val3)<((-1)*fjh1))
/*这里的val3是否要加绝对值呢？书上没加*/
{
k=k+1;
f0=z2;
for(i=0;i<N;i++)
{
x0[i]=xx[i];
g1[i]=g0[i];
g0[i]=g1[i];
p2[i]=p1[i];
p1[i]=p2[i];
}
goto LOOP;
}
else
{
k=k+1;
f0=z2;
for(i=0;i<N;i++)
{
x0[i]=xx[i];
g1[i]=g0[i];
g0[i]=g1[i];
p1[i]=(-1)*p1[i];
}
goto LOOP;
}
}
}
}
/*待计算的函数（正定五元）*/
double Fun(double a[N],double f[v])
{double sum;
sum=(f[0]*a[0]*a[0])+(f[1]*a[1]*a[1])+(f[2]*a[2]*a[2])+(f[3]*a[3]*a[3])+(f[4]*a[4]*a[4])+(f[5]*a[0]*a[1])+(f[6]*a[0]*a[2])+(f[7]*a[0]*a[3])+(f[8]*a[0]*a[4])+(f[9]*a[1]*a[2])+(f[10]*a[1]*a[3])+(f[11]*a[1]*a[4])+(f[12]*a[2]*a[3])+(f[13]*a[2]*a[4])+(f[14]*a[3]*a[4])+(f[15]*a[0])+(f[16]*a[1])+(f[17]*a[2])+(f[18]*a[3])+(f[19]*a[4])+f[20];
return sum;
}
/*计算函数的梯度值（正定五元）*/
void G(double x[N],double f[v],double g[N])
{
double *a=x;
g[0]=2*a[0]*f[0]+f[5]*a[1]+f[6]*a[2]+f[7]*a[3]+f[8]*a[4]+f[15];
g[1]=2*a[1]*f[1]+f[5]*a[0]+f[9]*a[2]+f[10]*a[3]+f[11]*a[4]+f[16];
g[2]=2*a[2]*f[2]+f[6]*a[0]+f[9]*a[1]+f[12]*a[3]+f[13]*a[4]+f[17];
g[3]=2*a[3]*f[3]+f[7]*a[0]+f[10]*a[1]+f[12]*a[2]+f[14]*a[4]+f[18];
g[4]=2*a[4]*f[4]+f[8]*a[0]+f[11]*a[1]+f[13]*a[2]+f[14]*a[3]+f[19];
}
/************不精确的一维搜索******************/
double line(double *x0,double *p,double *g0,double *f,double *xx)
{
double a,f0,z=0.1,t1,t11,t12,f1,f2,f22,f222,l=0.4,b=1e20; /*t1代表初始步长，f1是初始的f0*/
t1=1;
f0=Fun(x0,f);
f1=f0;
f2=Multiply1(p,g0);
/*这里的f2相当于g1,p149*/
loop1: Multiply5(t1,p,x0,xx);
f22=Fun(xx,f);
G(xx,f,g0);
f222=Multiply1(g0,p);
if (f222<(l*f2))
{
a=t1;
t11=(a+b)/2;
t12=2*a;
if (t11>=t12)
{t1=t12;}
else
{t1=t11;}
goto loop1;
}
else
{
if ((f1-f22)<((-1)*z*t1*f2))
{
b=t1;
t1=(a+b)/2;
goto loop1;
}
else
{return t1;}/*这里返回数组的值不很清楚*/
}
}
/*计算行向量与列向量相乘*/
double Multiply1(double m[N],double b[N])
{
double sum=0.0;
int i=0;
for(i=0;i<N;i++)
sum+=m[i]*b[i];
return(sum);
}
/*计算矩阵加法*/
void Add(double m[N][N],double b[N][N],double ADD[N][N])
{
int i=0;
int j=0;
for(i=0;i<N;i++)
for(j=0;j<N;j++)
ADD[i][j]=m[i][j]+b[i][j];
}
/*计算数与矩阵相乘*/
void Multiply5(double a,double b[N],double MUL5[N],double xx[N])
{
int i=0;
int j=0;
double xxx[N];
for(i=0;i<N;i++)
{
xxx[i]=a*b[i];
}
for(j=0;j<N;j++)
{
xx[j]=MUL5[j]+xxx[j];
}
}
/*计算向量的模*/
double Mod(double m[N])
{
int i=0;
double sum=0.0;
for(i=0;i<N;i++)
sum=sum+m[i]*m[i];
return(sqrt(sum));
}
目的是通过共轭梯度法去求一个五元二次函数的最优解。
