simplx.c

一個線性規劃的程式 可以解出需要解的方程式的最佳解 說明檔在壓縮檔中 manual is in the rar file

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "nrutil.h"
#define EPS 1.0e-6
#define FREEALL free_ivector(l3,1,m);free_ivector(l2,1,m);free_ivector(l1,1,n+1);
double **a;
    
void main()
{
	
	int i,j,m,n,m1,m2,m3,*izrov,*iposv,icase;
	double *x;
	double **dmatrix(int nrl,int nrh,int ncl,int nch);
	double *dvector(int nl,int nh);
	void simplx(int m,int n,int m1,int m2,int m3,int *icase,int izrov[],int iposv[]);
	FILE *fp;
	fp=fopen("input2.txt","r");
	fscanf(fp,"%i %i %i %i %i",&m,&n,&m1,&m2,&m3);
//	m=4;n=4;m1=2;m2=1;m3=1;
	a=dmatrix(1,m+2,1,n+1);
	x=dvector(1,n);
	izrov=ivector(1,n);
	iposv=ivector(1,m);
//	c=dvector(1,n);
	for(i=1;i<=m+1;i++){
		for(j=1;j<=n+1;j++)
			fscanf(fp,"%lf",&a[i][j]);
		fscanf(fp,"\n");
	}
	fclose(fp);
/*	 a[1][1]=0;a[1][2]=1;a[1][3]=1;a[1][4]=3;a[1][5]=-0.5;
	a[2][1]=740;a[2][2]=-1;a[2][3]=0;a[2][4]=-2;a[2][5]=0;
	a[3][1]=0;a[3][2]=0;a[3][3]=-2;a[3][4]=0;a[3][5]=7;
	a[4][1]=0.5;a[4][2]=0;a[4][3]=-1;a[4][4]=1;a[4][5]=-2;
	a[5][1]=9;a[5][2]=-1;a[5][3]=-1;a[5][4]=-1;a[5][5]=-1;
	for(i=1;i<=n;i++)
		c[i]=a[1][i+1]; 
	fp=fopen("A.txt","w");
	fprintf(fp,"%i %i %i %i %i\n",m,n,m1,m2,m3);
	for(i=1;i<=m+1;i++){	
		for(j=1;j<=n+1;j++)
			fprintf(fp,"%g ",a[i][j]);
		fprintf(fp,"\n");
	}
	fclose(fp);
	for(i=2;i<=m+1;i++)	
		for(j=2;j<=n+1;j++)
			a[i][j]=-a[i][j]; */
	simplx(m,n,m1,m2,m3,&icase,izrov,iposv);
	fp=fopen("output.txt","w");
	fprintf(fp,"icase=%i\n",icase);
/*	for(i=1;i<=m+2;i++){	
		for(j=1;j<=n+1;j++)
			fprintf(fp,"%g ",a[i][j]);
		fprintf(fp,"\n");
	}*/
	if(icase==0){
		for(i=1;i<=m;i++)	
			if(iposv[i]<=n)
				x[iposv[i]]=a[i+1][1];
		for(i=1;i<=n;i++)
			if(izrov[i]<=n)
				x[izrov[i]]=0;
		for(i=1;i<=n;i++){
			fprintf(fp,"x[%i]=%g\n",i,x[i]);
		}
		fprintf(fp,"The objective function =  %g\n",a[1][1]);
	}
	fclose(fp);
}

void simplx(int m,int n,int m1,int m2,int m3,int *icase,int izrov[],int iposv[])
{
	int i,ip,ir,is,k,kh,kp,m12,nl1,nl2;
	int *l1,*l2,*l3,*ivector(int nl,int nh);
	double q1,bmax;
	void simp1(int mm,int ll[],int nll,int iabf,int *kp,double *bmax);
	void simp2(int n,int l2[],int nl2,int *ip,int kp,double *q1);
	void simp3(int i1,int k1,int ip,int kp);
	void nrerror(char error_text[]);
	void free_ivector(int *v,int nl,int nh);
    
	if (m != (m1+m2+m3)) nrerror("Bad input constraint counts in SIMPLX");
	l1=ivector(1,n+1);
	l2=ivector(1,m);
	l3=ivector(1,m);
	nl1=n;
	for (k=1;k<=n;k++) l1[k]=izrov[k]=k;
	nl2=m;
	for (i=1;i<=m;i++) {
		if (a[i+1][1] < 0.0) nrerror("Bad input tableau in SIMPLX");
		l2[i]=i;
		iposv[i]=n+i;
	}
	for (i=1;i<=m2;i++) l3[i]=1;
	ir=0;
	if (m2+m3) {
		ir=1;
		for (k=1;k<=(n+1);k++) {
			q1=0.0;
			for (i=m1+1;i<=m;i++) q1 += a[i+1][k];
			a[m+2][k] = -q1;
		}
		do {
			simp1(m+1,l1,nl1,0,&kp,&bmax);
			if (bmax <= EPS && a[m+2][1] < -EPS) {
				*icase = -1;
				FREEALL return;
			} else if (bmax <= EPS && a[m+2][1] <= EPS) {
				m12=m1+m2+1;
				if (m12 <= m) {
					for (ip=m12;ip<=m;ip++) {
						if (iposv[ip] == (ip+n)) {
							simp1(ip,l1,nl1,1,&kp,&bmax);
							if (bmax > 0.0)
								goto one;
						}
					}
				}
				ir=0;
				--m12;
				if (m1+1 <= m12)
					for (i=m1+1;i<=m12;i++)
						if (l3[i-m1] == 1)
							for (k=1;k<=n+1;k++)
								a[i+1][k] = -a[i+1][k];
				break;
			}
			simp2(n,l2,nl2,&ip,kp,&q1);
			if (ip == 0) {
				*icase = -1;
				FREEALL return;
			}
one:		simp3(m+1,n,ip,kp);
			if (iposv[ip] >= (n+m1+m2+1)) {
				for (k=1;k<=nl1;k++)
					if (l1[k] == kp) break;
				--nl1;
				for (is=k;is<=nl1;is++) l1[is]=l1[is+1];
				a[m+2][kp+1] += 1.0;
				for (i=1;i<=m+2;i++) a[i][kp+1] = -a[i][kp+1];
			} else {
				if (iposv[ip] >= (n+m1+1)) {
					kh=iposv[ip]-m1-n;
					if (l3[kh]) {
						l3[kh]=0;
						a[m+2][kp+1] += 1.0;
						for (i=1;i<=m+2;i++)
							a[i][kp+1] = -a[i][kp+1];
					}
				}
			}
			is=izrov[kp];
			izrov[kp]=iposv[ip];
			iposv[ip]=is;
		} while (ir);
	}
	for (;;) {
		simp1(0,l1,nl1,0,&kp,&bmax);
		if (bmax <= 0.0) {
			*icase=0;
			FREEALL return;
		}
		simp2(n,l2,nl2,&ip,kp,&q1);
		if (ip == 0) {
			*icase=1;
			FREEALL return;
		}
		simp3(m,n,ip,kp);
		is=izrov[kp];
		izrov[kp]=iposv[ip];
		iposv[ip]=is;
	}
}
    
#undef EPS
#undef FREEALL

void simp1(int mm,int ll[],int nll,int iabf,int *kp,double *bmax)
{
       int k;
       double test;
       *kp=ll[1];
       *bmax=a[mm+1][*kp+1];
       for (k=2;k<=nll;k++) {
            if (iabf == 0)
                test=a[mm+1][ll[k]+1]-(*bmax);
            else
                test=fabs(a[mm+1][ll[k]+1])-fabs(*bmax);
                if (test > 0.0) {
                    *bmax=a[mm+1][ll[k]+1];
                    *kp=ll[k];
                           }
                }
}

#define EPS 1.0e-6
    
void simp2(int n,int l2[],int nl2,int *ip,int kp,double *q1)
{
	int k,ii,i;
	double qp,q0,q;
    
	*ip=0;
	for (i=1;i<=nl2;i++) {
		if (a[l2[i]+1][kp+1] < -EPS) {
			*q1 = -a[l2[i]+1][1]/a[l2[i]+1][kp+1];
			*ip=l2[i];
			for (i=i+1;i<=nl2;i++) {
				ii=l2[i];
				if (a[ii+1][kp+1] < -EPS) {
					q = -a[ii+1][1]/a[ii+1][kp+1];
					if (q < *q1) {
						*ip=ii;
						*q1=q;
					} else if (q == *q1) {
						for (k=1;k<=n;k++) {
							qp = -a[*ip+1][k+1]/a[*ip+1][kp+1];
							q0 = -a[ii+1][k+1]/a[ii+1][kp+1];
							if (q0 != qp) break;
						}
						if (q0 < qp) *ip=ii;
					}
				}
			}
		}
	}
}

void simp3(int i1,int k1,int ip,int kp)
{
	int kk,ii;
	double piv;
    
	piv=1.0/a[ip+1][kp+1];
	for (ii=1;ii<=i1+1;ii++)
		if (ii-1 != ip) {
			a[ii][kp+1] *= piv;
			for (kk=1;kk<=k1+1;kk++)
				if (kk-1 != kp)
					a[ii][kk] -= a[ip+1][kk]*a[ii][kp+1];
		}
	for (kk=1;kk<=k1+1;kk++)
		if (kk-1 != kp) a[ip+1][kk] *= -piv;
	a[ip+1][kp+1]=piv;
}

void nrerror(char error_text[])
{     
  void exit();    
  printf("Numerical Recipes run-time error...\n");        
  printf("%s\n",error_text);      
  printf("...now exiting to system...\n");        
  exit(1);
}

void free_ivector(int *v,int nl,int nh)
{  
  free((char*) (v+nl));
}

int *ivector(int nl,int nh)
{     
  int *v;
  void nrerror(char error_text[]);
  v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));      
  if (!v) 
    nrerror("allocation failure in ivector()");     
  return v-nl;
}

double **dmatrix(int nrl,int nrh,int ncl,int nch)
{        
  int i;  
  double **m;     
  m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));   
  if (!m) 
     nrerror("allocation failure 1 in dmatrix()");   
  m -= nrl;       
  for(i=nrl;i<=nrh;i++) {         
    m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));          
    if (!m[i]) 
      nrerror("allocation failure 2 in dmatrix()");                
    m[i] -= ncl;    }       
  return m;
}

double *dvector(int nl,int nh)
{  
  double *v;      
  v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));        
  if (!v) 
    nrerror("allocation failure in dvector()");    
  return v-nl;
}