乘子法,网上收集到的一个c程序
- #include <stdio.h>
- #include <math.h>
- #define L 10
- void FX(int m,int n,float x[L][L],float z,float s[L][L],float x0[L][L])
- {int i,j;
- for(i=0;i<m;i++) {for(j=0;j<n;j++)
- { x0[i][j]=x[i][j]+z*s[i][j];}}}
- void grad(int m,float input[L],float x[L][L],float g0[L][L])
- {int i,j;
- for(i=0;i<m;i++)
- { for(j=0;j<m;j++) { if(j!=i)
- g0[i][0]=2*input[i]*x[i][0]+input[2]*x[j][0]+input[i+3];}}}
- void qiudao(float input[L],float s[L][L],float x[L][L],float a,float b,float *z)
- {a=input[0]*s[0][0]*s[0][0]+input[1]*s[1][0]*s[1][0]+input[2]*s[0][0]*s[1][0];
- b=2*input[0]*x[0][0]*s[0][0]+2*input[1]*x[1][0]*s[1][0]+input[2]*(x[0][0]*s[1][0]
- +x[1][0]*s[0][0])+input[3]*s[0][0]+input[4]*s[1][0];
- a=2*a;*z=-1*b/a;}
- void gongetidu(float input[L],float x[L][L])
- {int i,j,k=0,stop=0;float s[L][L]={0};float g0[L][L]={0},g1[L][L]={0},a=0,b=0,z=0,w=0;
- while(stop==0)
- { grad(2,input,x,g0);
- if((g0[0][0]*g0[0][0]+g0[1][0]*g0[1][0])<0.000001)
- {stop=1;break;}
- for(i=0;i<2;i++)
- {for(j=0;j<1;j++) {s[i][j]=-1*g0[i][j];}}
- qiudao(input,s,x,a,b,&z);
- FX(2,1,x,z,s,x);
- grad(2,input,x,g1);
- w=(g1[0][0]*g1[0][0]+g1[1][0]*g1[1][0])/(g0[0][0]*g0[0][0]+g0[1][0]*g0[1][0]);
- for(i=0;i<2;i++)
- {for(j=0;j<1;j++){s[i][j]=-1*g1[i][j]+w*s[i][j];}}
- qiudao(input,s,x,a,b,&z);
- FX(2,1,x,z,s,x);
- k++;
- }
- }
- void main()
- {float f[L]={0},x[L][L]={0};
- float r=0.25,a=5,u=0,c=2,e=0.000001,h,h1,h2;
- int i,j,k=0,stop=1;
- while(stop==1)
- { f[0]=0.5+c/2;f[1]=0.166667+c/2;f[2]=c;f[3]=u-c;
- f[4]=u-c;f[5]=(u*u+u*c)/(2*c);h1=x[0][0]+x[1][0]-1;
- gongetidu(f,x);
- if(fabs(x[0][0]+x[1][0]-1)<e)stop=0;
- h2=x[0][0]+x[1][0]-1;
- h=h1/h2;
- if(h>r) c=a*c;
- u=u+c*h2;
- k++;}
- k--;
- printf("k=%d\n",k);
- printf("x*=");
- {for(i=0;i<2;i++)
- {for(j=0;j<1;j++)
- printf("%f ",x[i][j]);}
- printf(" \n");}
-
- }
复制代码 |