内点惩罚函数法子程.doc

#include hjfgf.c#include stdlib.hvoid gradient(double x,double g,int n)int i;double af,f1,f2,dltx=0.000001;for(i=0;in;i+)gi=0;f1=objf(x);for(i=0;in;i+)af=*(x+i); *(x+i)=af+dltx; f2=objf(x); gi=(f2-f1)/dltx; *(x+i)=af;double oneoptim(double x0,double s,double h0,double epsg,int n,double x)double *a,*b,ff;a=(double *)malloc(n*sizeof (double );b=(double *)malloc(n*sizeof (double );jtf(x0,h0,s,n,a,b);ff=gold(a,b,epsg,n,x);free(a);free(b);return(ff);double dfpopt(double xx,double h0,double eps,double epsg,int n)int i,j,k;double ae,zcc;double *s,*x,*ay2,*df2,*zd2,*zc2,*zh2;s=(double *)malloc(n*sizeof (double );x=(double *)malloc(n*sizeof (double );for(i=0;i2;i+)ayi=(double*)malloc(n*sizeof(double); dfi=(double*)malloc(n*sizeof(double); zdi=(double*)malloc(n*sizeof(double); zci=(double*)malloc(n*sizeof(double); zhi=(double*)malloc(n*n*sizeof(double);for(i=0;in;i+)*(ay0+i)=xxi;L1: k=0; for (i=0;in;i+) for (j=0;jn;j+) *(zh0+i*n+j)=0; if(i=j) *(zh0+i*n+j)=1; for(i=0;in;i+)*(x+i)=*(ay0+i); *(ay1+i)=*(ay0+i);gradient(x,df0,n);for(i=0;in;i+)*(s+i)=0; *(df1+i)=*(df0+i); for(j=0;jn;j+) *(s+i)=*(s+i)-(*(zh0+i*n+j)*(*(df0+j);L2:oneoptim(x,s,h0,epsg,n,ay0);for(i=0;in;i+) *(x+i)=*(ay0+i);gradient(x,df0,n);ae=0;for(i=0;in;i+) ae=ae+(*(df0+i)*(*(df0+i);if(ae=eps) for(i=0;in;i+) *(xx+i)=*(x+i); free(s); free(x); for(i=0;i2;i+) free(ayi); free(dfi); free(zdi); free(zci); free(zhi); return(objf(xx); if(k=n) goto L1; zcc=0; for(i=0;in;i+) *(zd0+i)=*(ay0+i)-(*(ay1+i); *(zd1+i)=*(df0+i)-(*(df1+i); *(df1+i)=*(df0+i); zcc=zcc+(*(zd0+i)*(*(zd1+i); for(i=0;in;i+) for(j=0;jn;j+) *(zh1+i*n+j)=*(zd0+i)*(*(zd0+j)/zcc; for(i=0;in;i+) *(zc0+i)=0; for(j=0;jn;j+) *(zc0+i)=*(zc0+i)+(*(zd1+j)*(*(zh0+j*n+i); zcc=0; for(i=0;in;i+) zcc=zcc+(*(zc0+i)*(*(zd1+i); for(i=0;in;i+) *(zc0+i)=0; *(zc1+i)=0; for (j=0;jn;j+) *(zc0+i)=*(zc0+i)+(*(zh0+i*n+j)*(*(zd1+j); *(zc1+i)=*(zc1+i)+(*(zh1+j)*(*(zh0+j*n+i); for(i=0;in;i+) for(j=0;jn;j+) *(zh0+i*n+j)=*(zh0+i*n+j)+(*(zh1+i*n+j)-(*(zc0+i)*(*(zc1