<>帮knight代发他改进后的代码:</P><>/*
Name: LinearProgramming
Copyright: Knight-Studio
Author: Ran Ni
Date: 12-04-04 14:29
Description: The Definition of Class LinearProgramming
*/
#include <vector>
#include <math.h>
#include <algorithm>
#include "LP.h"</P><>using namespace std;</P><P>void Unit(vector<double> &x)
{
for(int i=0;i<(int)x.size();++i)
{
x<i>=0;
}
}
double operator * (vector<double> a,vector<double> b)
{
double data(0);
for(int i=0;i<(int)a.size();++i)
{
data+=a<i>*b<i>;
}
return data;
}</P><P>vector<vector<double> > Inverse(vector<vector<double> > a)
{
vector<vector<double> > b(a.size());
for(int i=0;i<(int)b.size();++i)
{
vector<double> c(a.size(),0);
c<i>=1;
b<i>=c;
}
for(int k=0;k<(int)b.size();++k)
{
int g=k;
for(int i=k;i<(int)b.size();++i)
{
if(fabs(a[k]<i>)>fabs(a[k][g]))
{
g=i;
}
}
if(g!=k)
{
for(int i=0;i<(int)b.size();++i)
{
double temp=a<i>[k];
a<i>[k]=a<i>[g];
a<i>[g]=temp;
temp=b<i>[k];
b<i>[k]=b<i>[g];
b<i>[g]=temp;
}
}
if(a[k][k]!=0)
{
double x=a[k][k];
for(int i=0;i<(int)b.size();++i)
{
a<i>[k]=(a<i>[k]/x);
b<i>[k]=(b<i>[k]/x);
}
for(int i=0;i<(int)b.size();++i)
{
if(i!=k)
{
double y=a[k]<i>;
for(int j=k;j<(int)b.size();++j)
{
a[j]<i>=a[j]<i>-a[j][k]*y/a[k][k];
}
for(int j=0;j<(int)b.size();++j)
{
b[j]<i>=b[j]<i>-b[j][k]*y/a[k][k];
}
}
}
}
}
return b;
}
vector<double> operator * (vector<vector<double> > a,vector<double> b)
{
vector<double> x(b.size(),0);
for(int i=0;i<(int)x.size();++i)
{
double data(0);
for(int j=0;j<(int)b.size();++j)
{
data+=b[j]*a[j]<i>;
}
x<i>=data;
}
return x;
}</P><P>vector<double> CalculateW(vector<double> c,vector<vector<double> > b)
{
vector<double> data(b.size(),0);
for(int i=0;i<(int)data.size();++i)
{
data<i>=c*b<i>;
}
return data;
}</P><P>vector<double> operator - (vector<double> a,vector<double> b)
{
for(int i=0;i<(int)a.size();++i)
{
a<i>-=b<i>;
}
return a;
}</P><P>int max(vector<double> c)
{
int k(0);
for(int i=0;i<(int)c.size();++i)
{
if(c<i>>c[k])
{
k=i;
}
}
return k;
}</P><P>int min(vector<double> c)
{
int k(0);
for(int i=0;i<(int)c.size();++i)
{
if(c<i><c[k])
{
k=i;
}
}
return k;
}</P><P>bool Judge(vector<double> a)
{
int i=0;
while((a<i><=0)&&(i<(int)a.size()))
{
i++;
}
if(i==(int)a.size())
{
return true;
}
else
{
return false;
}
}</P><P>bool Jugdeb(vector<double> b)
{
int i=0;
while((b<i>>=0)&&(i<(int)b.size()))
{
i++;
}
if(i==(int)b.size())
{
return true;
}
else
{
return false;
}
}</P><P>bool LinearProgramming::Run()
{
vector<vector<double> > AA(A);
Num=0;
for(int i=0;i<(int)b.size();++i)
{
vector<double> tem(b.size(),0);
tem<i>=1;
AA.push_back(tem);
}
vector<double> cc(c);
int k;
for(int i=0;i<(int)b.size();++i)
{
cc.push_back(0);
}
vector<int> Base(b.size(),0);
for(int i=0;i<(int)b.size();++i)
{
Base<i>=A.size()+i;
}
vector<double> cb(b.size(),0);
vector<vector<double> > Br;
vector<double> xb;
vector<double> xx(AA.size(),0);
double f;
vector<double> w;
vector<double> xt(AA.size()+1,0);
vector<double> dd(AA.size(),0);
vector<double> bb(b);
vector<vector<double> > B(b.size());
if(!Jugdeb(bb))
{
vector<double> xa(b.size(),-1);
AA.push_back(xa);
vector<double> ct(cc.size(),0);
ct.push_back(1);
int rt(min(bb));
for(int i=0;i<(int)AA.size();++i)
{
AA<i>[rt]=-AA<i>[rt];
}
bb[rt]=-bb[rt];
for(int i=0;i<(int)bb.size();++i)
{
if(i!=rt)
{
for(int j=0;j<(int)AA.size();++j)
{
AA[j]<i>-=(AA[j][rt]*(AA[AA.size()-1]<i>/AA[AA.size()-1][rt]));
}
bb<i>+=(bb[rt]);
}
}
Base[rt]=AA.size()-1;
++Num;
step3: vector<vector<double> > Bs;
for(int i=0;i<(int)Base.size();++i)
{
Bs.push_back(AA[Base<i>]);
}
for(int i=0;i<(int)bb.size();++i)
{
cb<i>=(ct[Base<i>]);
}
Br=Inverse(Bs);
xb=(Br*bb);
Unit(xt);
for(int i=0;i<(int)Base.size();++i)
{
xt[Base<i>]=xb<i>;
}
w=(CalculateW(cb,Br));
Unit(dd);
for(int i=0;i<(int)AA.size();++i)
{
if(find(Base.begin(),Base.end(),i)==Base.end())
{
dd<i>=(w*AA<i>)-ct<i>;
}
}
int ks(max(dd));
if(dd[ks]<=0)
{
goto step1;
}
else
{
vector<double> yk(Br*AA[ks]);
int r(0);
while(!(yk[r]>0))
{
++r;
}
double dddd(bb[r]/yk[r]);
for(int i=r;i<(int)yk.size();++i)
{
if((yk<i>>0)&&((bb<i>/yk<i>)<dddd))
{
r=i;
dddd=bb<i>/yk<i>;
}
}
Base[r]=ks;
++Num;
goto step3;
}
}
step1: if(find(Base.begin(),Base.end(),AA.size()-1)!=Base.end())
{
if(xt[AA.size()-1]!=0)
{
Feasible=false;
goto step2;
}
else
{
vector<int> :: iterator it=find(Base.begin(),Base.end(),(int)AA.size()-1);
int i=0;
while(Base<i>!=(*it))
{
++i;
}
int j=0;
while(!((AA[j]<i>!=0)&&(find(Base.begin(),Base.end(),j)==Base.end())))
{
++j;
}
*it=j;
}
}
for(int i=0;i<(int)B.size();++i)
{
B<i>=(AA[Base<i>]);
}
for(int i=0;i<(int)bb.size();++i)
{
cb<i>=(cc[Base<i>]);
}
Br=Inverse(B);
xb=(Br*bb);
Unit(xx);
for(int i=0;i<(int)Base.size();++i)
{
xx[Base<i>]=xb<i>;
}
f=cb*xb;
w=(CalculateW(cb,Br));
Unit(dd);
for(int i=0;i<(int)AA.size()-1;++i)
{
if(find(Base.begin(),Base.end(),i)==Base.end())
{
dd<i>=(w*AA<i>)-cc<i>;
}
}
k=max(dd);
if(dd[k]<=0)
{
this->Feasible=true;
this->Limited=true;
BestSolution=xx;
Minimum=f;
goto step2;
}
else
{
vector<double> yk(Br*AA[k]);
if(Judge(yk))
{
this->Feasible=false;
this->Limited=false;
goto step2;
}
else
{
int r(0);
while(!(yk[r]>0))
{
++r;
}
double dddd(bb[r]/yk[r]);
for(int i=r;i<(int)yk.size();++i)
{
if((yk<i>>0)&&((bb<i>/yk<i>)<dddd))
{
r=i;
dddd=bb<i>/yk<i>;
}
}
Base[r]=k;
++Num;
goto step1;
}
}
step2: if(this->Feasible==true)
{
return true;
}
else
{
return false;
}
}
</P> |