Welcome, Guest! Sign Up RSS

Clever Space

Friday, 11.22.2024
Main » 2014 » August » 3 » 一份较详细的单纯形模板
3:28 PM
一份较详细的单纯形模板

float a[MAX][MAX], b[MAX], c[MAX], x[MAX], v;
int n, m;

void pivot(bool flag_row[], bool flag_column[], float c[MAX], int l, int e)
{
 //计算原等式变形系数
 b[e] = b[l] / a[l][e];

 for(int j = 0; j <= m+n; ++j)
  if(j != e && flag_column[j]) a[e][j] = a[l][j] / a[l][e];
 
 a[e][l] = 1 / a[l][e];

 //更新限制等式
 for(int i = 0; i <= m+n; ++i)
 {
  if(i != l && flag_row[i]) 
  {
   b[i] -= (a[i][e]*b[e]);

   for(int j = 0; j <= m+n; ++j)
    if(j != e && flag_column[j]) a[i][j] -= (a[i][e]*a[e][j]);

   a[i][l] = -a[i][e]*a[e][l];
  }
 }

 //更新目标等式
 v += (c[e]*b[e]);

 for(int j = 0; j <= m+n; ++j)
  if(j != e && flag_column[j]) c[j] -= (c[e]*a[e][j]);

 c[l] = -c[e]*a[e][l];

 //更新集合标记
 flag_row[l] = false;
 flag_row[e] = true;

 flag_column[l] = true;
 flag_column[e] = false;
}

bool Run(bool flag_row[], bool flag_column[], float c[])
{
 int i, j;
 
 while(true)
 {
  for(j = 0; j <= m+n; ++j)
  {
   if(flag_column[j] && c[j] > ESP)
   {
    break;
   }
  }

  if(j == m+n+1) break;

  float data = (float)INF;
  int e = j, l;

  for(i = 0; i <= m+n; ++i)
  {
   if(a[i][e] > ESP && flag_row[i])
   {
    float tmp = b[i]/a[i][e];
    if(data > tmp)
    {
     data = tmp;
     l = i;
    }
   }
  }

  if(data == INF) return false;
  else pivot(flag_row, flag_column, c, l, e);
 }

 for(j = 0; j <= n+m; ++j)
 {
  if(flag_column[j])
   x[j] = 0;
  else
   x[j] = b[j];
 }

 return true;
}

bool Initialize(bool flag_row[], bool flag_column[])
{
 float tmp = (float)INF;
 float tmpc[MAX];
 int l;

 memset(tmpc, 0, sizeof(tmpc));

 for(int i = n; i <= m+n; ++i)
 {
  if(tmp > b[i])
  {
   l = i;
   tmp = b[i];
  }
 }

 v = 0.0f;

 for(int j = 1; j <= n; ++j) flag_column[j] = true;
 for(int i = 1; i <= m; ++i) flag_row[n+i] = true;

 if(tmp < 0.0f)
 {
  tmpc[0] = -1;
  flag_column[0] = true;

  for(int i = 1; i <= m; ++i)
   a[n+i][0] = -1;
  
  pivot(flag_row, flag_column, tmpc, l, 0);

  if(!Run(flag_row, flag_column, tmpc)) return false;

  if(v < 0.0f) return false;
  else v = 0.0f;

  flag_row[0] = flag_column[0] = false;

  for(int i = 0; i <= m+n; ++i)
  {
   if(c[i] != 0.0f && flag_row[i])
   {
    v += c[i]*b[i];

    for(int j = 0; j <= m+n; ++j)
    {
     if(flag_column[j]) c[j] -= a[i][j]*c[i];
    }

    c[i] = 0.0f;
   }
  }
 }

 return true;
}

 

//样例///////////////////////////////////////////

//m = n = 3;

//c[1] = 3; c[2] = 1; c[3] = 2;
//b[4] = 30; b[5] = 24; b[6] = 36;

//a[4][1] = 1; a[4][2] = 1; a[4][3] = 3;
//a[5][1] = 2; a[5][2] = 2; a[5][3] = 5;
//a[6][1] = 4; a[6][2] = 1; a[6][3] = 2;
////////////////////////////////////////////////

 

bool Simplex()  //n个未知列,m个约束等式
{
 bool flag_row[MAX], flag_column[MAX];
 v = 0.0f;

 memset(flag_row, 0, sizeof(flag_row));
 memset(flag_column, 0, sizeof(flag_column));

 if(!Initialize(flag_row, flag_column))
  return false;

 return Run(flag_row, flag_column, c);
}

Views: 600 | Added by: dhy0077 | Rating: 5.0/1
Total comments: 0
Only registered users can add comments.
[ Sign Up | Login ]