From dbd941ed19b558edd09219a372b2b0832957b283 Mon Sep 17 00:00:00 2001 From: Karl Berry Date: Wed, 24 Feb 2021 18:33:44 +0000 Subject: asymptote 2.69 support files git-svn-id: svn://tug.org/texlive/trunk@57876 c570f23f-e606-0410-a88d-b1316a301751 --- Master/texmf-dist/asymptote/rationalSimplex.asy | 99 +++++++++++++++---------- 1 file changed, 61 insertions(+), 38 deletions(-) (limited to 'Master/texmf-dist/asymptote/rationalSimplex.asy') 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); + } } } -- cgit v1.2.3