//模型:
//	max   f＝CX
//	s.t.  AX=b (X>=0, b>=0)
//dim_m：矩陣A的行數
//dim_n：矩陣A的列數
//ptr_C_Array：C向量數組，數組長度為dim_n;
//ptr_b_Array：b向量數組，長度為dim_m;
//ptr_A_Matrix：矩陣A數組，長度為[dim_m*dim_n];
//ptr_X_Array：解向量輸出數組，長度為dim_n;
//value：返回值類型，0為唯一解，1為無界解，2為無窮多解。
//函數返回值為f

/*例子
	單純形法m*n (m= 3, n=5 )
	目標函數C＝[2,1,0,0,0]
	向量b＝[15,24,5]
	矩陣A=
	[ 0,5,1,0,0
	  6,2,0,1,0
	  1,1,0,0,1]
*/

double LP_WorkOut(double * ptr_C_Array,double * ptr_b_Array,
		double  * ptr_A_Matrix,int dim_m,int dim_n,double * & ptr_X_Array,int & value)
{
	int * ptrBaseArray=new int[dim_m];
	for( int tmpcount=0; tmpcount < dim_m; tmpcount++ )	ptrBaseArray[ tmpcount ] = dim_n - dim_m + tmpcount;
	double * ptrTestArray=new double[dim_n];
	for( tmpcount=0; tmpcount < dim_n; tmpcount++ )	ptrTestArray[ tmpcount ] = ptr_C_Array[tmpcount];

	do{
		int  find_k = -1;
		for(int index=0;index<dim_n;index++)
		{
			if( ptrTestArray[index] > 0 ){
				if( find_k == -1 )	find_k = index;
				else if( ptrTestArray[index] > ptrTestArray[find_k] )	find_k=index;
			}
		}
		if( find_k == -1 ) 
		{
			int ZeroNum=0;
			for( int index=0; index<dim_n; index++)	if(ptrTestArray[index]==0) ZeroNum++;
			if( ZeroNum != dim_m ) value=2;
			break;
		}
		
		int find_l=-1;
		for(index=0;index<dim_m;index++)
		{
			//由于已經對 Pj<=0 進行了檢驗,所以必存在某一index
			//使得Matrix[index][find_k]>0
			if( ptr_A_Matrix[ index * dim_n + find_k ] >0 )
			{
				if( find_l == -1 )
					find_l =index;
				else 
				{
					double comp_index=ptr_b_Array[index]/ptr_A_Matrix[index*dim_n+find_k];
					double comp_find_l=ptr_b_Array[find_l]/ptr_A_Matrix[find_l*dim_n+find_k];
					if( comp_index < comp_find_l )	find_l=index;
				}
			}
		}
		if(find_l==-1)//存在Pj<=0;為無界解
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=1;
			return 0;
		}
		//find_l,find_k都已經找到
		double * ptr_Matrix_tmp=new double[dim_m*dim_n];
		int icount=0,jcount=0;
		for(icount=0;icount<dim_m;icount++)
			for(jcount=0;jcount<dim_n;jcount++)
				ptr_Matrix_tmp[icount*dim_n+jcount]=ptr_A_Matrix[icount*dim_n+jcount];
		
		double find_bl=ptr_b_Array[find_l];
		double find_l_k=ptr_Matrix_tmp[find_l*dim_n+find_k];

		for(icount=0; icount<dim_m; icount++)
		{
			if( find_l==icount )
				ptr_b_Array[icount]/=find_l_k;
			else
			{
				double tmp=ptr_Matrix_tmp[icount*dim_n+find_k]*find_bl;
				tmp/=find_l_k;
				ptr_b_Array[icount]-=tmp;
			}

			for(int jcount=0; jcount<dim_n; jcount++)
			{
				if( find_l==icount )
					ptr_A_Matrix[icount*dim_n+jcount] /= find_l_k;
				else
				{
					double tmp = ptr_Matrix_tmp[icount*dim_n+find_k]*ptr_Matrix_tmp[find_l*dim_n+jcount];
					tmp /= find_l_k;
					ptr_A_Matrix[icount*dim_n+jcount]-=tmp;
				}
			}
		}
		double find_ck=ptrTestArray[find_k];
		for(index=0;index<dim_n;index++)
		{
			if(index==ptrBaseArray[find_l]) 
			{
				ptrTestArray[index]=find_ck/find_l_k;
				ptrTestArray[index]*=-1;
			}
			else
			{
				double tmp=ptr_Matrix_tmp[find_l*dim_n+index]/find_l_k;
				tmp*=find_ck;
				ptrTestArray[index]-=tmp;
			}
		}
		for(index=0;index<dim_m;index++)
		{
			if(index==find_l)
			{
				ptrBaseArray[index]=find_k;
				break;
			}
		}
		delete ptr_Matrix_tmp;
	}while(1);
	
	for(int icount=0;icount<dim_m;icount++)
			ptr_X_Array[ ptrBaseArray[icount] ] = ptr_b_Array[icount];

	double WorkOut=0;
	for(icount=0;icount<dim_n;icount++)
		WorkOut+=ptr_C_Array[icount]*ptr_X_Array[icount];

	delete ptrTestArray;
	delete ptrBaseArray;

	return WorkOut;
}