Двойственная задача/симплекс-метод [c] - C (СИ)

Формулировка задачи:

Дана была прямая задача, необходимо найти решение, а после по полученному решению найти решение двойственной ей задачи. при данной системе ограничений: Решить прямую задачу симплекс-методом получилось, вывод двойственной задачи, задачи в канонической форме тоже, но что делать дальше я не знаю, сказано, что решение двойственной задачи может быть найдено из формулы: где D - матрица, составленная из векторов, входящих в последний базис, а - вектор базисных переменных(последний получившийся в ходе преобразований) Проблема-то, собственно, в том, что я не знаю, как получить обратную матрицу D, кроме того верно ее составив. И не помню, как верно перемножить вектор на матрицу. И код не является универсальным, написан для конкретно этой задачи, и я не знаю, как теперь все пристроить.

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#define WIDTH 5
#define HIGH 3
 
void printCanonic(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2);
void printTask(const int arr1[], const double mat[][WIDTH], const double vect[], const int len); // вывод задачи
void inverseMatrix(double **, double **, int);
void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2); // вывод двойствеоон задачи
double simplex(double A[HIGH][WIDTH], double A0[]); // 3-5 столбец - базис, 0 - вектор свободных членов
void printSimplex(int [], int[], int [], int[], double[], double[], double[][HIGH], int, double); //вывод симлекс-таблицы
int nminpos(double [], int); // минимальный отрицательный, позиция, -1 в случае отсутвия отриц.
int pminpos(double [], double[], int); // позиция минимального неотрицательго частного 
double vectorMultiple(int [], double [], int); //скалярное
void swap(int *, int *);
 
const int f[5] = { 2, 3 }; // последние 3 для введенных переменных
 
int main()
{
    double A[HIGH][WIDTH] = { {-3,2,1,0,0}, { 1, -4, 0, 1, 0}, { 1, -1, 0, 0, 1} };
    double b[HIGH] = { 6, 2, 5 };
 
    printTask(f, A, b, 2);
    printCanonic(f,A,b,2,HIGH);
    printDuality(f,A,b,2,HIGH);
 
    simplex(A, b);
    return 0;
 
}
 
void printCanonic(const int arr1[], const double mat[][5], const double vect[], const int len1, const int len2)
{
    int i, j;
    printf("##################################\n");
    printf("#     Task in canonical form     #\n");
    printf("##################################\n\n");
 
    printf("f(x) = ");
    for (i = 0; i < len1; i++)
        printf("%c%dx%d", (i>0 && arr1[i]>0) ? '+' : ' ' ,arr1[i], i+1);
    printf("-> max\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < len2; i++) {
        for (j = 0; j < 5; j++)
            printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1);
        printf(" = %.lf\n", vect[i]);
    }
    putchar('\n');
}
 
void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2)
{
    int i, j;
    printf("##################################\n");
    printf("#           Twofold task         #\n");
    printf("##################################\n\n");
 
    printf("g(y) = ");
    for (i = 0; i < len2; i++)
        printf("%c%.lfy%d", (i>0 && vect[i]>0) ? '+' : ' ' ,vect[i], i+1);
    printf("-> min\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < len1; i++) {
        for (j = 0; j < len2; j++)
            printf("%c%.lfx%d", (j>0 && mat[j][i]>=0) ? '+' : ' ' ,mat[j][i], j+1);
        printf(" >= %d\n", arr1[i]);
    }
    putchar('\n');
}
 
void printTask(const int arr[], const double mat[][WIDTH], const double vect[], const int len1)
{
    int i, j;
    printf("##################################\n");
    printf("#           Straight task        #\n");
    printf("##################################\n\n");
 
    printf("f(x) = ");
    for (i = 0; i < len1; i++)
        printf("%c%dx%d", (i>0 && arr[i]>0 && arr[i] != 1) ? '+' : ' ' ,arr[i], i+1);
    printf("-> max\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < HIGH; i++) {
        for (j = 0; j < len1; j++)
            printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1);
        printf(" <= %.lf\n", vect[i]);
    }
    putchar('\n');
}
 
double simplex(double A[HIGH][WIDTH], double A0[])
{
    int Cbasis[HIGH] = { 0, 0, 0 };
    int basis[HIGH] = {3, 4, 5};
    int notbase[2] = {1,2};
    int Cj[2] = { f[0], f[1] };
    double fstr[2]/*, x[WIDTH] = { 0, 0, A0[0], A0[1], A0[2]}*/;
    double matr[WIDTH][HIGH], commonmatr[WIDTH][HIGH], Q, **D, **D1;
    int r, s, i, j, counter = 0;
 
    for(i=0; i<WIDTH; i++)
        for (j=0; j<HIGH; j++)
            matr[i][j] = A[j][i];
 
    for (i=0; i<2;i++)
        fstr[i] = vectorMultiple(Cbasis, matr[i], HIGH) - f[i];
    Q = vectorMultiple(Cbasis, A0, HIGH);
    
    printf("\nSimplex-matrix:\n");
    printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q);
 
    while ((r = nminpos(fstr, 2)) != -1) {
        counter++;
        if ((s = pminpos(matr[r], A0, HIGH)) == -1)
            return -1;
        
        for(i=0; i<WIDTH; i++)
            for (j=0; j<HIGH; j++)
                commonmatr[i][j] = matr[i][j];
 
        matr[r][s] = 1/commonmatr[r][s];
 
        for (i = 0; i < HIGH; i++)
            if (i!=s) matr[r][i] = -commonmatr[r][i] / commonmatr[r][s];
 
        for (i = 0; i < 2; i++)
            if (i!=r) matr[i][s] = commonmatr[i][s] / commonmatr[r][s];
 
        for (i = 0; i < HIGH; i++)
            for (j=0; j<2; j++)
                if ((i != r) && (j!=s))
                    matr[i][j] = (commonmatr[i][j] * commonmatr[r][s] - commonmatr[r][j] * commonmatr[i][s])/commonmatr[r][s];
        
        for (i = 0; i < HIGH; i++) {
            if (i!=s)
                A0[i] = (A0[i] * commonmatr[r][s] - commonmatr[r][i] * A0[s])/commonmatr[r][s]; }
        A0[s] /= commonmatr[r][s];
 
        //f-строка
        for (i = 0; i < 2; i++) {
            if (i!=r)
                fstr[i] = (fstr[i] * commonmatr[r][s] - commonmatr[i][s] * f[r])/commonmatr[r][s]; }
        fstr[r] /= (-commonmatr[r][s]);
 
        swap(&Cbasis[s], &Cj[r]);
        swap(&basis[s], &notbase[r]);
 
        Q = vectorMultiple(Cbasis, A0, HIGH);
        printf("\nSimplex-matrix after %d iteration:\n", counter);
        printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q);
    }
    return Q;
}
 
void swap(int *x, int *y)
{
    int temp = *x;
    *x = *y;
    *y = temp;
}
 
int pminpos(double A[], double A0[], int length)
{
    int i, minp = -1;
    double min;
    min = DBL_MAX;
    for (i = 0; i < length; i++) {
        if ((A[i] > 0) && (A0[i]/A[i] < min)) {
            min = A0[i]/A[i];
            minp = i; }
    }
    return minp;
}
 
int nminpos(double A[], int length)
{
    int minp = -1, i;
    double min = 0;
 
    for (i = 0; i < length; i++) {
        if (A[i] < min) {
            min = A[i];
            minp = i; }
    }
    return minp;
}
 
void printSimplex(int Cbas[], int bas[], int cj[], int notb[], double fs[], double a0[], double matrix[][HIGH], int clmn, double Q)
{
    int i;
    printf("_____|__Cj__|__%-5d|__%-5d|_____\n", cj[0], cj[1]);
    printf("Basis|______|_x%-5d|_x%-5d|__A0_\n", notb[0], notb[1]);
    for (i = 0; i < HIGH; i++)
        printf("__%-3d|_x%-4d|%7.2lf|%7.2lf|%5.2lf\n", Cbas[i], bas[i], matrix[0][i], matrix[1][i], a0[i]);
    printf("_____|___f__|%7.2lf|%7.2lf|%5.2lf\n", fs[0], fs[1], Q);
    printf("_____|______|_DELTA1|_DELTA2|__Q__\n");
}
 
double vectorMultiple(int a1[], double a2[], int length)
{
    int i;
    double sum = 0;
 
    for (i = 0; i < length; i++)
        sum += a1[i]*a2[i];
    
    return sum;
}

В общем я совсем запутался, вполне возможно, забыл что-то ещё сделать или сделал не так.

Решение задачи: «Двойственная задача/симплекс-метод [c]»

textual
                            Листинг программы 
                        

                            #include <stdio.h>
#include <stdlib.h>
#include <float.h>
#define WIDTH 5
#define HIGH 3
 
void printCanonic(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2);
void printTask(const int arr1[], const double mat[][WIDTH], const double vect[], const int len);
void inverseMatrix(double **, int);
void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2);
double simplex(double A[HIGH][WIDTH], double A0[]); // 3-5 базис, 0 - вектор свободных членов
void printSimplex(int [], int[], int [], int[], double[], double[], double[][HIGH], int, double);
int nminpos(double [], int);
int pminpos(double [], double[], int);
double vectorMultiple(int [], double [], int);
void swap(int *, int *);
void matrXvect(double **, double *, int);
 
const int f[5] = { 2, 3 };
 
 
int main()
{
    double A[HIGH][WIDTH] = { {-3,2,1,0,0}, { 1, -4, 0, 1, 0}, { 1, -1, 0, 0, 1} };
    double b[HIGH] = { 6, 2, 5 }, result;
 
    printTask(f, A, b, 2);
    printCanonic(f,A,b,2,HIGH);
    printDuality(f,A,b,2,HIGH);
 
    result = simplex(A, b);
    if (result!=-1)
        printf("\nMaximim is %.2lf\nResult is correct.\n", result);
    else
        printf("\nInclorrect result/undecidability problem\n");
    return 0;
 
}
 
void printCanonic(const int arr1[], const double mat[][5], const double vect[], const int len1, const int len2)
{
    int i, j;
    printf("##################################\n");
    printf("#     Task in canonical form     #\n");
    printf("##################################\n\n");
 
    printf("f(x) = ");
    for (i = 0; i < len1; i++)
        printf("%c%dx%d", (i>0 && arr1[i]>0) ? '+' : ' ' ,arr1[i], i+1);
    printf("-> max\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < len2; i++) {
        for (j = 0; j < 5; j++)
            printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1);
        printf(" = %.lf\n", vect[i]);
    }
    putchar('\n');
}
 
void printDuality(const int arr1[], const double mat[][WIDTH], const double vect[], const int len1, const int len2)
{
    int i, j;
    printf("##################################\n");
    printf("#           Twofold task         #\n");
    printf("##################################\n\n");
 
    printf("g(y) = ");
    for (i = 0; i < len2; i++)
        printf("%c%.lfy%d", (i>0 && vect[i]>0) ? '+' : ' ' ,vect[i], i+1);
    printf("-> min\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < len1; i++) {
        for (j = 0; j < len2; j++)
            printf("%c%.lfx%d", (j>0 && mat[j][i]>=0) ? '+' : ' ' ,mat[j][i], j+1);
        printf(" >= %d\n", arr1[i]);
    }
    putchar('\n');
}
 
void printTask(const int arr[], const double mat[][WIDTH], const double vect[], const int len1)
{
    int i, j;
    printf("##################################\n");
    printf("#           Straight task        #\n");
    printf("##################################\n\n");
 
    printf("f(x) = ");
    for (i = 0; i < len1; i++)
        printf("%c%dx%d", (i>0 && arr[i]>0 && arr[i] != 1) ? '+' : ' ' ,arr[i], i+1);
    printf("-> max\n");
 
    printf("System of restrictions: \n");
    for (i = 0; i < HIGH; i++) {
        for (j = 0; j < len1; j++)
            printf("%c%.lfx%d", (j>0 && mat[i][j]>=0) ? '+' : ' ' ,mat[i][j], j+1);
        printf(" <= %.lf\n", vect[i]);
    }
    putchar('\n');
}
 
double simplex(double A[HIGH][WIDTH], double A0[])
{
    int Cbasis[HIGH] = { 0, 0, 0 };
    int basis[HIGH] = {3, 4, 5};
    int notbase[2] = {1,2};
    int Cj[2] = { f[0], f[1] };
    double fstr[2]/*, x[WIDTH] = { 0, 0, A0[0], A0[1], A0[2]}*/;
    double matr[WIDTH][HIGH], commonmatr[WIDTH][HIGH], Q, Qy, D[3][3], *x[3], basdouble[HIGH], commonA0[HIGH];
    int r, s, i, j, counter = 0;
 
    for (j=0; j<HIGH; j++)
        commonA0[j] = A0[j];
 
    for(i=0; i<WIDTH; i++)
        for (j=0; j<HIGH; j++)
            matr[i][j] = A[j][i];
 
    for (i=0; i<2;i++)
        fstr[i] = vectorMultiple(Cbasis, matr[i], HIGH) - f[i];
    Q = vectorMultiple(Cbasis, A0, HIGH);
    
    printf("\nSimplex-matrix:\n");
    printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q);
 
    while ((r = nminpos(fstr, 2)) != -1) {
        counter++;
        if ((s = pminpos(matr[r], A0, HIGH)) == -1)
            return -1;
        
        for(i=0; i<WIDTH; i++)
            for (j=0; j<HIGH; j++)
                commonmatr[i][j] = matr[i][j];
 
        matr[r][s] = 1/commonmatr[r][s];
 
        for (i = 0; i < HIGH; i++)
            if (i!=s) matr[r][i] = -commonmatr[r][i] / commonmatr[r][s];
 
        for (i = 0; i < 2; i++)
            if (i!=r) matr[i][s] = commonmatr[i][s] / commonmatr[r][s];
 
        for (i = 0; i < HIGH; i++)
            for (j=0; j<2; j++)
                if ((i != r) && (j!=s))
                    matr[i][j] = (commonmatr[i][j] * commonmatr[r][s] - commonmatr[r][j] * commonmatr[i][s])/commonmatr[r][s];
        
        for (i = 0; i < HIGH; i++) {
            if (i!=s)
                A0[i] = (A0[i] * commonmatr[r][s] - commonmatr[r][i] * A0[s])/commonmatr[r][s]; }
        A0[s] /= commonmatr[r][s];
 
        //f-строка
        for (i = 0; i < 2; i++) {
            if (i!=r)
                fstr[i] = (fstr[i] * commonmatr[r][s] - commonmatr[i][s] * f[r])/commonmatr[r][s]; }
        fstr[r] /= (-commonmatr[r][s]);
 
        swap(&Cbasis[s], &Cj[r]);
        swap(&basis[s], ¬base[r]);
 
        Q = vectorMultiple(Cbasis, A0, HIGH);
        printf("\nSimplex-matrix after %d iteration:\n", counter);
        printSimplex(Cbasis, basis, Cj, notbase, fstr, A0, matr, 2, Q);
    }
 
    for (i=0; i<3; i++)
        for (j=0;j<3;j++)
            D[i][j] = A[j][basis[i]-1];
 
    printf("\nMatrix D of last basis columns:\n");
 
    for (i=0; i<3; i++) {
        x[i] = D[i];
        for (j=0;j<3;j++)
        printf("%.2lf ", x[i][j] = D[i][j]);
        putchar('\n');}
 
    inverseMatrix(x, 3);
 
    printf("\nInvert matrix D^-1:\n");
 
    for (i=0; i<3; i++) {
        for (j=0;j<3;j++)
            printf("%.2lf ", x[i][j]);
        putchar('\n');}
 
 
    for (i=0; i<3; i++)
        basdouble[i] = Cbasis[i];
 
    matrXvect(x, basdouble, 3);
 
    printf("\nVector y = (");
    for (i=0; i<3; i++)
        printf("%.2lf, ",basdouble[i]);
    printf(")\n");
 
    printf("\nValue of duality funclion:\n");
    printf("g(y) = ");
    Qy = 0;
    for (i = 0; i < HIGH; i++) {
        Qy+=commonA0[i]*basdouble[i];
        printf("%c%.lfy%d", (i>0 && commonA0[i]>0) ? '+' : ' ' ,commonA0[i], i+1); }
    printf(" = %.2lf\n", Qy);
 
    if (Q!=Qy)
        return -1;
 
    return Q;
}
 
void swap(int *x, int *y)
{
    int temp = *x;
    *x = *y;
    *y = temp;
}
 
int pminpos(double A[], double A0[], int length)
{
    int i, minp = -1;
    double min;
    min = DBL_MAX;
    for (i = 0; i < length; i++) {
        if ((A[i] > 0) && (A0[i]/A[i] < min)) {
            min = A0[i]/A[i];
            minp = i; }
    }
    return minp;
}
 
int nminpos(double A[], int length)
{
    int minp = -1, i;
    double min = 0;
 
    for (i = 0; i < length; i++) {
        if (A[i] < min) {
            min = A[i];
            minp = i; }
    }
    return minp;
}
 
void printSimplex(int Cbas[], int bas[], int cj[], int notb[], double fs[], double a0[], double matrix[][HIGH], int clmn, double Q)
{
    int i;
    printf("_____|__Cj__|__%-5d|__%-5d|_____\n", cj[0], cj[1]);
    printf("Basis|<br>|_x%-5d|_x%-5d|__A0_\n", notb[0], notb[1]);
    for (i = 0; i < HIGH; i++)
        printf("__%-3d|_x%-4d|%7.2lf|%7.2lf|%5.2lf\n", Cbas[i], bas[i], matrix[0][i], matrix[1][i], a0[i]);
    printf("_____|___f__|%7.2lf|%7.2lf|%5.2lf\n", fs[0], fs[1], Q);
    printf("_____|<br>|_DELTA1|_DELTA2|__Q__\n");
}
 
double vectorMultiple(int a1[], double a2[], int length)
{
    int i;
    double sum = 0;
 
    for (i = 0; i < length; i++)
        sum += a1[i]*a2[i];
    
    return sum;
}
 
void inverseMatrix(double **A, int n)
{
 
    ///Метод гауссса-йордана
    int i,j,k,m;
    double akk, aik, r, E[HIGH][HIGH] = { {1,0,0}, {0,1,0}, {0,0,1} };
    m = 0;
 
    for (k=0; k <= n-1; k++)
    {
        if ((A[k][k]) == 0) {
            m = k +1;
            while (m<=n && A[m][k]==0)
                m++;
            if (m > n)// {
                printf("Matrix is degenerate\n");
                //return -1 }
            for (i = 0; i<n; i++) {
                r = A[m][i];
                A[m][i] = A[k][i];
                A[k][i] = r;
                
                r = E[m][i];
                E[m][i] = E[k][i];
                E[k][i] = r; }
        }
        akk = A[k][k];
        //A[k][k]/=akk;
        E[k][k]/=akk;
        for (i = 0; i < n; i++)
            if (i!=k) {
                aik = A[i][k]/akk;
                for (j = k; j < n; j++) {
                    A[i][j] = A[i][j] - aik*A[k][j];
                    E[i][j] = (E[i][j] - aik*E[k][j])/A[i][i]; }
            };
 
    };
 
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            A[i][j] = E[i][j];
}
 
void matrXvect(double **A, double *c, int n)
{
    double ii;
    int i, j;
 
    for (i = 0; i < n; i++) {
        ii = 0;
        for (j = 0; j < n; j++)
            ii += A[j][i]*c[i];
        c[i] = ii; }
}
                        