シンプレックス・タブロー

シンプレックス・タブローで線形計画を解くコードです。
自分の理解のため作りました。見れ分かると思いますが、適当です。
制約条件を<=にする場合は適当に補助変数を足してください。
教科書に載っているような3×3程度の問題ならこのコードで解けますが、
大規模（１００００列とか）になると掃き出していくうちに数値誤差が累積するのか
縮退・巡回しまくって収束しないです。==0 の条件を緩めるとか、Bland則・ピボットを
ランダム選択するなどをしないと駄目なようです。
でも、そうやるのも大変なのではっきりいって、線形計画パッケージのlp_solveを使ったほうがいいです。
lp_solveは拡張シンプレックス法（良く分からないが）を使って、整数計画問題にも対応していて
分枝限定法を使うらしいです。lp_solveは初めて使いましたが、簡単だし、ちゃんと答えもでるし、
優れものです。大規模線形計画コードは自作できそうにないので助かります。

// 目的関数：z = c^T x, 制約条件：Ax = b　をタブロー法で解く
// A:M行N列
double solve_simplex_tableau(double *A, double *b, double * c, double *x, int N, int M)
{
    double z0;
    double *tmp;
 
    tmp = new double[(M + 1) * (N + 1)];
    memset(tmp, 0, sizeof(double) * (M + 1) * (N + 1));
    int i, j;
    for(i = 0;i < M;i++)
    {
        for(j = 0;j < N;j++)
        {
            tmp[(i+1) * (N+1) + j] = A[i * N + j];
        }
        tmp[(i+1) * (N+1) + j] = b[i];
    }
    // 最上列
    for(j = 0;j < N;j++)
        tmp[j] = -c[j];
 
    int it;
    for(it = 0;it < 1000;it++){
        // 最上列から最小要素となる列を選択
        double min_z, min_c;
        int    min_i, min_j;
        min_z = 0.0;
        for(j = 0;j < N;j++)
        {
            if(min_z > tmp[j])
            {
                min_z = tmp[j];
                min_j = j;
            }
        }
        // 最上列に負の要素がなくなれば修了
        if(min_z == 0)
        {
            z0 = tmp[N];
            for(i = 0;i < N;i++)
                x[i] = tmp[(i+1)*(N+1)+N];
            goto END_LOOP;
        }
        // 最小列の要素で定数項を割り、定数項が最小となる行を探す
        min_c = 100000.0;
        for(i = 0;i < M;i++)
        {
            double p = tmp[(i+1) * (N+1) + N] / tmp[(i+1) * (N+1) + min_j];
            if(p < min_c)
            {
                min_c = p;
                min_i = (i+1);
            }
        }
        // (min_i, min_j)をピボットとして掃き出し実行
        double a;
        a = tmp[min_i * (N+1) + min_j];
        for(j = 0;j <= N;j++)
            tmp[min_i * (N+1) + j] /= a;
        for(i = 0;i <= M;i++)
        {
            if(i == min_i)
                continue;
            a = tmp[i * (N+1) + min_j];
            for(j = 0;j <= N;j++)
                tmp[i * (N+1) + j] -= a * tmp[min_i * (N+1) + j];
        };
    };
    
END_LOOP:
    delete tmp;
 
    return z0;
}