lp_simplex.cpp

数值分析中

# include <iostream.h>
# include <stdlib.h>
# include "rat.h"
# define MAX 100
int n,m,i,j,s,t,L[MAX],R[MAX];
rat z[MAX],c[MAX],a[MAX][MAX],b[MAX],max,second_c[MAX],v,min,temp;
bool B[MAX],flag;
int digit(int a)
{
	if(a==0) return 1;
	int k=0;
	while(a){
		k++;
		a/=10;	
	}
	return k;	
}
int cal(rat x)
{
	int count=0;
	if(x<0) count++;
	if(x.Q!=1) count+=1+digit(x.Q);
	count+=digit(abs(x.P));
	return count;	
}
int output(rat x)
{
	int tmd;
	for(tmd=cal(x);tmd<=6;tmd++) cout<<' ';
	cout<<x;
	return 0;	
}
int table()
{
	for(i=0;i<m;i++){
		for(j=0;j<n;j++) output(a[i][j]);
		cout<<'|';output(b[i]);
		cout<<endl;
	}
	for(i=0;i<n;i++) cout<<"_______";
	cout<<"+_______"<<endl;
	for(i=0;i<n;i++) output(z[i]);
	cout<<'|';output(v);cout<<endl;	
	return 0;
}
int main()
{
	while(cin>>n>>m) {
		//输入a,b,c:系数，目标，约束
		for(i=0;i<m;i++)
		for(j=0;j<n;j++) cin>>a[i][j];
		for(i=0;i<n;i++) cin>>second_c[i];
		for(i=0;i<m;i++) cin>>b[i];
		
		//每行补s,构成第一阶段完整的表
		for(i=0;i<m;i++)
		for(j=0;j<m;j++) 
			if(i==j) a[i][j+n]=1;
			else a[i][j+n]=0;
		//当前基及目标:
		for(i=0;i<n;i++) {
			B[i]=0;
			c[i]=0;
		}
		for(i=n;i<n+m;i++) {
			B[i]=1;
			c[i]=1;
		}
		for(i=0;i<m;i++) {
			R[i]=n+i;
			L[n+i]=i;
		}
		n+=m;
				
		//计算了z 和 v
		//如果是第一阶段，直接加就可以了，可是这样移植性差
		for(i=0;i<n;i++){
			z[i]=0;
			for(j=0;j<m;j++) z[i]=z[i]+c[R[j]]*a[j][i];
			z[i]=z[i]-c[i];
		}
		v=0;
		for(i=0;i<m;i++){
			v=v+c[R[i]]*b[i];
		}
		
		cout<<"第一阶段要解决的问题是:"<<endl;
		cout<<"min z = ";
		for(i=0;i<m-1;i++) cout<<'s'<<i+1<<" + ";
		cout<<'s'<<i+1<<endl;
		table();
		int flag2=0;
		while(1){
			max=z[0];s=0;
			for(i=0;i<n;i++)if(z[i]>max) {
				s=i;
				max=z[i];	
			}
			if(max<0||max==0){
				if(v!=0) {
					flag2=1;
					cout<<"没有可行解"<<endl;	
				}
				break;
			}
			flag=0;min=1000000;
			for(i=0;i<m;i++)if(a[i][s]>0){
				flag=1;
				if((b[i]/a[i][s])<min) {
					min=b[i]/a[i][s];
					t=i;	
				}
			}
			if(!flag) {
				cout<<"unbound"<<endl;
				break;
			}
			//新表
			cout<<"枢元素为：L"<<t+1<<s+1<<" = "<<a[t][s]<<endl;
			cout<<"得新表:"<<endl;
			temp=a[t][s];
			b[t]=b[t]/temp;
			for(i=0;i<n;i++) a[t][i]=a[t][i]/temp;// 系数变1
			for(i=0;i<m;i++)if(i!=t){
				temp=a[i][s];
				b[i]=b[i]-temp*b[t];
				for(j=0;j<n;j++)
					a[i][j]=a[i][j]-temp*a[t][j];
			}	
			temp=z[s];
			for(i=0;i<n;i++) z[i]=z[i]-a[t][i]*temp;
			v=v-b[t]*temp;
			B[s]=1;B[R[t]]=0;L[s]=L[R[t]];R[t]=s;
			table();
			
		}
		if(flag2) break;
		cout<<"下面是第二阶段的运算:"<<endl;
		n-=m;
		//重算z,c和v
		for(i=0;i<n;i++) c[i]=second_c[i];
		for(i=0;i<n;i++){
			z[i]=0;
			for(j=0;j<m;j++) z[i]=z[i]+c[R[j]]*a[j][i];
			z[i]=z[i]-c[i];
		}
		v=0;
		for(i=0;i<m;i++){
			v=v+c[R[i]]*b[i];
		}
		table();
		while(1){
			max=z[0];s=0;
			for(i=0;i<n;i++)if(z[i]>max) {
				s=i;
				max=z[i];	
			}
			if(max<0||max==0){
				cout<<"最优解是："<<v<<endl;
				for(i=0;i<n;i++) if(B[i]) 
					cout<<'x'<<i+1<<" = "<<b[L[i]]<<endl;
				cout<<"其余变量值为0"<<endl;
				break;
			}
			flag=0;min=1000000;
			for(i=0;i<m;i++)if(a[i][s]>0){
				flag=1;
				if((b[i]/a[i][s])<min) {
					min=b[i]/a[i][s];
					t=i;	
				}
			}
			if(!flag) {
				cout<<"没有最优解，即趋于负无穷"<<endl;
				break;
			}
			//新表
			cout<<"枢元素为：L"<<t+1<<s+1<<" = "<<a[t][s]<<endl;
			cout<<"得新表:"<<endl;
			temp=a[t][s];
			b[t]=b[t]/temp;
			for(i=0;i<n;i++) a[t][i]=a[t][i]/temp;// 系数变1
			for(i=0;i<m;i++)if(i!=t){
				temp=a[i][s];
				b[i]=b[i]-temp*b[t];
				for(j=0;j<n;j++)
					a[i][j]=a[i][j]-temp*a[t][j];
			}	
			temp=z[s];
			for(i=0;i<n;i++) z[i]=z[i]-a[t][i]*temp;
			v=v-b[t]*temp;
			B[s]=1;B[R[t]]=0;L[s]=L[R[t]];R[t]=s;
			table();
			
		}	
	}
	return 0;
}