单纯形解决线形规划问题的代码（c）

单纯形解决线形规划问题的代码（c）

总算完成了，下面是代码。
又：因为运筹学并没有学过，这次凑巧帮管院一个同学写单纯形的代码，
所以就有了这个。是从昨天开始看书的，如果有错误，还请见谅------觉得那个运筹学还真不容易呢，虽比不上咱们的离散：）~~~~~~~~

测试结果：
Successful:
X[3]=77.142860 X[1]=7.142854 X[2]=28.571430
max Z=585.714294
Press any key to continue

/*
*data.txt
*/

5 2 1 0 0
2 3 0 1 0
1 5 0 0 1

170 100 150

10 18 0 0 0

1 0 0
0 1 0
0 0 1

1 0 0
0 1 0
0 0 1

0 0 0

2 3 4

/*
本题的例子是：

max Z=10*x1+18*x2;
s.t. 5*x1+2*x2<=170
2*x1+3*x2<=100
x1+5*x2<=150
x1,x2>=0

本文件的数据是上式经过加入松弛变量后的标准形式的数据

*/

/*
*fmain.cpp
*/
#include <stdio.h>
#include <stdlib.h>
#define M 3
#define N 5

/* 解决形如
max z=cx
s.t. Ax=b,
x>=0
的线形规划问题

 _B[M][M]是B[M][M]的逆阵
*/
float A[M][N],b[M],
c[N],B[M][M],
  cn[M],X[M],Y[M],
  Kao[N],P[M],
  E[M][M],_B[M][M];
/*
解序号
*/
int Xcount[M];
float sita;
/*
进基变量，出基变量
*/
int xin,xout;

int init(FILE *pf);
void computeE();
void compute_B();
void computeXY();
int computeKao();
int computePsita();
void computeBcn();
int succ();
int fall();

void main()
{
 FILE *pf;
float result=0.0;
if(!init(pf))
 {
  printf("Error:can't open data file,please check.../n");
  exit(-1);
 }
 computeXY();
 xin=computeKao();
 while(!succ())
 {
xout=computePsita();
  if(fall())
  {
   printf("Fall...already exit.../n");
   exit(0);
  }
  computeBcn();
  computeE();
  compute_B();
  computeXY();
xin=computeKao();
 }
 printf("Successful:/n");
 for(int i=0;i<M;i++)
  printf("X[%d]=%f/t",Xcount[i]+1,X[i]);

 for(int i=0;i<M;i++)
result+=c[Xcount[i]]*X[i];
 printf("/nmax Z=%f/n",result);
}
/*
用data.txt文件初始化，包括：
A,b,c,B,_B,cn,Xcount
*/
int init(FILE *pf)
{
 pf=fopen("data.txt","r");
 if(pf==NULL){
  fclose(pf);
  return 0;
 }
 /*input A[M][N]*/
for(int i=0;i<M;i++)
  for(int j=0;j<N;j++){
   fscanf(pf,"%f",&A[i][j]);
  }

/*input b[M]*/
 for(int i=0;i<M;i++){
  fscanf(pf,"%f",&b[i]);
 }

/*input c[N]*/
 for(int i=0;i<N;i++){
  fscanf(pf,"%f",&c[i]);
 }

/*input B[M][M]*/
 for(int i=0;i<M;i++)
  for(int j=0;j<M;j++)
   fscanf(pf,"%f",&B[i][j]);

/*input _B[M][M]*/
 for(int i=0;i<M;i++)
  for(int j=0;j<M;j++){
   fscanf(pf,"%f",&_B[i][j]);
  }

/*input cn[M]*/
 for(int i=0;i<M;i++)
  fscanf(pf,"%f",&cn[i]);

 /* input Xcount[M] */
 for(int i=0;i<M;i++){
  fscanf(pf,"%d",&Xcount[i]);
 }

 /*其它*/
 for(int i=0;i<M;i++)
  P[i]=X[i]=Y[i]=0;

 for(int i=0;i<N;i++)
  Kao[i]=0;

  fclose(pf);
 return 1;
}
/*
计算初等变换矩阵E
*/
void computeE()
{
for(int i=0;i<M;i++)
  for(int j=0;j<M;j++)
  {
   E[i][j]=0;
   if(i==j)
    E[i][j]=1;
  }

 for(int i=0;i<M;i++)
  E[i][xout]=-P[i]/P[xout];
 E[xout][xout]=1/P[xout];
}
/*
计算B逆矩阵_B
*/
void compute_B()
{
 float __B[M][M];
 for(int i=0;i<M;i++)
  for(int j=0;j<M;j++)
   __B[i][j]=0;

 for(int i=0;i<M;i++)
for(int j=0;j<M;j++)
  for(int k=0;k<M;k++)
__B[i][j]+=E[i][k]*_B[k][j];

 for(int i=0;i<M;i++)
  for(int j=0;j<M;j++){
   _B[i][j]=__B[i][j];
  }

}

/*
计算 X=_B*b ,Y=_B*cn
*/
void computeXY()
{
 for(int i=0;i<M;i++)
  P[i]=X[i]=Y[i]=0;

 for(int i=0;i<M;i++){
  for(int j=0;j<M;j++)
   X[i]+=_B[i][j]*b[j];
  }
 for(int i=0;i<M;i++){
  for(int j=0;j<M;j++)
   Y[i]+=cn[j]*_B[j][i];
 }

}
/* 计算检验数 kao=c-cn*_B*A
xin<-return
*/
int computeKao()
{
 float midval=0.0;
 int k;
 for(int i=0;i<N;i++){
  for(int j=0;j<M;j++)
   midval+=Y[j]*A[j][i];
  Kao[i]=c[i]-midval;
  midval=0;
 }
 midval=Kao[0];
 k=0;
for(int j=0;j<N;j++)
  if(Kao[j]>midval){
   midval=Kao[j];
   k=j;
 }

 return k;
}
/*计算进基列向量P,判断出基变量
xout<-return
*/
int computePsita()
{
 int k=0;
 for(int i=0;i<M;i++){
  for(int j=0;j<M;j++)
   P[i]+=_B[i][j]*A[j][xin];
 }

 sita=X[0]/P[0];
 for(int j=0;j<M;j++)
  if(sita>X[j]/P[j])
  {
   sita=X[j]/P[j];
   k=j;
  }

 return k;
}
/*
计算矩阵B,cn
*/
void computeBcn()
{
 for(int i=0;i<M;i++)
  B[i][xout]=A[i][xin];
 cn[xout]=c[xin];

 Xcount[xout]=xin;
}
/*
判断是否成功
*/
int succ()
{
int k=1;
 for(int i=0;i<N;i++)
  k=k&&(Kao[i]<=0.0);

 return k;
}
/*
算法失败，退出
*/
int fall()
{
int k=1;
for(int i=0;i<N;i++)
  k=k&&(P[i]<=0.0);
return k;
}