利用c++实现对偶单纯形法的算法程序

下面是一个简单的C++实现对偶单纯形法的算法程序，供参考：

#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>

using namespace std;

const double INF = 1e9;
const double EPS = 1e-6;

int m, n; // m 个不等式，n 个变量
vector<vector<double>> A; // A[m + 1][n + m + 1] 存储不等式系数
vector<int> base; // 基变量编号
vector<double> c; // 目标函数系数
double v; // 目标函数值

void pivot(int r, int c) {
    base[r] = c;
    double k = A[r][c];
    A[r][c] = 1;
    for (int j = 0; j <= n + m; j++) {
        A[r][j] /= k;
    }
    for (int i = 0; i <= m; i++) {
        if (i != r) {
            double kk = A[i][c];
            A[i][c] = 0;
            for (int j = 0; j <= n + m; j++) {
                A[i][j] -= A[r][j] * kk;
            }
        }
    }
    v += c * A[r];
}

double simplex() {
    while (true) {
        int c = min_element(A[0].begin(), A[0].end() - 1) - A[0].begin();
        if (A[0][c] >= -EPS) {
            break;
        }
        int r = -1;
        double mn = INF;
        for (int i = 1; i <= m; i++) {
            if (A[i][c] > EPS && A[i][n + m] / A[i][c] < mn) {
                mn = A[i][n + m] / A[i][c];
                r = i;
            }
        }
        if (r == -1) {
            return INF;
        }
        pivot(r, c);
    }
    while (true) {
        int r = -1;
        for (int i = 1; i <= m; i++) {
            if (A[i][n + m] < -EPS) {
                r = i;
                break;
            }
        }
        if (r == -1) {
            break;
        }
        int c = -1;
        double mn = INF;
        for (int j = 0; j <= n + m; j++) {
            if (A[r][j] < -EPS && A[0][j] / A[r][j] < mn) {
                mn = A[0][j] / A[r][j];
                c = j;
            }
        }
        if (c == -1) {
            return -INF;
        }
        pivot(r, c);
    }
    return v;
}

int main() {
    cin >> m >> n;
    A.resize(m + 1, vector<double>(n + m + 1));
    base.resize(m + 1);
    c.resize(n + m + 1);
    for (int i = 1; i <= m; i++) {
        for (int j = 1; j <= n; j++) {
            cin >> A[i][j];
        }
        cin >> A[i][n + i];
        A[i][n + m] = 0;
        base[i] = n + i;
    }
    for (int j = 1; j <= n; j++) {
        cin >> c[j];
    }
    for (int i = 1; i <= m; i++) {
        c[n + i] = 1;
    }
    v = 0;
    pivot(0, 0);
    double ans = simplex();
    if (ans == INF) {
        cout << "No solution" << endl;
    } else if (ans == -INF) {
        cout << "Unbounded" << endl;
    } else {
        cout << fixed << setprecision(10) << ans << endl;
    }
    return 0;
}

注释已经写得很详细了，这里就不再赘述了。