summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/asymptote/rationalSimplex.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/asymptote/rationalSimplex.asy')
-rw-r--r--Master/texmf-dist/asymptote/rationalSimplex.asy99
1 files changed, 61 insertions, 38 deletions
diff --git a/Master/texmf-dist/asymptote/rationalSimplex.asy b/Master/texmf-dist/asymptote/rationalSimplex.asy
index 823d91c2d0b..3a376ba2ddb 100644
--- a/Master/texmf-dist/asymptote/rationalSimplex.asy
+++ b/Master/texmf-dist/asymptote/rationalSimplex.asy
@@ -1,8 +1,8 @@
// Rational simplex solver written by John C. Bowman and Pouria Ramazi, 2018.
import rational;
-void simplexStandard(rational[] c, rational[][] A, int[] s=new int[],
- rational[] b) {}
+void simplexInit(rational[] c, rational[][] A, int[] s=new int[],
+ rational[] b, int count) {}
void simplexTableau(rational[][] E, int[] Bindices, int I=-1, int J=-1) {}
void simplexPhase1(rational[] c, rational[][] A, rational[] b,
int[] Bindices) {}
@@ -35,7 +35,10 @@ struct simplex {
int case;
rational[] x;
+ rational[] xStandard;
rational cost;
+ rational[] d;
+ bool dual=false;
int m,n;
int J;
@@ -70,7 +73,7 @@ struct simplex {
int iterate(rational[][] E, int N, int[] Bindices) {
while(true) {
- // Find first negative entry in bottom (reduced cost) row
+ // Bland's rule: first negative entry in reduced cost (bottom) row enters
rational[] Em=E[m];
for(J=1; J <= N; ++J)
if(Em[J] < 0) break;
@@ -94,7 +97,7 @@ struct simplex {
rational r=E[i][0]/u;
if(r <= t && (r < t || Bindices[i] < Bindices[I])) {
t=r; I=i;
- } // Bland's rule: exiting variable has smallest minimizing index
+ } // Bland's rule: exiting variable has smallest minimizing subscript
}
}
if(I == -1)
@@ -111,8 +114,7 @@ struct simplex {
int iterateDual(rational[][] E, int N, int[] Bindices) {
while(true) {
- // Find first negative entry in zeroth (basic variable) column
- rational[] Em=E[m];
+ // Bland's rule: negative variable with smallest subscript exits
int I;
for(I=0; I < m; ++I) {
if(E[I][0] < 0) break;
@@ -121,23 +123,30 @@ struct simplex {
if(I == m)
break;
+ for(int i=I+1; i < m; ++i) {
+ if(E[i][0] < 0 && Bindices[i] < Bindices[I])
+ I=i;
+ }
+
+ rational[] Em=E[m];
+ rational[] EI=E[I];
int J=0;
rational t;
for(int j=1; j <= N; ++j) {
- rational u=E[I][j];
+ rational u=EI[j];
if(u < 0) {
- t=-E[m][j]/u;
+ t=-Em[j]/u;
J=j;
break;
}
}
for(int j=J+1; j <= N; ++j) {
- rational u=E[I][j];
+ rational u=EI[j];
if(u < 0) {
- rational r=-E[m][j]/u;
+ rational r=-Em[j]/u;
if(r <= t && (r < t || j < J)) {
t=r; J=j;
- } // Bland's rule: exiting variable has smallest minimizing index
+ } // Bland's rule: smallest minimizing subscript enters
}
}
if(J == 0)
@@ -157,8 +166,7 @@ struct simplex {
// b is a vector of length m, and c is a vector of length n.
// Can set phase1=false if the last m columns of A form the identity matrix.
void operator init(rational[] c, rational[][] A, rational[] b,
- bool phase1=true, bool dual=false) {
- if(dual) phase1=false;
+ bool phase1=true) {
// Phase 1
m=A.length;
if(m == 0) {case=INFEASIBLE; return;}
@@ -243,7 +251,7 @@ struct simplex {
simplexPhase1(c,A,b,Bindices);
iterate(E,n+k,Bindices);
-
+
if(Em[0] != 0) {
simplexTableau(E,Bindices);
case=INFEASIBLE;
@@ -278,7 +286,7 @@ struct simplex {
for(int i=0; i < m; ++i) {
int k=Bindices[i];
if(k > n) continue;
- Bindices[ip]=k;
+ Bindices[ip]=k;
cB[ip]=c[k-1];
rational[] Dip=D[ip];
rational[] Ei=E[i];
@@ -319,14 +327,27 @@ struct simplex {
case=(dual ? iterateDual : iterate)(D,n,Bindices);
simplexTableau(D,Bindices);
- if(case != OPTIMAL)
- return;
- for(int j=0; j < n; ++j)
- x[j]=0;
+ if(case != INFEASIBLE) {
+ x=new rational[n];
+ for(int j=0; j < n; ++j)
+ x[j]=0;
- for(int k=0; k < m; ++k)
- x[Bindices[k]-1]=D[k][0];
+ for(int k=0; k < m; ++k)
+ x[Bindices[k]-1]=D[k][0];
+ }
+
+ if(case == UNBOUNDED) {
+ d=new rational[n];
+ for(int j=0; j < n; ++j)
+ d[j]=0;
+ d[J-1]=1;
+ for(int k=0; k < m; ++k)
+ d[Bindices[k]-1]=-D[k][J];
+ }
+
+ if(case != OPTIMAL)
+ return;
cost=-Dm[0];
}
@@ -353,21 +374,21 @@ struct simplex {
ai[j]=Ai[j];
}
}
-
+
int k=0;
bool phase1=false;
- bool dual=count == m && all(c >= 0);
+ dual=count == m && all(c >= 0);
for(int i=0; i < m; ++i) {
rational[] ai=a[i];
for(int j=0; j < k; ++j)
ai[n+j]=0;
+ int si=s[i];
if(k < count)
- ai[n+k]=-s[i];
+ ai[n+k]=-si;
for(int j=k+1; j < count; ++j)
ai[n+j]=0;
- int si=s[i];
if(si == 0) phase1=true;
else {
++k;
@@ -378,23 +399,25 @@ struct simplex {
for(int j=0; j < n+count; ++j)
ai[j]=-ai[j];
}
- } else if(si*bi > 0) {
- if(dual && si == 1) {
- b[i]=-bi;
- s[i]=-1;
- for(int j=0; j < n+count; ++j)
- ai[j]=-ai[j];
- } else
- phase1=true;
- }
+ } else if(dual && si == 1) {
+ b[i]=-bi;
+ s[i]=-1;
+ for(int j=0; j < n+count; ++j)
+ ai[j]=-ai[j];
+ } else if(si*bi > 0)
+ phase1=true;
}
}
+ if(dual) phase1=false;
rational[] C=concat(c,array(count,rational(0)));
- if(count > 0) simplexStandard(C,a,b);
- operator init(C,a,b,phase1,dual);
+ simplexInit(C,a,b,count);
+ operator init(C,a,b,phase1);
- if(case == OPTIMAL && count > 0)
- x.delete(n,n+count-1);
+ if(case != INFEASIBLE) {
+ xStandard=copy(x);
+ if(count > 0)
+ x.delete(n,n+count-1);
+ }
}
}