dfpopt.c

运筹学非线性规划的经典c语言编写的DFP变尺度法子程序

 #include "hjfgf.c"
 #include "stdlib.h"
 void gradient(double x[],double g[],int n)
 {int i;
  double af,f1,f2,dltx=0.000001;
  for(i=0;i<n;i++)
  g[i]=0;
  f1=objf(x);
  for(i=0;i<n;i++)
  {af=*(x+i);
   *(x+i)=af+dltx;
   f2=objf(x);
   g[i]=(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,*ay[2],*df[2],*zd[2],*zc[2],*zh[2];
  s=(double *)malloc(n*sizeof (double ));
  x=(double *)malloc(n*sizeof (double ));
  for(i=0;i<2;i++)
  {ay[i]=(double*)malloc(n*sizeof(double));
   df[i]=(double*)malloc(n*sizeof(double));
   zd[i]=(double*)malloc(n*sizeof(double));
   zc[i]=(double*)malloc(n*sizeof(double));
   zh[i]=(double*)malloc(n*n*sizeof(double));
  }
  for(i=0;i<n;i++)
  *(ay[0]+i)=xx[i];
  L1:
   k=0;
   for (i=0;i<n;i++)
   for (j=0;j<n;j++)
    {*(zh[0]+i*n+j)=0;
     if(i==j)
      *(zh[0]+i*n+j)=1;
    }
  for(i=0;i<n;i++)
  {*(x+i)=*(ay[0]+i);
   *(ay[1]+i)=*(ay[0]+i);
  }
  gradient(x,df[0],n);
  for(i=0;i<n;i++)
  {*(s+i)=0;
   *(df[1]+i)=*(df[0]+i);
   for(j=0;j<n;j++)
   *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
  }
  L2:
  oneoptim(x,s,h0,epsg,n,ay[0]);
  for(i=0;i<n;i++)
   *(x+i)=*(ay[0]+i);
  gradient(x,df[0],n);
  ae=0;
  for(i=0;i<n;i++)
   ae=ae+(*(df[0]+i))*(*(df[0]+i));
  if(ae<=eps)
   {for(i=0;i<n;i++)
    *(xx+i)=*(x+i);
    free(s);
    free(x);
    for(i=0;i<2;i++)
    {free(ay[i]);
     free(df[i]);
     free(zd[i]);
     free(zc[i]);
     free(zh[i]);
    }
    return(objf(xx));
   }
   if(k==n) goto L1;
   zcc=0;
   for(i=0;i<n;i++)
   {*(zd[0]+i)=*(ay[0]+i)-(*(ay[1]+i));
    *(zd[1]+i)=*(df[0]+i)-(*(df[1]+i));
    *(df[1]+i)=*(df[0]+i);
    zcc=zcc+(*(zd[0]+i))*(*(zd[1]+i));
   }
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh[1]+i*n+j)=*(zd[0]+i)*(*(zd[0]+j))/zcc;
   for(i=0;i<n;i++)
    { *(zc[0]+i)=0;
      for(j=0;j<n;j++)
       *(zc[0]+i)=*(zc[0]+i)+(*(zd[1]+j))*(*(zh[0]+j*n+i));
    }
     zcc=0;
    for(i=0;i<n;i++)
     zcc=zcc+(*(zc[0]+i))*(*(zd[1]+i));
     for(i=0;i<n;i++)
     {*(zc[0]+i)=0;
      *(zc[1]+i)=0;
      for (j=0;j<n;j++)
       {*(zc[0]+i)=*(zc[0]+i)+(*(zh[0]+i*n+j))*(*(zd[1]+j));
        *(zc[1]+i)=*(zc[1]+i)+(*(zh[1]+j))*(*(zh[0]+j*n+i));
       }
     }
    for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    *(zh[0]+i*n+j)=*(zh[0]+i*n+j)+(*(zh[1]+i*n+j))-(*(zc[0]+i))*(*(zc[1]+j))/zcc;
    for(i=0;i<n;i++)
     {*(s+i)=0;
      for(j=0;j<n;j++)
       *(s+i)=*(s+i)-(*(zh[0]+i*n+j))*(*(df[0]+j));
     }
     k=k+1;
     goto L2;
 }