Source

`scripts/bench/integral_threads.purr`

2# Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.

3# integral_threads.purr — the trapezoid integral of integral.purr, parallelized with the

4# `thread` module: N workers each integrate a contiguous chunk into a LOCAL sum, then make ONE

5# lease write into a shared memory.Owner<float> — coordination at the edges, never in the hot

6# loop. Prints, for 1/2/4/8 workers: the integral (must match the single-thread value) and the

7# wall-clock seconds. The single-file baseline lives in integral.py / integral.purr.

8#

9# STATISTICS AND STRIATION. This used to run each worker count exactly once, in the order

10# 1, 2, 4, 8 — and docs/performance.md then described the resulting table as the "median of

11# repeated runs", which it could not have been. Two things are wrong with one shot per

12# configuration. There is no dispersion, so a single descheduled run is indistinguishable

13# from a real effect; and running the four configurations as consecutive blocks means the

14# 8-worker row is measured last, on the hottest silicon, which biases exactly the number the

15# table is about.

16#

17# So: ROUNDS passes, and each pass runs all four configurations back to back. Any thermal

18# drift across the session now moves all four rows together instead of tilting the scaling

19# curve, and each row reports a median with its spread rather than one sample.

20import io

21import math

22import memory

23import statistics

24import stamp

25import sys

26import thread

27import time

29# One pass runs all four configurations; ROUNDS passes are taken. 7 is the smallest count

30# that gives a median a real middle and still leaves a visible [min, max] range.

31fn rounds() { return 7 }

33struct Sample {

34 integral: float

35 secs: float

36}

38fn points() { return 20000000 }

39fn width() { return 50.0 / points() }

41# Integrate f(x) = sin(x)·e^(−x/100) over grid points [lo, hi) into a local sum; one lease

42# write at the end (an Owner lease in a 20M-iteration loop would time the lock, not the math).

43fn chunk(o : memory.Owner<float>, lo : int, hi : int) {

44 let h = width()

45 let s = 0.0

46 for i in range(lo, hi) {

47 let x = h * i

48 s = s + math.sin(x) * math.exp(-0.01 * x)

49 }

50 with o.rwrite().acquire() as w {

51 w.write(w.read() + s)

52 }

53}

55# The trapezoid endpoint correction (the half-weight first/last samples of integral.purr).

56fn edges() -> float {

57 let b = 50.0

58 return 0.5 * (math.sin(0.0) * math.exp(0.0) + math.sin(b) * math.exp(-0.01 * b))

59}

61fn integral_of(total : float) -> float { return (total + edges()) * width() }

63# median, spread and range for one configuration's samples. `statistics` has no quantiles(),

64# and there is no sort builtin to compute an IQR by hand, so the spread reported is the

65# standard deviation plus the observed [min, max] — enough for a reader to see whether a

66# gap between two rows is larger than the noise that produced it.

67fn report(tag : str, integral : float, samples : list<float>) {

68 let lo = samples[0]

69 let hi = samples[0]

70 for v in samples {

71 if v < lo { lo = v }

72 if v > hi { hi = v }

73 }

74 io.print(tag, integral, statistics.median(samples), statistics.stdev(samples), lo, hi)

75}

77fn ms(v : float) -> str { return io.fixed(v * 1000.0, 0) + " ms" }

79# One row of the published table. `speedup` is against the single-worker median; the spread is

80# the standard deviation, because `statistics` has no quantiles() and there is no sort builtin

81# to compute an IQR by hand — reported as `σ` so the reader knows which it is.

82fn row(workers : str, samples : list<float>, base : float, integral : float) -> str {

83 let med = statistics.median(samples)

84 let sp = "—"

85 if base > 0.0 and workers != "1" { sp = "**" + io.fixed(base / med, 2) + "×**" }

86 let out = "| " + workers + " | " + ms(med) + " | ±" + ms(statistics.stdev(samples))

87 out = out + " | " + sp + " | " + io.str(integral) + " |\n"

88 return out

89}

91fn run_1() {

92 let o = memory.own(0.0)

93 let t0 = time.monotonic()

94 chunk(o, 1, points())

95 let t1 = time.monotonic()

96 return Sample(integral_of(o.rread().acquire().read()), t1 - t0)

97}

99fn run_2() {

100 let o = memory.own(0.0)

101 let n = points()

102 let t0 = time.monotonic()

103 with thread.spawn(chunk, o, 1, n / 2) {

104 with thread.spawn(chunk, o, n / 2, n) {

105 }

106 }

107 let t1 = time.monotonic()

108 return Sample(integral_of(o.rread().acquire().read()), t1 - t0)

109}

110

111fn run_4() {

112 let o = memory.own(0.0)

113 let n = points()

114 let q = n / 4

115 let t0 = time.monotonic()

116 with thread.spawn(chunk, o, 1, q) {

117 with thread.spawn(chunk, o, q, 2 * q) {

118 with thread.spawn(chunk, o, 2 * q, 3 * q) {

119 with thread.spawn(chunk, o, 3 * q, n) {

120 }

121 }

122 }

123 }

124 let t1 = time.monotonic()

125 return Sample(integral_of(o.rread().acquire().read()), t1 - t0)

126}

127

128fn run_8() {

129 let o = memory.own(0.0)

130 let n = points()

131 let e = n / 8

132 let t0 = time.monotonic()

133 with thread.spawn(chunk, o, 1, e) {

134 with thread.spawn(chunk, o, e, 2 * e) {

135 with thread.spawn(chunk, o, 2 * e, 3 * e) {

136 with thread.spawn(chunk, o, 3 * e, 4 * e) {

137 with thread.spawn(chunk, o, 4 * e, 5 * e) {

138 with thread.spawn(chunk, o, 5 * e, 6 * e) {

139 with thread.spawn(chunk, o, 6 * e, 7 * e) {

140 with thread.spawn(chunk, o, 7 * e, n) {

141 }

142 }

143 }

144 }

145 }

146 }

147 }

148 }

149 let t1 = time.monotonic()

150 return Sample(integral_of(o.rread().acquire().read()), t1 - t0)

151}

152

153fn main() {

154 let s1: list<float> = []

155 let s2: list<float> = []

156 let s4: list<float> = []

157 let s8: list<float> = []

158 let i1 = 0.0

159 let i2 = 0.0

160 let i4 = 0.0

161 let i8 = 0.0

162 # Striated: one round touches every configuration before any configuration is repeated.

163 for r in range(0, rounds()) {

164 let a = run_1()

165 let b = run_2()

166 let c = run_4()

167 let d = run_8()

168 s1.append(a.secs)

169 s2.append(b.secs)

170 s4.append(c.secs)

171 s8.append(d.secs)

172 i1 = a.integral

173 i2 = b.integral

174 i4 = c.integral

175 i8 = d.integral

176 }

177 io.print("workers integral median_s stdev_s min_s max_s")

178 report("1:", i1, s1)

179 report("2:", i2, s2)

180 report("4:", i4, s4)

181 report("8:", i8, s8)

182

183 # A published table is GENERATED or it is not published. sys.argv[1] is where

184 # scripts/bench_table.purr expects the artifact: docs/bench/thread-scaling.md.

185 if len(sys.argv) > 1 {

186 let base = statistics.median(s1)

188 t = t + "|--------:|-------------------:|-------:|--------:|----------|\n"

189 t = t + row("1", s1, base, i1)

190 t = t + row("2", s2, base, i2)

191 t = t + row("4", s4, base, i4)

192 t = t + row("8", s8, base, i8)

193 let st = stamp.build("thread-scaling",

194 "stdlib/thread/, stdlib/memory/, scripts/bench/integral_threads.purr",

195 "none — cheatah against itself at 1/2/4/8 workers",

196 rounds(),

197 "median wall clock; ± is the sample standard deviation",

198

                             "purrc scripts/bench/integral_threads.purr -o /tmp/it.so --import-root scripts && cheatah /tmp/it.so docs/bench/thread-scaling.md")

199 stamp.write_region(sys.argv[1], st, t)

200 }

201}

202

203main()