线性规划单纯形法源码.txt

线性规划单纯形法源码,线性规划单纯形法源码.

#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;  
}