数学建模社区-数学中国

标题: Unconstrained optimization C 源码 [打印本页]

作者: 一生有你 时间: 2005-8-30 10:24
标题: Unconstrained optimization C 源码

#include <stdlib.h>
#define Abs(x) ( ((x)<0.)?-(x)x) )
static double fev(double *x,double *py,double *ps,
double c,double (*func)());
int optmiz(double *x,int n,double (*func)(),double de,double test,int max)
{ double fs,fp,fa,fb,fc,s,sa,sb,sc; int k,m;
double *pd,*ps,*py,*pg,*ph,*p,*q,*r;
pd=(double *)calloc(n*(n+4),sizeof(double));
ps=pd+n; py=ps+n; pg=py+n; ph=pg+n;
for(p=ph,q=ph+n*n; p<q ;p+=n+1) *p=1.;
for(p=x,q=pg,fb=(*func)(x); q<ph ;){
*p+=de; *q++ =((*func)(x)-fb)/de; *p++ -=de;}
for(m=0; m<max ;++m){
for(p=ps,r=ph; p<py ;++p){ *p=0.;
for(q=pg; q<ph ;) *p-= *r++ * *q++; }
fp=fa=fb; sa=sb=0.; sc=1.;
for(;;){ if((fc=fev(x,py,ps,sc,func))>fb) break;
fa=fb; sa=sb; fb=fc; sb=sc; sc*=2.; }
if(sc==1.){ sb=.5;
for(;;){ if((fb=fev(x,py,ps,sb,func))<fa||sb<1.e-3) break;
fc=fb; sc=sb; sb/=2.;} }
for(k=0; k<3 ;++k){ s=(fc-fa)/(sc-sa);
if((fs=(s-(fb-fa)/(sb-sa))/(sc-sb))<0.) break;
if((s=(sa+sc-s/fs)/2.)==sb) s-=(sb-sa)/5.;
fs=fev(x,py,ps,s,func);
if(fs<fb){ if(s<sb){ sc=sb; fc=fb;} else{ sa=sb; fa=fb;}
sb=s; fb=fs; }
else{ if(s<sb){ sa=s; fa=fs;} else{ sc=s; fc=fs;} }
}
for(p=x,r=ps; r<py ;){ *r*=sb; *p++ += *r++;}
if(Abs(fp-fb)<test){ free(pd); return (m+1);}
for(p=x,q=pg,r=pd; q<ph ;){
*p+=de; fa=((*func)(x)-fb)/de; *p++ -=de;
*r++ =fa- *q; *q++ =fa; }
for(p=py,r=ph; p<pg ;++p){ *p=0.;
for(q=pd; q<ps ;) *p+= *r++ * *q++; }
for(p=py,q=ps,r=pd,sa=sb=0.; p<pg ;){
sa+= *r* *p++; sb+= *q++ * *r++; }
sa=1.+sa/sb;
for(p=ps,q=py,r=ph; p<py ;++p,++q){
for(k=0; k<n ;++k)
*r++ +=(*(ps+k)* *p*sa- *(py+k)* *p- *(ps+k)* *q)/sb;
}
}
free(pd); return 0;
}
static double fev(double *x,double *py,double *ps,double c,double (*func)())
{ double *p,*q,*r;
for(p=x,q=py,r=ps; r<py ;) *q++ = *p++ +c* *r++;
return (*func)(py);
}

说明 ,测试环镜 rh9：

#include "stdio.h"

#include "stdlib.h"

#include "math.h"

#define ND 2

char nm[]="Rosenbrock";
void main(void)
{ double x[ND],de,test;
double func(double *x);
int i,n=ND,max=60;
test=1.e-12; de=1.e-9;
printf(" Test of Optimizer (%s function)\n",nm);
printf(" max iterations = %d\n",max);
printf(" convergence threshold= %8.2e\n",test);
printf(" derivative interval= %8.2e\n",de);
fprintf(stderr," input initial vector (%d components)\n",n);
for(i=0; i<n ;) scanf("%lf",x+i++);
printf(" input vector: ");
for(i=0; i<n ;) printf(" %f",x[i++]); printf("\n");
if(optmiz(x,n,func,de,test,max)) printf(" optimal solution\n");
else printf(" convergence failure\n");
printf(" x-vector: ");
for(i=0; i<n ;) printf(" %f",x[i++]); printf("\n");
printf(" minimum F(x)= %e\n",func(x));
}
/* Rosenbrock criteria function */
double func(double *x)
{ double f,y;
y=x[0]; f=1.-y; y=x[1]-y*y;
return 100.*y*y+f*f;
}

////////////////////////////////////////////////////////////////

作者: 一生有你 时间: 2005-8-31 08:52

超强

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)