線形計画の単純アルゴリズム

問題の定義:
問題の定義は複雑で、「アルゴリズム導論」の線形計画の一章を見ることを提案した.単純アルゴリズムは、次のような問題を解くために使用されます.
例:
式の最小値を求めます:-2 X 1–3 X 2
引数は次の制約を満たします.
 X1 + X2 = 7
X1  –  2X2 <= 4
X1 >=  0
コンストレイントを標準に変換するには、次の手順に従います.
標準型の条件:
1.目標関数の最大値を求める
2.各引数がゼロ以上(負の制約なし)
3.制約不等式、最小化制約のみ
変換結果は次のとおりです.
max     2X1 – 3X2 + 3X3
そして、以下を満たす.
X1 + X2- X3 <= 7
-X1 – X2+ X3 <= -7
X1 – 2X2+ 2X3 <= 4
X1,  X2, X3 >= 0
 
標準型をリラックス型に変換するには、次の手順に従います.
z = 2X1– 3X2 + 3X3
X4 = 7- X1 - X2 + X3
X5 = -7+ X1 + X2 - X3
X6 = 4- X1 + 2X2 - 2X3
基本変数の下付き集合B={4,5,6}
非基本変数の下付き集合N={1,2,3}
C = [ 2   3  3 ]T
b = [ 7   -7 4 ]T
A        
1
1
-1
-1
-1
1
1   
-2   
2   
 
v = 0
具体的には,A,b,C,vを1つの配列に格納する.
data
    0v   
     2c   
    -3c     
       3c    
7b
1A
1A
-1A
-7b
-1A
-1A
1A
4b
1A
-2A
2A
 
コードは次のとおりです.

/*
*Copyright(c) Computer Science Department of XiaMen University 
*
*Authored by laimingxing on: 2012  03  03      00:14:35 CST
*
* @desc:
*
* @history
*/
#include <iostream>
#include<algorithm>
#include <vector>
#include <string.h>
using namespace std;

const double unbounded = 10000000;
const int MAX = 10;
int n;
double data[MAX][MAX];//A,b,c,v    data  

void PIVOT( vector<int> &N, vector<int> &B, int l ,int e);
void simple( vector<int> &N, vector<int>&B);

int main(int argc, char* argv[])
{
	
//          
	//     
	const int nX = 3;
	//      ,        
	const int nEquation = 4;
	int i = 0, j = 0, k = 0;
	vector<int>B,N;
	double arr[ nEquation + 1 ][nX + 1 ] = {
	{ 2, -3, 3,0},{1, 1, -1, 7}, {-1, -1, 1, -7}, {1, -2, 2, 4}
				};
	
	//N
	for( k = 1; k <= nX; k++)
		N.push_back( k);
	//B
	for( ; k < nX + nEquation; k++)
		B.push_back( k);

	n = nX + nEquation ; 

	memset( data, 0, sizeof(int) * n * n);

	//c
	for( i = 1; i < n; i++)
		if( i <= nX) data[0][i] = arr[0][i - 1];
		else data[0][i] = 0;
	//A
	for( i = nX + 1; i < n; i++)	
	{
		for( j = 1; j <= nX ; j++)
		{
			data[i][j] = arr[ i - nX][j -1];
			//cout << data[i][j] << "\t";
		}
	}
	//b
	for( i = nX + 1; i < n; i++)
		data[i][0] = arr[ i - nX][nX];
 
	simple( N,B);
	return 0;
}
void simple( vector<int> &N, vector<int>&B)
{
	for( int j = 1; j <= n; j++)
	{
		if(data[0][j] <= 0) continue;

		double minBound = unbounded;
		int minIndex = 0;
		for( vector<int>::iterator i = B.begin(); i != B.end(); i++)	
		{
			if( data[*i][j] <= 0)continue;
			double temp = data[*i][0] / data[*i][j];
			if( temp < minBound)
			{
				minBound = temp;
				minIndex = *i;
			}
		}
		if( minBound > unbounded - 1)
			cout <<"Unbounded" << endl;
		else
			PIVOT( N, B,  minIndex, j);
	}	

	for( int i = 1 ; i <= n ; i++)
		if( find(B.begin(), B.end(),i) != B.end() )
			cout << "x"<<i << "  = " << data[i][0] << endl;
		else
			cout << "x"<<i << "  = " << 0 <<endl;
}

//N      ，B     
void PIVOT( vector<int> &N, vector<int> &B, int l ,int e)
{
	data[e][0] = data[l][0] / data[l][e];
	vector<int>::iterator  j;

	for( j = N.begin(); j != N.end(); j++)
	{
		if( *j == e)continue;
		data[e][*j] = data[l][*j] / data[l][e];
	}
	data[e][l] = 1/ data[l][e];

	for( vector<int>::iterator iter = B.begin(); iter != B.end(); iter++)
	{
		if( *iter == l) continue;
		data[*iter][0] -= data[*iter][e] * data[e][0];
		for( vector<int>::iterator it = N.begin(); it != N.end(); it++)
		{
			if( *it == e)continue;
			data[*iter][*it] -= data[*iter][e] * data[e][*it];
		}	
		data[*iter][l] = -1 * data[*iter][l] * data[e][l];
	}
	//Computer the objective function
	data[0][0] += data[0][e] * data[e][0];
	for( j = N.begin(); j != N.end(); j++)
	{
		if( *j == e) continue;
		data[0][*j] -= data[0][e] * data[e][*j];
	}
	data[0][l] = -1 * data[0][e] * data[e][l];
	//Compute new sets of basic and nonbasic variables.
	vector<int>::iterator temp_it = find( N.begin(), N.end(), e);
	if( temp_it != N.end() ) N.erase( temp_it);
	N.push_back(l);

	temp_it = find( B.begin(), B.end(), l);
	if( temp_it != B.end() ) B.erase( temp_it);
	B.push_back(e);

}