对偶单纯形算法.txt

运筹学单纯形算法

//模型:
//	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 LD_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_m;index++)
		{
			if( ptr_b_Array[index] < 0 ){
				if( find_k == -1 )	find_k = index;
				else if( ptr_b_Array[index] < ptr_b_Array[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_n;index++)
		{
			if( ptr_A_Matrix[ find_k * dim_n + index ] < 0 )
			{
				if( find_l == -1 )
					find_l =index;
				else 
				{
					double comp_index=ptrTestArray[index]/ptr_A_Matrix[find_k*dim_n+index];
					double comp_find_l=ptrTestArray[find_l]/ptr_A_Matrix[find_k*dim_n+find_l];
					if( comp_index < comp_find_l )	find_l=index;
				}
			}
		}
		if(find_l==-1)//存在Pj<=0;无可行解
		{
			delete ptrTestArray;
			delete ptrBaseArray;
			value=1;
			return 0;
		}

		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_bk=ptr_b_Array[find_k];
		double find_k_l=ptr_Matrix_tmp[find_k*dim_n+find_l];

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

			for(int jcount=0; jcount<dim_n; jcount++)
			{
				if( find_k==icount )
					ptr_A_Matrix[icount*dim_n+jcount] /= find_k_l;
				else
				{
					double tmp = ptr_Matrix_tmp[icount*dim_n+find_l]*ptr_Matrix_tmp[find_k*dim_n+jcount];
					tmp /= find_k_l;
					ptr_A_Matrix[icount*dim_n+jcount]-=tmp;
				}
			}
		}
		double find_ck=ptrTestArray[find_l];
		for(index=0;index<dim_n;index++)
		{
			if(index==ptrBaseArray[find_k]) 
			{
				ptrTestArray[index]=find_ck/find_k_l;
				ptrTestArray[index]*=-1;
			}
			else
			{
				double tmp=ptr_Matrix_tmp[find_k*dim_n+index]/find_k_l;
				tmp*=find_ck;
				ptrTestArray[index]-=tmp;
			}
		}
		for(index=0;index<dim_m;index++)
		{
			if(index==find_k)
			{
				ptrBaseArray[index]=find_l;
				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;
}