feasible method is stucked in infinite loop

Question

I have run the following simplex method which is written in C++ .The code is as follows :

#include<bits/stdc++.h>
#include<cstdio>
using namespace std;

#define maxm 500
#define maxn 500
double inf = 1e100;
double eps = 1e-13;

int row,col;
double A[maxm][maxn];
double B[maxn];
///////////////////////////////////////////////////////////////////////////////////////////
// Simon Lo's
// Simplex algorithm on augmented matrix a of dimension (m+1)x(n+1)
// returns 1 if feasible, 0 if not feasible, -1 if unbounded
// returns solution in b[] in original var order, max(f) in ret
// form: maximize sum_j(a_mj*x_j)-a_mn s.t. sum_j(a_ij*x_j)<=a_in
// in standard form.
// To convert into standard form:
// 1. if exists equality constraint, then replace by both >= and <=
// 2. if variable x doesn't have nonnegativity constraint, then replace by
// difference of 2 variables like x1-x2, where x1>=0, x2>=0
// 3. for a>=b constraints, convert to -a<=-b
// note: watch out for -0.0 in the solution, algorithm may cycle
// eps = 1e-7 may give wrong answer, 1e-10 is better



void pivot(int m, int n, double a[maxm][maxn], int B[maxm], int N[maxn], int r, int c) {
        int i, j;
        swap(N[c], B[r]);
        a[r][c]=1/a[r][c];
        for (j=0; j<=n; j++)if (j!=c) a[r][j]*=a[r][c];
        for (i=0; i<=m; i++)if (i!=r) {
                for (j=0; j<=n; j++)if (j!=c)
                        a[i][j]-=a[i][c]*a[r][j];
                a[i][c] = -a[i][c]*a[r][c];
        }
}
int feasible(int m, int n, double a[maxm][maxn], int B[maxm], int N[maxn]) {
        int r, c, i; double p, v;
        while (1) {
                for (p=inf, i=0; i<m; i++) if (a[i][n]<p) p=a[r=i][n];
                if (p>-eps) return 1;
                for (p=0, i=0; i<n; i++) if (a[r][i]<p) p=a[r][c=i];
                if (p>-eps) return 0;
                cout<<"Sultan"<<endl;
                p = a[r][n]/a[r][c];
                for (i=r+1; i<m; i++) if (a[i][c]>eps) {
                        v = a[i][n]/a[i][c];
                        if (v<p) r=i, p=v;
                }
                pivot(m, n, a, B, N, r, c);
        }
}
int simplex(int m, int n, double a[maxm][maxn], double b[maxn], double& ret)
 {
        int B[maxm], N[maxn], r, c, i; double p, v;
        for (i=0; i<n; i++) N[i]=i;
        for (i=0; i<m; i++) B[i]=n+i;
        if (!feasible(m, n, a, B, N)) return 0;
        while (1) {
                for (p=0, i=0; i<n; i++) if (a[m][i]>p)
                        p=a[m][c=i];
                if (p<eps) {
                        for (i=0; i<n; i++) if (N[i]<n)
                                b[N[i]]=0;
                        for (i=0; i<m; i++) if (B[i]<n)
                                b[B[i]]=a[i][n];
                        ret = -a[m][n];
                        return 1;
                }
                for (p=inf, i=0; i<m; i++) if (a[i][c]>eps) {
                        v = a[i][n]/a[i][c];
                        if (v<p) p=v, r=i;
                }
                if (p==inf) return -1;
                pivot(m, n, a, B, N, r, c);
        }
}
//////////////////////////////////////////////////////////////////////////////////////////////

void read_file()
{
    freopen("Dimen.txt","r",stdin);

    scanf("%d",&row);
    scanf("%d",&col);

    row = 30 ;
    col = 100;

    cout << row << " "<<col<<endl;

    freopen("A1.txt","r",stdin);
    for(int i=0;i<row;i++)
       for(int j=0;j<col;j++)
         scanf("%lf",&A[i][j]);
    cout<"Completed A";

    freopen("B1.txt","r",stdin);
    for(int i=0;i<row;i++){
       scanf("%lf",&A[i][col]);
    }
    cout<"Completed B";

    freopen("F1.txt","r",stdin);
    for(int j=0;j<col;j++){
        scanf("%lf",&A[row][j]);
        //B[i]=-B[i];
    }
    cout<"Completed F";


}

int main()
{
    read_file();
    double value,opt=0;
    int flag = simplex(row,col, A,B, value);
    freopen("opt.txt","r",stdin);
    scanf("%lf",&opt);
    if(flag != -1)
    {
        cout<<"The result is "<<(value+opt);
    }
    else
    {
        cout<<"The result is infeasible";
    }

    return 0;
}

But this code runs in infinite loop . The feasible method stucks in infinite loop . How can I solve this error ? Please help me .

Answer 1

Lets start by fixing the unclear bits of the code:

#include<bits/stdc++.h>  

You should not use any includes that starts with bits/ include the standard header that will include it.

#include<cstdio>
using namespace std;

Using std is bad practice, it might confuse the compiler in large projects and select a function/template/whatever that you did not mean to.

#define maxm 500
#define maxn 500

As minimum change these to

#define maxm (500)
#define maxn (500)

So macro substitution doesn't get confused. Better use c++11 and use

constexpr int maxm = 500; 
constexpr int maxn = 500; 

double inf = 1e100;
double eps = 1e-13;

There are standard ways of finding the correct values for inf and eps

double inf = std::numeric_limits<double>::max();

From you code I presume you really meant max double. For eps I wonder if you need

double eps = std::numeric_limits<double>::epsilon();

or

double eps = std::numeric_limits<double>::round_error();