ch10.h

将C语言的常用程序集移植到VC开发环境的源代码

/************************************************
 Expect bugs!
 Please use and enjoy, and let me know of any bugs/mods/improvements 
 that you have found/implemented and I will fix/incorporate them into 
 this file. Thank Mr. Xushiliang once again.

					hujinshan@2002.11.3
				Airforce Engineering University
************************************************/

/*****  #include "CH10.h"  极值问题*****/
#ifndef CH10_H_
#define CH10_H_

#include "stdlib.h"
#include "math.h"
#include "stdio.h"

#include "ch2.h"
//*******************************************************************
//作参数的函数指针指向用户自编函数，在具体应用中必须予以确定，函数定义参见文件末
//调用例子: jmax1(x,eps,10,js,mjmax1f);直接将函数名作实参即可
//*******************************************************************
void jmax1(double x[2],double eps,int k,int js[],
		   void (*jmax1f)(double x,double y[2]));//一维极值连分式法

void jmaxn(double x[],int n,double eps,int k,int js[],
		   double (*jmaxnf)(double x[],int n,int j));// n维极值连分式法

int jlplq(double a[],double b[],double c[],int m,int n,double x[]);//为等式约束线性规划问题

int jjsim(int n,double d,double u,double v,double x[],double eps,int k,double xx[],
		  double f[], double (*jjsimf)(double x[],int n));//求n维极值的单形调优法

int jcplx(int n,int m,double a[],double b[],double alpha,double eps,
		  double x[],double xx[],int k,double (*jcplxf)(double x[],int n),
		  void (*jcplxs)(int n,int m,double x[],double c[],double d[],double w[]));//求约束条件下n维极值的复形调优法

//*******************************************************************
void jmax1(double x[2],double eps,int k,int js[],
		   void (*jmax1f)(double x,double y[2]))//一维极值连分式法
{ 
	//extern void jmax1f();
    int i,j,m,jt;
    double xx,h1,h2,dx,y[10],b[10],z[2];
    js[0]=0; jt=1; h2=0.0;
    while (jt==1)
      { j=0;
        while (j<=7)
          { if (j<=2) xx=x[0]+0.01*j;
            else xx=h2;
            jmax1f(xx,z);
            if (fabs(z[1])<eps)
              { jt=0; j=10;}
            else
              { h1=z[1]; h2=xx;
                if (j==0) { y[0]=h1; b[0]=h2;}
                else
                  { y[j]=h1; m=0; i=0;
                    while ((m==0)&&(i<=j-1))
			{ if (fabs(h2-b[i])+1.0==1.0) m=1;
			  else h2=(h1-y[i])/(h2-b[i]);
                          i=i+1;
			}
                    b[j]=h2;
                    if (m!=0) b[j]=1.0e+35;
                    h2=0.0;
                    for (i=j-1; i>=0; i--)
                      h2=-y[i]/(b[i+1]+h2);
                    h2=h2+b[0];
                  }
                j=j+1;
              }
          }
        x[0]=h2;
        js[0]=js[0]+1;
        if (js[0]==k)  jt=0;
      }
    xx=x[0];
    jmax1f(xx,z); x[1]=z[0];
    if (fabs(x[0])<=1.0) dx=1.0e-05;
    else dx=fabs(x[0]*1.0e-05);
    xx=x[0]-dx;
    jmax1f(xx,z); h1=z[0];
    xx=x[0]+dx;
    jmax1f(xx,z); h2=z[0];
    js[1]=-1;
    if ((h1+h2-2.0*x[1])<=0.0) js[1]=1;
    return;
}
/////////////////////////////////////////////////////////////
void jmaxn(double x[],int n,double eps,int k,int js[],
		   double (*jmaxnf)(double x[],int n,int j))// n维极值连分式法
{
	//extern double jmaxnf();
    int i,j,m,l,jt,il;
    double y[10],b[10],p,z,t,h1,h2,f,dx;
    js[0]=0; jt=1; h2=0.0;
    while (jt==1)
      { t=0.0; js[0]=js[0]+1;
        for (i=1; i<=n; i++)
          { f=jmaxnf(x,n,i);
            t=t+fabs(f);
          }
        if (t<eps) jt=0;
        else
          { for (i=0; i<=n-1; i++)
              { il=5;
                while (il!=0)
                  { j=0; t=x[i]; il=il-1;
                    while (j<=7)
                      { if (j<=2) z=t+j*0.01;
                        else z=h2;
                        x[i]=z;
                        f=jmaxnf(x,n,i+1);
                        if (fabs(f)+1.0==1.0)
                          { j=10; il=0;}
                        else
                          { h1=f; h2=z;
                            if (j==0)
                              { y[0]=h1; b[0]=h2;}
                            else
                              { y[j]=h1; m=0; l=0;
                                while ((m==0)&&(l<=j-1))
                                    { p=h2-b[l];
                                      if (fabs(p)+1.0==1.0) m=1;
                                      else h2=(h1-y[l])/p;
                                      l=l+1;
                                    }
                                b[j]=h2;
                                if (m!=0) b[j]=1.0e+35;
                                h2=0.0;
                                for (l=j-1; l>=0; l--)
                                  h2=-y[l]/(b[l+1]+h2);
                                h2=h2+b[0];
                              }
                            j=j+1;
                          }
                      }
                    x[i]=h2; 
                  }
                x[i]=z;
              }
            if (js[0]==k) jt=0;
          }
      }
    js[1]=1;
    f=jmaxnf(x,n,0); x[n]=f;
    dx=0.00001; t=x[0];
    x[0]=t+dx; h1=jmaxnf(x,n,0);
    x[0]=t-dx; h2=jmaxnf(x,n,0);
    x[0]=t;
    t=h1+h2-2.0*f;
    if (t>0.0) js[1]=-1;
    j=1; jt=1;
    while (jt==1)
      { j=j+1; dx=0.00001; jt=0;
        t=x[j-1];
        x[j-1]=t+dx; h2=jmaxnf(x,n,0);
        x[j-1]=t-dx; h1=jmaxnf(x,n,0);
        x[j-1]=t; t=h1+h2-2.0*f;
        if ((t*js[1]<0.0)&&(j<n)) jt=1;
      }
    if (t*js[1]>0.0) js[1]=0;
    return;
}
/////////////////////////////////////////////////////////////
int jlplq(double a[],double b[],double c[],int m,int n,double x[])
{ 
	int i,mn,k,j,l,*js;
    double s,z,dd,y,*p,*d;
    //extern int brinv();
    //extern void brmul();
    js=(int*)malloc(m*sizeof(int));
    p=(double*)malloc(m*m*sizeof(double));
    d=(double*)malloc(m*(m+n)*sizeof(double));
    for (i=0; i<=m-1; i++) js[i]=n+i;
    mn=m+n; s=0.0;
    while (1==1)
      { for (i=0; i<=m-1; i++)
          for (j=0; j<=m-1; j++)
            p[i*m+j]=a[i*mn+js[j]];
        l=brinv(p,m);
        if (l==0)
          { x[n]=s; free(js); free(p);
            free(d); return(l);
          }
        brmul(p,a,m,m,mn,d);
        for (i=0; i<=mn-1; i++) x[i]=0.0;
        for (i=0; i<=m-1; i++)
          { s=0.0;
            for (j=0; j<=m-1; j++)
              s=s+p[i*m+j]*b[j];
            x[js[i]]=s;
          }
        k=-1; dd=1.0e-35;
        for (j=0; j<=mn-1; j++)
          { z=0.0;
            for (i=0; i<=m-1; i++)
              z=z+c[js[i]]*d[i*mn+j];
            z=z-c[j];
            if (z>dd) { dd=z; k=j;}
          }
        if (k==-1)
          { s=0.0;
            for (j=0; j<=n-1; j++)
              s=s+c[j]*x[j];
	    x[n]=s; free(js); free(p);
            free(d); return(1);
          }
        j=-1;
        dd=1.0e+20;
        for (i=0; i<=m-1; i++)
          if (d[i*mn+k]>=1.0e-20)
            { y=x[js[i]]/d[i*mn+k];
              if (y<dd) { dd=y; j=i;}
            }
        if (j==-1) 
          { x[n]=s; free(js); free(p);
            free(d); return(0);
          }
        js[j]=k;
      }
}
/////////////////////////////////////////////////////////////

int jjsim(int n,double d,double u,double v,double x[],double eps,int k,double xx[],
		  double f[], double (*jjsimf)(double x[],int n))//求n维极值的单形调优法
{ 
	int r,g,i,j,l,kk;
    double nn,fe,fr,fl,fg,ft,ff,fs;
    double *xt,*xf,*xe;
    //extern double jjsimf();
    xt=(double*)malloc(n*sizeof(double));
    xf=(double*)malloc(n*sizeof(double));
    xe=(double*)malloc(n*sizeof(double));
    kk=0; nn=1.0*n;
    fr=sqrt(nn+1.0);
    fl=d*(fr-1.0)/(1.414*nn);
    fg=d*(fr+nn-1.0)/(1.414*nn);
    for(i=0; i<=n-1; i++)
    for(j=0; j<=n; j++)
       xx[i*(n+1)+j]=0.0;
    for(i=1; i<=n; i++)
    for(j=0; j<=n-1; j++)
       xx[j*(n+1)+i]=fl;
    for(i=1; i<=n; i++)
       xx[(i-1)*(n+1)+i]=fg;
    for(i=0; i<=n; i++)
      { for(j=0; j<=n-1; j++)
           { xt[j]=xx[j*(n+1)+i]; }
        f[i]=jjsimf(xt,n);
      }
    ft=1.0+eps;
    while ((kk<k)&&(ft>eps))
      { kk=kk+1;
        fr=f[0]; fl=f[0]; r=0; l=0;
        for (i=1; i<=n; i++)
          { if (f[i]>fr) { r=i; fr=f[i];}
            if (f[i]<fl) { l=i; fl=f[i];}
          }
        g=0; fg=f[0];
        j=0;
        if (r==0) { g=1; fg=f[1]; j=1;}
        for (i=j+1; i<=n; i++)