求线性规划“单纯形”法源程序

[canjian]

实在太急，明天就要

[LeeMaRS]

以前收集的一个源代码, 作者已经记不清楚了...

#include <iostream.h>
#include <conio.h>
#include <math.h>

typedef int BOOL;
#define TRUE 1
#define FALSE 0
typedef double REAL;
#define ZERO 1e-10

//矩阵求逆
BOOL Inv( REAL ** a, int n );
BOOL Inv(REAL * a,int n);

//矩阵相乘
void Damul(REAL * a,REAL * b,size_t m,size_t n,size_t k,REAL * c);

//线形规划--修正的单纯形法,返回真表示计算成功，返回假表示计算失败
BOOL Line_Optimize(REAL * A, REAL * B, REAL * C, int m, int n,
REAL * Result, REAL * X, int * Is);
//A-----------等式约束的系数矩阵
//B-----------等式约束右端常数
//C-----------所求最小值表达式的各项系数
//m-----------系数矩阵行数
//n-----------系数矩阵列
//Result------返回最小值
//X-----------返回计算结果
//Is----------作为基底的X索引,计算结束后返回各X值所在索引

template <class T>
inline void ExChange( T& a, T& b )
{
T temp=a;
a = b;
b = temp;
}


BOOL Inv( REAL ** a, int n )
{
REAL d;
int i,j,k;
int success = FALSE;
int * is=new int [n];
int * js=new int [n];

for (k=0;k<n;k++){
d=0.0;
for(i=k;i<n;i++){
for(j=k;j<n;j++){
if (fabs(a[i][j])>d){
d=fabs(a[i][j]);
is[k]=i;
js[k]=j;
}
}
}

if( d<ZERO ) goto Clear;

for( j=0; j<n; j++ )ExChange( a[k][j], a[is[k]][j] );

for( i=0; i<n; i++ )ExChange( a[i][k], a[i][js[k]] );

a[k][k]=1/a[k][k];

for( j=0; j<n; j++ ){
if( j!=k )a[k][j]*=a[k][k];
}

for( i=0; i<n; i++){
if ( i!=k ){
for( j=0; j<n; j++ )
if( j!=k )a[i][j]-=a[i][k]*a[k][j];
}
}

for( i=0; i<n; i++ ){
if ( i!=k ){
a[i][k]*=((-1.0)*a[k][k]);
}
}

} //end for

for( k=n-1; k>=0; k-- ){
for( j=0; j<n; j++ )ExChange( a[k][j], a[js[k]][j] );

for( i=0; i<n; i++ )ExChange( a[i][k], a[i][is[k]] );
}
success = TRUE;
Clear:
delete [] is;
delete [] js;
return success;
}


BOOL Inv(REAL * a,int n)
{
REAL **kma = new REAL*[n];
for(int i=0;i<n;i++){
kma[i]=a+i*n;
}
BOOL ret=Inv( kma,n);
delete [] kma;
return ret;
}

void Damul(REAL * a,REAL * b,size_t m,size_t n,size_t k,REAL * c)
{
int i,j,l,u;
for (i=0; i<=m-1; i++){
for (j=0; j<=k-1; j++){
u=i*k+j;
c[u]=0.0;
for (l=0; l<=n-1; l++){
c[u] += a[i*n+l]*b[l*k+j];
}
}
}
return;
}



BOOL Line_Optimize(REAL * A, REAL * B, REAL * C, int m, int n,
REAL * Result, REAL * X, int * Is)
{
REAL r;
int i,j,k;
int Success = FALSE;
REAL* b = new REAL [m*m];
REAL* MatTmp = new REAL [m*m];
REAL* Mat1 = new REAL [m];
REAL* Mat2 = new REAL [m];
REAL* E = new REAL [m*m];

for (i=0; i<m; i++){
for (j=0; j<m; j++){
b[i*m+j] = A[i*n+Is[j]];
}
}

if (!Inv(b,m)){
goto Release;
}

Damul(b,B,m,m,1,X);


for(;;){
for (i=0; i<m; i++){
Mat2[i] = C[Is[i]];
}
Damul(Mat2,b,1,m,m,Mat1);

for (i=0; i<n; i++){
for (j=0; j<m; j++){
Mat2[j] = A[j*n+i];
}
Damul(Mat1,Mat2,1,m,1,&r);
r = C[i] - r;
if ( r<-ZERO ) {
break;
}
}

if (i >= n){
*Result = 0;
for (i=0; i<m; i++){
*Result += C[Is[i]]*X[i];
}
Success = TRUE;
goto Release;
}

Damul(b, Mat2, m, m, 1, Mat1);

r = 1E10;
j = -1;
for (k=0; k<m; k++){
if (Mat1[k]>ZERO){
REAL temp = X[k]/Mat1[k];
if (temp<r){
r = temp;
j = k;
}
}
}

if (j < 0) {
Success = FALSE;
goto Release;
}

for (k=0; k<m*m; k++){
E[k] = 0;
}

for (k=0; k<m; k++){
E[k*m+k] = 1;
}

for (k=0; k<m; k++){
E[k*m+j] = -Mat1[k]/Mat1[j];
}
E[j*m+j] = 1/Mat1[j];

Is[j] = i;

Damul(E,b,m,m,m,MatTmp);
Damul(E,X,m,m,1,Mat2);

for (i=0; i<m*m; i++){
b[i] = MatTmp[i];
}

for (i=0; i<m; i++){
X[i] = Mat2[i];
}
}

Release:
delete [] E;
delete [] Mat2;
delete [] Mat1;
delete [] MatTmp;
delete [] b;
return Success;
}

//Example:
//min(z)=-(3X1+X2+3X3)
//2X1+X2+X3+X4=2
//X1+2X2+3X3+X5=5
//2X1+2X2+X3+X6=6
//Xi>=0,i=1,2,...,6

void main()
{
REAL A[]={
2,1,1,1,0,0,
1,2,3,0,1,0,
2,2,1,0,0,1
};
REAL B[]={
2,5,6
};
REAL C[]={
-3,-1,-3,0,0,0
};
REAL RESULT;
REAL x[3];
int Is[]={
1,2,3
};

if (!Line_Optimize(A,B,C,3,6,&RESULT,x,Is)){
cout<<"Calculate Wrong!"<<endl;
return;
}

for (int i=0; i<3; i++){
cout<<"The No."<<Is[i]<<" Root is:"<<x[i]<<endl;
}

cout << "The other root is Zero."<<endl;
cout << "The Minmize value is: "<<RESULT<<endl;
}

[lihaogao]

强 学习!!!

up!

[kwest]

想“小菜虎”学习，多收藏些经典的算法代码。

[kwest]

向“小菜虎”学习，多收藏些经典的算法代码。

[wangjinbao]

收了，能不能说说算法的思想！！！