diff options
Diffstat (limited to 'graphics/asymptote/base/simplex.asy')
-rw-r--r-- | graphics/asymptote/base/simplex.asy | 54 |
1 files changed, 29 insertions, 25 deletions
diff --git a/graphics/asymptote/base/simplex.asy b/graphics/asymptote/base/simplex.asy index c27b3236be..e6d0410b23 100644 --- a/graphics/asymptote/base/simplex.asy +++ b/graphics/asymptote/base/simplex.asy @@ -8,6 +8,7 @@ struct simplex { int case; real[] x; real cost; + bool dual=false; int m,n; int J; @@ -43,7 +44,7 @@ struct simplex { int iterate(real[][] 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 real[] Em=E[m]; for(J=1; J <= N; ++J) if(Em[J] < 0) break; @@ -67,7 +68,7 @@ struct simplex { real 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) @@ -82,8 +83,7 @@ struct simplex { int iterateDual(real[][] E, int N, int[] Bindices) { while(true) { - // Find first negative entry in zeroth (basic variable) column - real[] 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; @@ -92,23 +92,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; + } + + real[] Em=E[m]; + real[] EI=E[I]; int J=0; real t; for(int j=1; j <= N; ++j) { - real u=E[I][j]; + real u=EI[j]; if(u < -EpsilonA) { - t=-E[m][j]/u; + t=-Em[j]/u; J=j; break; } } for(int j=J+1; j <= N; ++j) { - real u=E[I][j]; + real u=EI[j]; if(u < -EpsilonA) { - real r=-E[m][j]/u; - if(r <= t && (r < t || j < J)) { + real r=-Em[j]/u; + if(r < t) { t=r; J=j; - } // Bland's rule: exiting variable has smallest minimizing index + } // Bland's rule: smallest minimizing subscript enters } } if(J == 0) @@ -125,9 +132,7 @@ struct simplex { // where A is an m x n matrix, x is a vector of n non-negative numbers, // 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(real[] c, real[][] A, real[] b, bool phase1=true, - bool dual=false) { - if(dual) phase1=false; + void operator init(real[] c, real[][] A, real[] b, bool phase1=true) { static real epsilon=sqrt(realEpsilon); real normA=norm(A); real epsilonA=100.0*realEpsilon*normA; @@ -324,11 +329,11 @@ struct simplex { real[] 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; @@ -339,19 +344,18 @@ 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; } } - operator init(concat(c,array(count,0.0)),a,b,phase1,dual); + if(dual) phase1=false; + operator init(concat(c,array(count,0.0)),a,b,phase1); if(case == OPTIMAL && count > 0) x.delete(n,n+count-1); |