Unconstrained optimization C 源码
<P>#include <stdlib.h><BR>#define Abs(x) ( ((x)<0.)?-(x):(x) )<BR>static double fev(double *x,double *py,double *ps,<BR> double c,double (*func)());<BR>int optmiz(double *x,int n,double (*func)(),double de,double test,int max)<BR>{ double fs,fp,fa,fb,fc,s,sa,sb,sc; int k,m;<BR> double *pd,*ps,*py,*pg,*ph,*p,*q,*r;<BR> pd=(double *)calloc(n*(n+4),sizeof(double));<BR> ps=pd+n; py=ps+n; pg=py+n; ph=pg+n;<BR> for(p=ph,q=ph+n*n; p<q ;p+=n+1) *p=1.;<BR> for(p=x,q=pg,fb=(*func)(x); q<ph ;){<BR> *p+=de; *q++ =((*func)(x)-fb)/de; *p++ -=de;}<BR> for(m=0; m<max ;++m){<BR> for(p=ps,r=ph; p<py ;++p){ *p=0.;<BR> for(q=pg; q<ph ;) *p-= *r++ * *q++; }<BR> fp=fa=fb; sa=sb=0.; sc=1.;<BR> for(;;){ if((fc=fev(x,py,ps,sc,func))>fb) break;<BR> fa=fb; sa=sb; fb=fc; sb=sc; sc*=2.; }<BR> if(sc==1.){ sb=.5;<BR> for(;;){ if((fb=fev(x,py,ps,sb,func))<fa||sb<1.e-3) break;<BR> fc=fb; sc=sb; sb/=2.;} }<BR> for(k=0; k<3 ;++k){ s=(fc-fa)/(sc-sa);<BR> if((fs=(s-(fb-fa)/(sb-sa))/(sc-sb))<0.) break;<BR> if((s=(sa+sc-s/fs)/2.)==sb) s-=(sb-sa)/5.;<BR> fs=fev(x,py,ps,s,func);<BR> if(fs<fb){ if(s<sb){ sc=sb; fc=fb;} else{ sa=sb; fa=fb;}<BR> sb=s; fb=fs; }<BR> else{ if(s<sb){ sa=s; fa=fs;} else{ sc=s; fc=fs;} }<BR> }<BR> for(p=x,r=ps; r<py ;){ *r*=sb; *p++ += *r++;}<BR> if(Abs(fp-fb)<test){ free(pd); return (m+1);}<BR> for(p=x,q=pg,r=pd; q<ph ;){<BR> *p+=de; fa=((*func)(x)-fb)/de; *p++ -=de;<BR> *r++ =fa- *q; *q++ =fa; }<BR> for(p=py,r=ph; p<pg ;++p){ *p=0.;<BR> for(q=pd; q<ps ;) *p+= *r++ * *q++; }<BR> for(p=py,q=ps,r=pd,sa=sb=0.; p<pg ;){<BR> sa+= *r* *p++; sb+= *q++ * *r++; }<BR> sa=1.+sa/sb;<BR> for(p=ps,q=py,r=ph; p<py ;++p,++q){<BR> for(k=0; k<n ;++k)<BR> *r++ +=(*(ps+k)* *p*sa- *(py+k)* *p- *(ps+k)* *q)/sb;<BR> }<BR> }<BR> free(pd); return 0;<BR>}<BR>static double fev(double *x,double *py,double *ps,double c,double (*func)())<BR>{ double *p,*q,*r;<BR> for(p=x,q=py,r=ps; r<py ;) *q++ = *p++ +c* *r++;<BR> return (*func)(py);<BR>}<BR></P><P>说明 ,测试环镜 rh9：</P>
<P>#include "stdio.h"</P>
<P>#include "stdlib.h"</P>
<P>#include "math.h"</P>
<P>#define ND 2</P>
<P><BR>char nm[]="Rosenbrock";<BR>void main(void)<BR>{ double x,de,test;<BR> double func(double *x);<BR> int i,n=ND,max=60;<BR> test=1.e-12; de=1.e-9;<BR> printf(" Test of Optimizer (%s function)\n",nm);<BR> printf(" max iterations = %d\n",max);<BR> printf(" convergence threshold= %8.2e\n",test);<BR> printf(" derivative interval= %8.2e\n",de);<BR> fprintf(stderr," input initial vector (%d components)\n",n);<BR> for(i=0; i<n ;) scanf("%lf",x+i++);<BR> printf(" input vector: ");<BR> for(i=0; i<n ;) printf(" %f",x); printf("\n");<BR> if(optmiz(x,n,func,de,test,max)) printf(" optimal solution\n");<BR> else printf(" convergence failure\n");<BR> printf(" x-vector: ");<BR> for(i=0; i<n ;) printf(" %f",x); printf("\n"); <BR> printf(" minimum F(x)= %e\n",func(x));<BR>}<BR>/* Rosenbrock criteria function */<BR>double func(double *x)<BR>{ double f,y;<BR> y=x; f=1.-y; y=x-y*y;<BR> return 100.*y*y+f*f;<BR>}</P>
<P>////////////////////////////////////////////////////////////////</P> <P>超强</P>
页:
[1]