运筹学上机实验 - 单纯形方法的两阶段法

//运行环境GCC6.4.0 C++11  非VC++
//实验一 单纯性方法
#include "cmath"
#include "cstdio"
#include "vector"
#include "algorithm"
#include "iostream"
using namespace std;

int N, M;
double **A;
int *mark;

bool Simplex(int ROW, int COL) {//两阶段法
    //mark标记进基变量
    mark = (int *)malloc(sizeof(int)*N);
    for (int i = 0; i < M; i++) mark[i] = M+i;
    //使得添加的变量的检验数为0
    for (int i = 2; i < N+2; i++) {
        for (int j = 0; j < COL; j++) {
            A[1][j] += A[i][j];
        }
    }
    //按照单纯形方法进行迭代
    double maxx = -1;
    int pos = 0;
    for (int i = 0; i < N+M; i++) {
        if (A[1][i] > maxx) {
            maxx = A[1][i]; pos = i;
        }
    }
    while (maxx > 0) {
        double minn = -1;
        int pos1 = 0;
        for (int i = 0; i < N; i++) {
            if (A[i+2][pos] > 0) {
                if (minn < 0) { pos1 = i+2; minn = A[i+2][M+N]/A[i+2][pos];}
                else if (A[i+2][M+N]/A[i+2][pos] < minn) {
                    pos1 = i+2;
                    minn = A[i+2][N+M]/A[i+2][pos];
                }
            }
        }
        if (minn == -1.0) return false;
        double s = A[pos1][pos];
        for (int i = 0; i < COL; i++) {
            A[pos1][i] /= s;
        }
        mark[pos1-2] = pos;
        for (int i = 0; i < ROW; i++) {
            if (i == pos1) continue;
            s = A[i][pos];
            for (int j = 0; j < COL; j++) {
                A[i][j] -= s*A[pos1][j];
            }
        }
        maxx = -1;
        pos = 0;
        for (int i = 0; i < N+M; i++) {
            if (A[1][i] > maxx) {
                maxx = A[1][i]; pos = i;
            }
        }
    }
    return true;
}
int main(int argc, char const *argv[])
{
    printf("请输入自变量的个数和方程组的个数:");
    scanf("%d%d", &M, &N);
    A = (double **)malloc(sizeof(double *)*(N+3)); //A为单纯形表
    for (int i = 0; i < N+3; i++) {
        *(A+i) = (double *)malloc(sizeof(double)*(M+N+2));
    }
    printf("请输入约束矩阵:\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < M; j++) {
            scanf("%lf", &A[i+2][j]);
        }
    }
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (i == j) A[i+2][j+M] = 1.0;
            else A[i+2][j+M] = 0.0;
        }
    }
    printf("请输入价值向量:\n");
    for (int i = 0; i < M; i++) {
        double t; scanf("%lf", &t);
        A[0][i] = -t;
    }
    for (int i = 0; i < N + M; i++) {
        if (i < M) A[1][i] = 0.0;
        else A[1][i] = -1.0;
    }
    printf("请输入右端向量:\n");
    A[0][N+M] = A[1][N+M] = 0.0;
    for (int i = 0; i < N; i++) {
        scanf("%lf", &A[i+2][N+M]);
    }
    //向量R中保存的是基本解向量对应的值
    int *R = (int *)malloc(sizeof(int)*(N+M));
    for (int i = 0; i < N+M; i++) R[i] = -1;
    if (Simplex(N+2, M+N+1)) {
        double g = 0.0;//辅助问题的g
        for (int i = 0; i < N; i++) {
            R[mark[i]] = i;
            if (mark[i] >= M) g += A[i][M+N];
        }
        // 如果min g != 0 方程无解
        if (g > 0) printf("该线性方程无解！\n");
        else {
            printf("利用单纯形方法得到的解为%.6lf\n", A[0][M+N]);
            printf("该线性规划的一个基本可行解为:\n");
            for (int i = 0; i < M; i++) {
                if (R[i] != -1) printf("%lf\n", A[R[i]+2][M+N]);
                else printf("%lf\n", 0.0);
            }
        }
    }
    else printf("该线性方程无解！\n");
    return 0;
}