xianxingguihua.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]
{
           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;
}