1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
|
const Basic = 1
const AtLower = 2
const AtUpper = 3
const dualTol = 1e-7
# Adapted from https://github.com/mlubin/SimplexBenchmarks.
# See that project for a suite of cross-language benchmarks.
# This code performs a two-pass "Harris" ratio test used within the
# dual simplex algorithm for linear programming.
# See the thesis of Achim Koberstein,
# "The dual simplex method, techniques for a fast and stable implementation"
# 2005, Section 6.2.2.2.
# Data taken from sample greenbea iteration
function doTwoPassRatioTest()
n = 8000
candidates = zeros(Int,n) # don't count allocation time, assume reuse
thetaMax = 1e25
pivotTol = 1e-7
red1 = [144.96711464551225,1.1799262694258228,15.258941028685781,0.0,74.21362490887168,8.771366223499031,0.0,0.0,6.819776136699029,0.0,-25.47923898089684,2.5194527603132846,74.99857065834307,0.0,0.0,88.42327610251432,0.6820330826047487,3.3715482573829685,0.0,0.0]
var1 = [3,2,2,1,2,2,1,1,2,1,3,2,2,1,1,2,2,2,1,2]
tab1 = [-3.823733185508287,0.7857013769778555,9.347700223333298,-0.0,-24.57958470726409,3.549760868834472,-0.0,-0.0,-0.0,-0.0,2.276570692386853,-3.1513940897258808,6.600120188297597,-0.0,-0.0,2.483577554811755,-0.5411982936821893,-0.20714710316669951,-0.0,-0.0]
redcost = repmat(red1,400,1)
varstate = repmat(var1,400,1)
tabrow = repmat(tab1,400,1)
t = time()
for k in 1:1000
ncandidates = 0
for i in 1:n
thisState = varstate[i]
pivotElt = tabrow[i]
if (thisState == AtLower && pivotElt > pivotTol) || (thisState == AtUpper && pivotElt < -pivotTol)
candidates[ncandidates += 1] = i
ratio = 0.
if (pivotElt < 0.)
ratio = (redcost[i] - dualTol)/pivotElt
else
ratio = (redcost[i] + dualTol)/pivotElt
end
if (ratio < thetaMax)
thetaMax = ratio
end
end
end
# pass 2
enter = -1
maxAlpha = 0.
for k in 1:ncandidates
i = candidates[k]
ratio = redcost[i]/tabrow[i]
if (ratio <= thetaMax)
absalpha = abs(tabrow[i])
if (absalpha > maxAlpha)
maxAlpha = absalpha
enter = i
end
end
end
end
return time() - t
end
|