DFP变尺度法优化子程序

风花雪月

本程序包含5个C文件
mdfp.c
dfpopt.c
funct.c
jtf.c
hjfgf.c

本程序由heizi友情提供

风花雪月

mdfp.c代码如下：

1. #include "dfpopt.c"
   
2. main()
   
3. {double x[]={1,2};
   
4. double ff;
   
5.   ff=dfpopt(x,0.2,0.0001,0.000001,2);
   
6.   printf("\nx[0]=%f,x[1]=%f,ff=%f",x[0],x[1],ff);
   
7.   getch();
   
8. }

复制代码

dfpopt.c代码如下：

1. #include "hjfgf.c"
   
2. #include "stdlib.h"
   
3. void gradient(double x[],double g[],int n)
   
4. {int i;
   
5.   double af,f1,f2,dltx=0.000001;
   
6.   for(i=0;i<n;i++)
   
7.   g[i]=0;
   
8.   f1=objf(x);
   
9.   for(i=0;i<n;i++)
   
10.   {af=*(x+i);
    
11.    *(x+i)=af+dltx;
    
12.    f2=objf(x);
    
13.    g[i]=(f2-f1)/dltx;
    
14.    *(x+i)=af;
    
15.   }
    
16. }
    
17. double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
    
18. {double *a,*b,ff;
    
19.   a=(double *)malloc(n*sizeof (double ));
    
20.   b=(double *)malloc(n*sizeof (double ));
    
21.   jtf(x0,h0,s,n,a,b);
    
22.   ff=gold(a,b,epsg,n,x);
    
23.   free(a);
    
24.   free(b);
    
25.   return(ff);
    
26. }
    
27. double dfpopt(double xx[],double h0,double eps,double epsg,int n)
    
28. {int i,j,k;
    
29.   double ae,zcc;
    
30.   double *s,*x,*ay[2],*df[2],*zd[2],*zc[2],*zh[2];
    
31.   s=(double *)malloc(n*sizeof (double ));
    
32.   x=(double *)malloc(n*sizeof (double ));
    
33.   for(i=0;i<2;i++)
    
34.   {ay[i]=(double*)malloc(n*sizeof(double));
    
35.    df[i]=(double*)malloc(n*sizeof(double));
    
36.    zd[i]=(double*)malloc(n*sizeof(double));
    
37.    zc[i]=(double*)malloc(n*sizeof(double));
    
38.    zh[i]=(double*)malloc(n*n*sizeof(double));
    
39.   }
    
40.   for(i=0;i<n;i++)
    
41.   *(ay[0]+i)=xx[i];
    
42.   L1:
    
43.    k=0;
    
44.    for (i=0;i<n;i++)
    
45.    for (j=0;j<n;j++)
    
46.     {*(zh[0]+i*n+j)=0;
    
47.      if(i==j)
    
48.       *(zh[0]+i*n+j)=1;
    
49.     }
    
50.   for(i=0;i<n;i++)
    
51.   {*(x+i)=*(ay[0]+i);
    
52.    *(ay[1]+i)=*(ay[0]+i);
    
53.   }
    
54.   gradient(x,df[0],n);
    
55.   for(i=0;i<n;i++)
    
56.   {*(s+i)=0;
    
57.    *(df[1]+i)=*(df[0]+i);
    
58.    for(j=0;j<n;j++)
    
59.    *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
    
60.   }
    
61.   L2:
    
62.   oneoptim(x,s,h0,epsg,n,ay[0]);
    
63.   for(i=0;i<n;i++)
    
64.    *(x+i)=*(ay[0]+i);
    
65.   gradient(x,df[0],n);
    
66.   ae=0;
    
67.   for(i=0;i<n;i++)
    
68.    ae=ae+(*(df[0]+i))*(*(df[0]+i));
    
69.   if(ae<=eps)
    
70.    {for(i=0;i<n;i++)
    
71.     *(xx+i)=*(x+i);
    
72.     free(s);
    
73.     free(x);
    
74.     for(i=0;i<2;i++)
    
75.     {free(ay[i]);
    
76.      free(df[i]);
    
77.      free(zd[i]);
    
78.      free(zc[i]);
    
79.      free(zh[i]);
    
80.     }
    
81.     return(objf(xx));
    
82.    }
    
83.    if(k==n) goto L1;
    
84.    zcc=0;
    
85.    for(i=0;i<n;i++)
    
86.    {*(zd[0]+i)=*(ay[0]+i)-(*(ay[1]+i));
    
87.     *(zd[1]+i)=*(df[0]+i)-(*(df[1]+i));
    
88.     *(df[1]+i)=*(df[0]+i);
    
89.     zcc=zcc+(*(zd[0]+i))*(*(zd[1]+i));
    
90.    }
    
91.    for(i=0;i<n;i++)
    
92.     for(j=0;j<n;j++)
    
93.     *(zh[1]+i*n+j)=*(zd[0]+i)*(*(zd[0]+j))/zcc;
    
94.    for(i=0;i<n;i++)
    
95.     { *(zc[0]+i)=0;
    
96.       for(j=0;j<n;j++)
    
97.        *(zc[0]+i)=*(zc[0]+i)+(*(zd[1]+j))*(*(zh[0]+j*n+i));
    
98.     }
    
99.      zcc=0;
    
100.     for(i=0;i<n;i++)
     
101.      zcc=zcc+(*(zc[0]+i))*(*(zd[1]+i));
     
102.      for(i=0;i<n;i++)
     
103.      {*(zc[0]+i)=0;
     
104.       *(zc[1]+i)=0;
     
105.       for (j=0;j<n;j++)
     
106.        {*(zc[0]+i)=*(zc[0]+i)+(*(zh[0]+i*n+j))*(*(zd[1]+j));
     
107.         *(zc[1]+i)=*(zc[1]+i)+(*(zh[1]+j))*(*(zh[0]+j*n+i));
     
108.        }
     
109.      }
     
110.     for(i=0;i<n;i++)
     
111.     for(j=0;j<n;j++)
     
112.     *(zh[0]+i*n+j)=*(zh[0]+i*n+j)+(*(zh[1]+i*n+j))-(*(zc[0]+i))*(*(zc[1]+j))/zcc;
     
113.     for(i=0;i<n;i++)
     
114.      {*(s+i)=0;
     
115.       for(j=0;j<n;j++)
     
116.        *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
     
117.      }
     
118.      k=k+1;
     
119.      goto L2;
     
120. }

复制代码

风花雪月

funct.c代码如下：

1. #include "stdio.h"
   
2. #include "stdlib.h"
   
3. #include "math.h"
   
4. double objf(double x[])
   
5. {double ff;
   
6. ff=x[0]*x[0]+x[1]*x[1]-x[0]*x[1]-10*x[0]-4*x[1]+60;
   
7. return(ff);
   
8. }

复制代码

jtf.c代码如下：

1. #include "funct.c"
   
2. void jtf(double x0[],double h0,double s[],int n,double a[],double b[])
   
3. {int i;
   
4. double *x[3],h,f1,f2,f3;
   
5. for(i=0;i<3;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7.   h=h0;
   
8. for(i=0;i<n;i++)
   
9.   *(x[0]+i)=x0[i];
   
10. f1=objf(x[0]);
    
11. for(i=0;i<n;i++)
    
12.   *(x[1]+i)=*(x[0]+i)+h*s[i];
    
13.   f2=objf(x[1]);
    
14. if(f2>=f1)
    
15.   {h=-h0;
    
16.     for(i=0;i<n;i++)
    
17.     *(x[2]+i)=*(x[0]+i);
    
18.    f3=f1;
    
19.     for(i=0;i<n;i++)
    
20.     {*(x[0]+i)=*(x[1]+i);
    
21.      *(x[1]+i)=*(x[2]+i);
    
22.     }
    
23.    f1=f2;
    
24.    f2=f3;
    
25.    }
    
26.    for(;;)
    
27.    {h=2*h;
    
28.      for(i=0;i<n;i++)
    
29.      *(x[2]+i)=*(x[1]+i)+h*s[i];
    
30.    f3=objf(x[2]);
    
31.    if(f2<f3) break;
    
32.    else
    
33.     { for(i=0;i<n;i++)
    
34.        {*(x[0]+i)=*(x[1]+i);
    
35.         *(x[1]+i)=*(x[2]+i);
    
36.        }
    
37.       f1=f2;
    
38.       f2=f3;
    
39.     }
    
40.    }
    
41.    if(h<0)
    
42.     for(i=0;i<n;i++)
    
43.     {a[i]=*(x[2]+i);
    
44.      b[i]=*(x[0]+i);
    
45.     }
    
46.    else
    
47.     for(i=0;i<n;i++)
    
48.     {a[i]=*(x[0]+i);
    
49.      b[i]=*(x[2]+i);
    
50.      }
    
51.    for(i=0;i<3;i++)
    
52.    free(x[i]);
    
53. }

复制代码

hjfgf.c代码如下：

1. #include "jtf.c"
   
2. double gold(double a[],double b[],double eps,int n,double xx[])
   
3. {int i;
   
4. double f1,f2,*x[2],ff,q,w;
   
5. for(i=0;i<2;i++)
   
6.   x[i]=(double *)malloc(n*sizeof(double));
   
7. for(i=0;i<n;i++)
   
8.   {*(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
   
9.    *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
   
10.   }
    
11.   f1=objf(x[0]);
    
12.   f2=objf(x[1]);
    
13.   do
    
14.    {if(f1>f2)
    
15.      {for(i=0;i<n;i++)
    
16.       {b[i]=*(x[0]+i);
    
17.        *(x[0]+i)=*(x[1]+i);
    
18.        }
    
19.      f1=f2;
    
20.      for(i=0;i<n;i++)
    
21.       *(x[1]+i)=a[i]+0.382*(b[i]-a[i]);
    
22.      f2=objf(x[1]);
    
23.      }
    
24.     else
    
25.      { for(i=0;i<n;i++)
    
26.        {a[i]=*(x[1]+i);
    
27.        *(x[1]+i)=*(x[0]+i);}
    
28.      f2=f1;
    
29.     for(i=0;i<n;i++)
    
30.      *(x[0]+i)=a[i]+0.618*(b[i]-a[i]);
    
31.     f1=objf(x[0]);
    
32.      }
    
33.   q=0;
    
34.   for(i=0;i<n;i++)
    
35.    q=q+(b[i]-a[i])*(b[i]-a[i]);
    
36.   w=sqrt(q);
    
37.   }while(w>eps);
    
38.   for(i=0;i<n;i++)
    
39.    xx[i]=0.5*(a[i]+b[i]);
    
40.   ff=objf(xx);
    
41.   for(i=0;i<2;i++)
    
42.   free(x[i]);
    
43.   return(ff);
    
44. }

复制代码

nicoleyong

xiexie

lzjld502

xiexiele  hehe:@D

jiaqh

thank you

bosiqi

谢谢分享