# Sum of Squares

Streams and lambdas, especially the limited support offered for primitive types, are a fantastic addition to the Java language. They’re not supposed to be fast, but how do these features compare to a good old `for` loop? For a simple calculation amenable to instruction level parallelism, I compare modern and traditional implementations and observe the differences in instructions generated.

### Sum of Squares

The sum of squares is the building block of a linear regression analysis so is ubiquitous in statistical computing. It is associative and therefore data-parallel. I compare four implementations: a sequential stream wrapping an array, a parallel stream wrapping an array, a generative sequential stream and a traditional `for` loop. The benchmark code is on github.

``````
@Param({"1024", "8192"})
int size;

private double[] data;

@Setup(Level.Iteration)
public void init() {
this.data = createDoubleArray(size);
}

@Benchmark
public double SS_SequentialStream() {
return DoubleStream.of(data)
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0D);
}

@Benchmark
public double SS_ParallelStream() {
return DoubleStream.of(data)
.parallel()
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0);
}

@Benchmark
public double SS_ForLoop() {
double result = 0D;
for (int i = 0; i < data.length; ++i) {
result += data[i] * data[i];
}
return result;
}

@Benchmark
public double SS_GenerativeSequentialStream() {
return IntStream.iterate(0, i -> i < size, i -> i + 1)
.mapToDouble(i -> data[i])
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0);
}
``````

I must admit I prefer the readability of the stream versions, but let’s see if there is a comedown after the syntactic sugar rush.

### Running a Benchmark

I compare the four implementations on an array of one million doubles. I am using `JDK 9.0.1, VM 9.0.1+11` on a fairly powerful laptop with 8 processors:

```\$ cat /proc/cpuinfo
processor       : 0
vendor_id       : GenuineIntel
cpu family      : 6
model           : 94
model name      : Intel(R) Core(TM) i7-6700HQ CPU @ 2.60GHz
stepping        : 3
cpu MHz         : 2592.000
cache size      : 256 KB
physical id     : 0
siblings        : 8
core id         : 0
cpu cores       : 4
apicid          : 0
initial apicid  : 0
fpu             : yes
fpu_exception   : yes
cpuid level     : 22
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi mmx fxsr sse sse2 ss ht tm pbe pni dtes64 monitor ds_cpl vmx est tm2 ssse3 fma cx16 xtpr pdcm sse4_1 sse4_2 x2apic movbe popcnt aes xsave osxsave avx f16c rdrand lahf_lm ida arat epb xsaveopt pln pts dtherm fsgsbase tsc_adjust bmi1 hle avx2 smep bmi2 erms invpcid rtm mpx rdseed adx smap clflushopt
clflush size    : 64
cache_alignment : 64
address sizes   : 39 bits physical, 48 bits virtual
power management:
```

Before running the benchmark we might expect the `for` loop and stream to have similar performance, and the parallel version to be about eight times faster (though remember that the arrays aren’t too big). The generative version is very similar to the `for` loop so a slow down might not be expected.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
SS_ForLoop thrpt 1 10 258351.774491 39797.567968 ops/s 1024
SS_ForLoop thrpt 1 10 29463.408428 4814.826388 ops/s 8192
SS_GenerativeSequentialStream thrpt 1 10 219699.607567 9095.569546 ops/s 1024
SS_GenerativeSequentialStream thrpt 1 10 28351.900454 828.513989 ops/s 8192
SS_ParallelStream thrpt 1 10 22827.821827 2826.577213 ops/s 1024
SS_ParallelStream thrpt 1 10 23230.623610 273.415352 ops/s 8192
SS_SequentialStream thrpt 1 10 225431.985145 9051.538442 ops/s 1024
SS_SequentialStream thrpt 1 10 29123.734157 1333.721437 ops/s 8192

The `for` loop and stream are similar. The parallel version is a long way behind (yes that’s right: more threads less power), but exhibits constant scaling (incidentally, a measurement like this is a good way to guess the minimum unit of work in a parallelised implementation). If the data is large it could become profitable to use it. The generative stream is surprisingly good, almost as good as the version that wraps the array, though there is a fail-safe way to slow it down: add a limit clause to the method chain (try it…).

Profiling with perfasm, it is clear that the `for` loop body is being vectorised, but only the loads and multiplications are done in parallel – the complicated string of SSE instructions is the reduction, which must be done in order.

```<-- unrolled load -->
0.01%    0x00000243d8969170: vmovdqu ymm1,ymmword ptr [r11+r8*8+0f0h]
0.07%    0x00000243d896917a: vmovdqu ymm2,ymmword ptr [r11+r8*8+0d0h]
0.75%    0x00000243d8969184: vmovdqu ymm3,ymmword ptr [r11+r8*8+0b0h]
0.01%    0x00000243d896918e: vmovdqu ymm4,ymmword ptr [r11+r8*8+90h]
0.02%    0x00000243d8969198: vmovdqu ymm5,ymmword ptr [r11+r8*8+70h]
0.03%    0x00000243d896919f: vmovdqu ymm6,ymmword ptr [r11+r8*8+50h]
0.77%    0x00000243d89691a6: vmovdqu ymm10,ymmword ptr [r11+r8*8+30h]
0.02%    0x00000243d89691ad: vmovdqu ymm7,ymmword ptr [r11+r8*8+10h]
<-- multiplication starts -->
0.01%    0x00000243d89691b4: vmulpd  ymm1,ymm1,ymm1
0.02%    0x00000243d89691b8: vmovdqu ymmword ptr [rsp+28h],ymm1
0.76%    0x00000243d89691be: vmulpd  ymm15,ymm7,ymm7
0.00%    0x00000243d89691c2: vmulpd  ymm12,ymm2,ymm2
0.01%    0x00000243d89691c6: vmulpd  ymm7,ymm3,ymm3
0.02%    0x00000243d89691ca: vmulpd  ymm8,ymm4,ymm4
0.72%    0x00000243d89691ce: vmulpd  ymm9,ymm5,ymm5
0.00%    0x00000243d89691d2: vmulpd  ymm11,ymm6,ymm6
0.01%    0x00000243d89691d6: vmulpd  ymm13,ymm10,ymm10
<-- multiplication ends here, scalar reduction starts -->
0.03%    0x00000243d89691db: vaddsd  xmm0,xmm0,xmm15
0.72%    0x00000243d89691e0: vpshufd xmm5,xmm15,0eh
0.01%    0x00000243d89691e6: vaddsd  xmm0,xmm0,xmm5
2.14%    0x00000243d89691ea: vextractf128 xmm6,ymm15,1h
0.03%    0x00000243d89691f0: vaddsd  xmm0,xmm0,xmm6
3.21%    0x00000243d89691f4: vpshufd xmm5,xmm6,0eh
0.02%    0x00000243d89691f9: vaddsd  xmm0,xmm0,xmm5
2.81%    0x00000243d89691fd: vaddsd  xmm0,xmm0,xmm13
2.82%    0x00000243d8969202: vpshufd xmm5,xmm13,0eh
0.03%    0x00000243d8969208: vaddsd  xmm0,xmm0,xmm5
2.87%    0x00000243d896920c: vextractf128 xmm6,ymm13,1h
0.01%    0x00000243d8969212: vaddsd  xmm0,xmm0,xmm6
3.03%    0x00000243d8969216: vpshufd xmm5,xmm6,0eh
0.03%    0x00000243d896921b: vaddsd  xmm0,xmm0,xmm5
2.94%    0x00000243d896921f: vaddsd  xmm0,xmm0,xmm11
2.70%    0x00000243d8969224: vpshufd xmm5,xmm11,0eh
0.03%    0x00000243d896922a: vaddsd  xmm0,xmm0,xmm5
2.98%    0x00000243d896922e: vextractf128 xmm6,ymm11,1h
0.01%    0x00000243d8969234: vaddsd  xmm0,xmm0,xmm6
3.11%    0x00000243d8969238: vpshufd xmm5,xmm6,0eh
0.03%    0x00000243d896923d: vaddsd  xmm0,xmm0,xmm5
2.95%    0x00000243d8969241: vaddsd  xmm0,xmm0,xmm9
2.61%    0x00000243d8969246: vpshufd xmm5,xmm9,0eh
0.02%    0x00000243d896924c: vaddsd  xmm0,xmm0,xmm5
2.89%    0x00000243d8969250: vextractf128 xmm6,ymm9,1h
0.04%    0x00000243d8969256: vaddsd  xmm0,xmm0,xmm6
3.13%    0x00000243d896925a: vpshufd xmm5,xmm6,0eh
0.01%    0x00000243d896925f: vaddsd  xmm0,xmm0,xmm5
2.96%    0x00000243d8969263: vaddsd  xmm0,xmm0,xmm8
2.83%    0x00000243d8969268: vpshufd xmm4,xmm8,0eh
0.01%    0x00000243d896926e: vaddsd  xmm0,xmm0,xmm4
3.00%    0x00000243d8969272: vextractf128 xmm10,ymm8,1h
0.02%    0x00000243d8969278: vaddsd  xmm0,xmm0,xmm10
3.13%    0x00000243d896927d: vpshufd xmm4,xmm10,0eh
0.01%    0x00000243d8969283: vaddsd  xmm0,xmm0,xmm4
3.01%    0x00000243d8969287: vaddsd  xmm0,xmm0,xmm7
2.95%    0x00000243d896928b: vpshufd xmm1,xmm7,0eh
0.02%    0x00000243d8969290: vaddsd  xmm0,xmm0,xmm1
3.06%    0x00000243d8969294: vextractf128 xmm2,ymm7,1h
0.01%    0x00000243d896929a: vaddsd  xmm0,xmm0,xmm2
3.07%    0x00000243d896929e: vpshufd xmm1,xmm2,0eh
0.02%    0x00000243d89692a3: vaddsd  xmm0,xmm0,xmm1
3.07%    0x00000243d89692a7: vaddsd  xmm0,xmm0,xmm12
2.92%    0x00000243d89692ac: vpshufd xmm3,xmm12,0eh
0.02%    0x00000243d89692b2: vaddsd  xmm0,xmm0,xmm3
3.11%    0x00000243d89692b6: vextractf128 xmm1,ymm12,1h
0.01%    0x00000243d89692bc: vaddsd  xmm0,xmm0,xmm1
3.02%    0x00000243d89692c0: vpshufd xmm3,xmm1,0eh
0.02%    0x00000243d89692c5: vaddsd  xmm0,xmm0,xmm3
2.97%    0x00000243d89692c9: vmovdqu ymm1,ymmword ptr [rsp+28h]
0.02%    0x00000243d89692cf: vaddsd  xmm0,xmm0,xmm1
3.05%    0x00000243d89692d3: vpshufd xmm2,xmm1,0eh
0.03%    0x00000243d89692d8: vaddsd  xmm0,xmm0,xmm2
2.97%    0x00000243d89692dc: vextractf128 xmm14,ymm1,1h
0.01%    0x00000243d89692e2: vaddsd  xmm0,xmm0,xmm14
2.99%    0x00000243d89692e7: vpshufd xmm2,xmm14,0eh
0.02%    0x00000243d89692ed: vaddsd  xmm0,xmm0,xmm2
```

The sequential stream code is not as good – it is scalar – but the difference in performance is not as stark as it might be because of the inefficient scalar reduction in the `for` loop: this is JLS floating point semantics twisting C2’s arm behind its back.

```  0.00%    0x0000021a1df54c24: vmovsd  xmm0,qword ptr [rbx+r9*8+48h]
0.00%    0x0000021a1df54c2b: vmovsd  xmm2,qword ptr [rbx+r9*8+18h]
0.02%    0x0000021a1df54c32: vmovsd  xmm3,qword ptr [rbx+r9*8+40h]
2.93%    0x0000021a1df54c39: vmovsd  xmm4,qword ptr [rbx+r9*8+38h]
0.00%    0x0000021a1df54c40: vmovsd  xmm5,qword ptr [rbx+r9*8+30h]
0.01%    0x0000021a1df54c47: vmovsd  xmm6,qword ptr [rbx+r9*8+28h]
0.02%    0x0000021a1df54c4e: vmovsd  xmm7,qword ptr [rbx+r9*8+20h]
2.99%    0x0000021a1df54c55: vmulsd  xmm8,xmm0,xmm0
0.00%    0x0000021a1df54c59: vmulsd  xmm0,xmm7,xmm7
0x0000021a1df54c5d: vmulsd  xmm6,xmm6,xmm6
0.01%    0x0000021a1df54c61: vmulsd  xmm5,xmm5,xmm5
2.91%    0x0000021a1df54c65: vmulsd  xmm4,xmm4,xmm4
0.00%    0x0000021a1df54c69: vmulsd  xmm3,xmm3,xmm3
0.00%    0x0000021a1df54c6d: vmulsd  xmm2,xmm2,xmm2
0.02%    0x0000021a1df54c71: vaddsd  xmm1,xmm2,xmm1
6.10%    0x0000021a1df54c75: vaddsd  xmm0,xmm0,xmm1
5.97%    0x0000021a1df54c79: vaddsd  xmm0,xmm6,xmm0
16.22%    0x0000021a1df54c7d: vaddsd  xmm0,xmm5,xmm0
7.86%    0x0000021a1df54c81: vaddsd  xmm0,xmm4,xmm0
11.16%    0x0000021a1df54c85: vaddsd  xmm1,xmm3,xmm0
11.90%    0x0000021a1df54c89: vaddsd  xmm0,xmm8,xmm1
```

The same code can be seen in `SS_ParallelStream`. `SS_GenerativeSequentialStream` is much more interesting because it hasn’t been unrolled – see the interleaved control statements. It is also not vectorised.

```           0x0000013c1a639c17: vmovsd  xmm0,qword ptr [rbp+r9*8+10h]
0.01%    0x0000013c1a639c1e: vmulsd  xmm2,xmm0,xmm0
0.01%    0x0000013c1a639c22: test    r8d,r8d
0x0000013c1a639c25: jne     13c1a639e09h
0x0000013c1a639c2b: mov     r10d,dword ptr [r12+rax*8+8h]
0x0000013c1a639c30: cmp     r10d,0f8022d85h
0x0000013c1a639c37: jne     13c1a639e3bh
0.01%    0x0000013c1a639c3d: vaddsd  xmm2,xmm1,xmm2
0.01%    0x0000013c1a639c41: vmovsd  qword ptr [rdi+10h],xmm2
0.00%    0x0000013c1a639c46: movsxd  r10,r9d
0x0000013c1a639c49: vmovsd  xmm0,qword ptr [rbp+r10*8+18h]
0.01%    0x0000013c1a639c50: vmulsd  xmm0,xmm0,xmm0
0.01%    0x0000013c1a639c54: mov     r10d,dword ptr [r12+rax*8+8h]
0.00%    0x0000013c1a639c59: cmp     r10d,0f8022d85h
0x0000013c1a639c60: jne     13c1a639e30h
0.02%    0x0000013c1a639c6a: vmovsd  qword ptr [rdi+10h],xmm0
0.02%    0x0000013c1a639c6f: mov     r10d,r9d
0x0000013c1a639c76: cmp     r10d,r11d
0x0000013c1a639c79: jnl     13c1a639d96h
0.01%    0x0000013c1a639c7f: add     r9d,4h
0.02%    0x0000013c1a639c83: vmovsd  xmm1,qword ptr [rbp+r10*8+10h]
0.00%    0x0000013c1a639c8a: movzx   r8d,byte ptr [rdi+0ch]
0.00%    0x0000013c1a639c8f: vmulsd  xmm1,xmm1,xmm1
0.01%    0x0000013c1a639c93: test    r8d,r8d
0x0000013c1a639c96: jne     13c1a639dfbh
0.01%    0x0000013c1a639c9c: vaddsd  xmm1,xmm0,xmm1
0.01%    0x0000013c1a639ca0: vmovsd  qword ptr [rdi+10h],xmm1
0.02%    0x0000013c1a639ca5: movsxd  r8,r10d
0.00%    0x0000013c1a639ca8: vmovsd  xmm0,qword ptr [rbp+r8*8+18h]
0x0000013c1a639caf: vmulsd  xmm0,xmm0,xmm0
0.06%    0x0000013c1a639cb7: vmovsd  qword ptr [rdi+10h],xmm0
```

So it looks like streams don’t vectorise like good old `for` loops, and you won’t gain from `Stream.parallelStream` unless you have humungous arrays (which you might be avoiding for other reasons). This was actually a very nice case for the `Stream` because optimal code can’t be generated for floating point reductions. What happens with sum of squares for `int`s? Generating data in an unsurprising way:

``````
@Benchmark
public int SS_SequentialStream_Int() {
return IntStream.of(intData)
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0);
}

@Benchmark
public int SS_ParallelStream_Int() {
return IntStream.of(intData)
.parallel()
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0);
}

@Benchmark
public int SS_ForLoop_Int() {
int result = 0;
for (int i = 0; i < intData.length; ++i) {
result += intData[i] * intData[i];
}
return result;
}

@Benchmark
public int SS_GenerativeSequentialStream_Int() {
return IntStream.iterate(0, i -> i < size, i -> i + 1)
.map(i -> intData[i])
.map(x -> x * x)
.reduce((x, y) -> x + y)
.orElse(0);
}
``````

The landscape has completely changed, thanks to the exploitation of associative arithmetic and the `VPHADDD` instruction which simplifies the reduction in the for loop.

```<-- load -->
0.00%    0x000001f5cdd8cd30: vmovdqu ymm0,ymmword ptr [rdi+r10*4+0f0h]
1.93%    0x000001f5cdd8cd3a: vmovdqu ymm1,ymmword ptr [rdi+r10*4+0d0h]
0.10%    0x000001f5cdd8cd44: vmovdqu ymm2,ymmword ptr [rdi+r10*4+0b0h]
0.07%    0x000001f5cdd8cd4e: vmovdqu ymm3,ymmword ptr [rdi+r10*4+90h]
0.05%    0x000001f5cdd8cd58: vmovdqu ymm4,ymmword ptr [rdi+r10*4+70h]
1.75%    0x000001f5cdd8cd5f: vmovdqu ymm5,ymmword ptr [rdi+r10*4+50h]
0.08%    0x000001f5cdd8cd66: vmovdqu ymm6,ymmword ptr [rdi+r10*4+30h]
0.07%    0x000001f5cdd8cd6d: vmovdqu ymm7,ymmword ptr [rdi+r10*4+10h]
<-- multiply -->
0.01%    0x000001f5cdd8cd74: vpmulld ymm0,ymm0,ymm0
1.81%    0x000001f5cdd8cd79: vmovdqu ymmword ptr [rsp+28h],ymm0
0.02%    0x000001f5cdd8cd7f: vpmulld ymm15,ymm7,ymm7
1.79%    0x000001f5cdd8cd84: vpmulld ymm11,ymm1,ymm1
0.06%    0x000001f5cdd8cd89: vpmulld ymm8,ymm2,ymm2
1.82%    0x000001f5cdd8cd8e: vpmulld ymm9,ymm3,ymm3
0.06%    0x000001f5cdd8cd93: vpmulld ymm10,ymm4,ymm4
1.79%    0x000001f5cdd8cd98: vpmulld ymm12,ymm5,ymm5
0.08%    0x000001f5cdd8cd9d: vpmulld ymm6,ymm6,ymm6
<-- vectorised reduce -->
1.83%    0x000001f5cdd8cda2: vphaddd ymm4,ymm15,ymm15
0.04%    0x000001f5cdd8cda7: vphaddd ymm4,ymm4,ymm7
1.85%    0x000001f5cdd8cdac: vextracti128 xmm7,ymm4,1h
0.07%    0x000001f5cdd8cdb2: vpaddd  xmm4,xmm4,xmm7
1.78%    0x000001f5cdd8cdb6: vmovd   xmm7,r8d
0.01%    0x000001f5cdd8cdbb: vpaddd  xmm7,xmm7,xmm4
0.11%    0x000001f5cdd8cdbf: vmovd   r11d,xmm7
0.05%    0x000001f5cdd8cdc4: vphaddd ymm4,ymm6,ymm6
1.84%    0x000001f5cdd8cdc9: vphaddd ymm4,ymm4,ymm7
5.43%    0x000001f5cdd8cdce: vextracti128 xmm7,ymm4,1h
0.13%    0x000001f5cdd8cdd4: vpaddd  xmm4,xmm4,xmm7
4.34%    0x000001f5cdd8cdd8: vmovd   xmm7,r11d
0.36%    0x000001f5cdd8cddd: vpaddd  xmm7,xmm7,xmm4
1.40%    0x000001f5cdd8cde1: vmovd   r8d,xmm7
0.01%    0x000001f5cdd8cde6: vphaddd ymm6,ymm12,ymm12
2.89%    0x000001f5cdd8cdeb: vphaddd ymm6,ymm6,ymm4
3.25%    0x000001f5cdd8cdf0: vextracti128 xmm4,ymm6,1h
0.87%    0x000001f5cdd8cdf6: vpaddd  xmm6,xmm6,xmm4
6.36%    0x000001f5cdd8cdfa: vmovd   xmm4,r8d
0.01%    0x000001f5cdd8cdff: vpaddd  xmm4,xmm4,xmm6
1.69%    0x000001f5cdd8ce03: vmovd   r8d,xmm4
0.03%    0x000001f5cdd8ce08: vphaddd ymm4,ymm10,ymm10
1.83%    0x000001f5cdd8ce0d: vphaddd ymm4,ymm4,ymm7
0.10%    0x000001f5cdd8ce12: vextracti128 xmm7,ymm4,1h
3.29%    0x000001f5cdd8ce18: vpaddd  xmm4,xmm4,xmm7
0.72%    0x000001f5cdd8ce1c: vmovd   xmm7,r8d
0.23%    0x000001f5cdd8ce21: vpaddd  xmm7,xmm7,xmm4
4.42%    0x000001f5cdd8ce25: vmovd   r11d,xmm7
0.12%    0x000001f5cdd8ce2a: vphaddd ymm5,ymm9,ymm9
1.69%    0x000001f5cdd8ce2f: vphaddd ymm5,ymm5,ymm1
0.12%    0x000001f5cdd8ce34: vextracti128 xmm1,ymm5,1h
3.28%    0x000001f5cdd8ce3a: vpaddd  xmm5,xmm5,xmm1
0.22%    0x000001f5cdd8ce3e: vmovd   xmm1,r11d
0.14%    0x000001f5cdd8ce43: vpaddd  xmm1,xmm1,xmm5
3.81%    0x000001f5cdd8ce47: vmovd   r11d,xmm1
0.22%    0x000001f5cdd8ce4c: vphaddd ymm0,ymm8,ymm8
1.58%    0x000001f5cdd8ce51: vphaddd ymm0,ymm0,ymm3
0.22%    0x000001f5cdd8ce56: vextracti128 xmm3,ymm0,1h
2.82%    0x000001f5cdd8ce5c: vpaddd  xmm0,xmm0,xmm3
0.36%    0x000001f5cdd8ce60: vmovd   xmm3,r11d
0.20%    0x000001f5cdd8ce65: vpaddd  xmm3,xmm3,xmm0
4.55%    0x000001f5cdd8ce69: vmovd   r8d,xmm3
0.10%    0x000001f5cdd8ce6e: vphaddd ymm2,ymm11,ymm11
1.71%    0x000001f5cdd8ce73: vphaddd ymm2,ymm2,ymm1
0.09%    0x000001f5cdd8ce78: vextracti128 xmm1,ymm2,1h
2.91%    0x000001f5cdd8ce7e: vpaddd  xmm2,xmm2,xmm1
1.57%    0x000001f5cdd8ce82: vmovd   xmm1,r8d
0.05%    0x000001f5cdd8ce87: vpaddd  xmm1,xmm1,xmm2
4.84%    0x000001f5cdd8ce8b: vmovd   r11d,xmm1
0.06%    0x000001f5cdd8ce90: vmovdqu ymm0,ymmword ptr [rsp+28h]
0.03%    0x000001f5cdd8ce96: vphaddd ymm13,ymm0,ymm0
1.83%    0x000001f5cdd8ce9b: vphaddd ymm13,ymm13,ymm14
2.16%    0x000001f5cdd8cea0: vextracti128 xmm14,ymm13,1h
0.14%    0x000001f5cdd8cea6: vpaddd  xmm13,xmm13,xmm14
0.09%    0x000001f5cdd8ceab: vmovd   xmm14,r11d
0.51%    0x000001f5cdd8ceb0: vpaddd  xmm14,xmm14,xmm13
```

If you’re the guy replacing all the `for` loops with streams because it’s 2018, you may be committing performance vandalism! That nice declarative API (as opposed to language feature) is at arms length and it really isn’t well optimised yet.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
SS_ForLoop_Int thrpt 1 10 1021725.018981 74264.883362 ops/s 1024
SS_ForLoop_Int thrpt 1 10 129250.855026 5764.608094 ops/s 8192
SS_GenerativeSequentialStream_Int thrpt 1 10 55069.227826 1111.903102 ops/s 1024
SS_GenerativeSequentialStream_Int thrpt 1 10 6769.176830 684.970867 ops/s 8192
SS_ParallelStream_Int thrpt 1 10 20970.387258 719.846643 ops/s 1024
SS_ParallelStream_Int thrpt 1 10 19621.397202 1514.374286 ops/s 8192
SS_SequentialStream_Int thrpt 1 10 586847.001223 22390.512706 ops/s 1024
SS_SequentialStream_Int thrpt 1 10 87620.959677 3437.083075 ops/s 8192

Parallel streams might not be the best thing to reach for.

# The Much Aligned Garbage Collector

A power of two is often a good choice for the size of an array. Sometimes you might see this being exploited to replace an integer division with a bitwise intersection. You can see why with a toy benchmark of a bloom filter, which deliberately folds in a representative cost of a hash function and array access to highlight the significance of the differential cost of the division mechanism to a method that does real work:

``````@State(Scope.Thread)
@OutputTimeUnit(TimeUnit.MICROSECONDS)
public class BloomFilter {

private long[] bitset;

@Param({"1000", "1024"})
int size;

@Setup(Level.Trial)
public void init() {
bitset = DataUtil.createLongArray(size);
}

@Benchmark
public boolean containsAnd() {
int hash = hash();
int pos = hash & (size - 1);
return (bitset[pos >>> 6] & (1L << pos)) != 0;
}

@Benchmark
public boolean containsAbsMod() {
int hash = hash();
int pos = Math.abs(hash % size);
return (bitset[pos >>> 6] & (1L << pos)) != 0;
}

private int hash() {
return ThreadLocalRandom.current().nextInt(); // a stand in for a hash function;
}
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
containsAbsMod thrpt 1 10 104.063744 4.068283 ops/us 1000
containsAbsMod thrpt 1 10 103.849577 4.991040 ops/us 1024
containsAnd thrpt 1 10 161.917397 3.807912 ops/us 1024

Disregarding the case which produces an incorrect result, you can do two thirds as many lookups again in the same period of time if you just use a 1024 element bloom filter. Note that the compiler clearly won’t magically transform cases like `AbsMod 1024`; you need to do this yourself. You can readily see this property exploited in any open source bit set, hash set, or bloom filter you care to look at. This is boring, at least, we often get this right by accident. What is quite interesting is a multiplicative decrease in throughput of DAXPY as a result of this same choice of lengths:

``````@OutputTimeUnit(TimeUnit.MICROSECONDS)
public class DAXPYAlignment {

@Param({"250", "256", "1000", "1024"})
int size;

double s;
double[] a;
double[] b;

@Setup(Level.Trial)
public void init() {
a = createDoubleArray(size);
b = createDoubleArray(size);
}

@Benchmark
public void daxpy(Blackhole bh) {
for (int i = 0; i < a.length; ++i) {
a[i] += s * b[i];
}
bh.consume(a);
}
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
daxpy thrpt 1 10 23.499857 0.891309 ops/us 250
daxpy thrpt 1 10 22.425412 0.989512 ops/us 256
daxpy thrpt 1 10 2.420674 0.098991 ops/us 1000
daxpy thrpt 1 10 6.263005 0.175048 ops/us 1024

1000 and 1024 are somehow very different, yet 250 and 256 are almost equivalent. The placement of the second array, which, being allocated on the same thread, will be next to the first array in the TLAB (thread-local allocation buffer) happens to be very unlucky. Let’s allocate an array in between the two we want to loop over, to vary the offsets between the two arrays:

``````  @Param({"0", "6", "12", "18", "24"})
int offset;

double s;
double[] a;
double[] b;

@Setup(Level.Trial)
public void init() {
a = createDoubleArray(size);
padding = new double[offset];
b = createDoubleArray(size);
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: offset Param: size
daxpy thrpt 1 10 2.224875 0.247778 ops/us 0 1000
daxpy thrpt 1 10 6.159791 0.441525 ops/us 0 1024
daxpy thrpt 1 10 2.350425 0.136992 ops/us 6 1000
daxpy thrpt 1 10 6.047009 0.360723 ops/us 6 1024
daxpy thrpt 1 10 3.332370 0.253739 ops/us 12 1000
daxpy thrpt 1 10 6.506141 0.155733 ops/us 12 1024
daxpy thrpt 1 10 6.621031 0.345151 ops/us 18 1000
daxpy thrpt 1 10 6.827635 0.970527 ops/us 18 1024
daxpy thrpt 1 10 7.456584 0.214229 ops/us 24 1000
daxpy thrpt 1 10 7.451441 0.104871 ops/us 24 1024

The pattern is curious (pay attention to the offset parameter) – the ratio of the throughputs for each size ranging from 3x throughput degradation through to parity:

The loop in question is vectorised, which can be disabled by setting `-XX:-UseSuperWord`. Doing so is revealing, because the trend is still present but it is dampened to the extent it could be waved away as noise:

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: offset Param: size
daxpy thrpt 1 10 1.416452 0.079905 ops/us 0 1000
daxpy thrpt 1 10 1.806841 0.200231 ops/us 0 1024
daxpy thrpt 1 10 1.408526 0.085147 ops/us 6 1000
daxpy thrpt 1 10 1.921026 0.049655 ops/us 6 1024
daxpy thrpt 1 10 1.459186 0.076427 ops/us 12 1000
daxpy thrpt 1 10 1.809220 0.199885 ops/us 12 1024
daxpy thrpt 1 10 1.824435 0.169680 ops/us 18 1000
daxpy thrpt 1 10 1.842230 0.204414 ops/us 18 1024
daxpy thrpt 1 10 1.934717 0.229822 ops/us 24 1000
daxpy thrpt 1 10 1.964316 0.039893 ops/us 24 1024

The point is, you may not have cared about alignment much before because it’s unlikely you would have noticed unless you were really looking for it. Decent autovectorisation seems to raise the stakes enormously.

### Analysis with Perfasm

It’s impossible to know for sure what the cause of this behaviour is without profiling. Since I observed this effect on my Windows development laptop, I use xperf via `WinPerfAsmProfiler`, which is part of JMH.

I did some instruction profiling. The same code is going to get generated in each case, with a preloop, main loop and post loop, but by looking at the sampled instruction frequency we can see what’s taking the most time in the vectorised main loop. From now on, superword parallelism is never disabled. The full output of this run can be seen at github. Here is the main loop for size=1024, offset=0, which is unrolled, spending most time loading and storing data (`vmovdqu`) but spending a decent amount of time in the multiplication:

```  0.18%    0x0000020dddc5af90: vmovdqu ymm0,ymmword ptr [r10+r8*8+10h]
9.27%    0x0000020dddc5af97: vmulpd  ymm0,ymm0,ymm2
0.22%    0x0000020dddc5af9b: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+10h]
7.48%    0x0000020dddc5afa2: vmovdqu ymmword ptr [r11+r8*8+10h],ymm0
10.16%    0x0000020dddc5afa9: vmovdqu ymm0,ymmword ptr [r10+r8*8+30h]
0.09%    0x0000020dddc5afb0: vmulpd  ymm0,ymm0,ymm2
3.62%    0x0000020dddc5afb4: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+30h]
10.60%    0x0000020dddc5afbb: vmovdqu ymmword ptr [r11+r8*8+30h],ymm0
0.26%    0x0000020dddc5afc2: vmovdqu ymm0,ymmword ptr [r10+r8*8+50h]
3.76%    0x0000020dddc5afc9: vmulpd  ymm0,ymm0,ymm2
0.20%    0x0000020dddc5afcd: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+50h]
13.23%    0x0000020dddc5afd4: vmovdqu ymmword ptr [r11+r8*8+50h],ymm0
9.46%    0x0000020dddc5afdb: vmovdqu ymm0,ymmword ptr [r10+r8*8+70h]
0.11%    0x0000020dddc5afe2: vmulpd  ymm0,ymm0,ymm2
4.63%    0x0000020dddc5afe6: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+70h]
9.78%    0x0000020dddc5afed: vmovdqu ymmword ptr [r11+r8*8+70h],ymm0
```

In the worst performer (size=1000, offset=0) a lot more time is spent on the stores, a much smaller fraction of observed instructions are involved with multiplication or addition. This indicates either a measurement bias (perhaps there’s some mechanism that makes a store/load easier to observe) or an increase in load/store cost.

```  0.24%    0x000002d1a946f510: vmovdqu ymm0,ymmword ptr [r10+r8*8+10h]
3.61%    0x000002d1a946f517: vmulpd  ymm0,ymm0,ymm2
4.63%    0x000002d1a946f51b: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+10h]
9.73%    0x000002d1a946f522: vmovdqu ymmword ptr [r11+r8*8+10h],ymm0
4.34%    0x000002d1a946f529: vmovdqu ymm0,ymmword ptr [r10+r8*8+30h]
2.13%    0x000002d1a946f530: vmulpd  ymm0,ymm0,ymm2
7.77%    0x000002d1a946f534: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+30h]
13.46%    0x000002d1a946f53b: vmovdqu ymmword ptr [r11+r8*8+30h],ymm0
3.37%    0x000002d1a946f542: vmovdqu ymm0,ymmword ptr [r10+r8*8+50h]
0.47%    0x000002d1a946f549: vmulpd  ymm0,ymm0,ymm2
1.47%    0x000002d1a946f54d: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+50h]
13.00%    0x000002d1a946f554: vmovdqu ymmword ptr [r11+r8*8+50h],ymm0
4.24%    0x000002d1a946f55b: vmovdqu ymm0,ymmword ptr [r10+r8*8+70h]
2.40%    0x000002d1a946f562: vmulpd  ymm0,ymm0,ymm2
8.92%    0x000002d1a946f566: vaddpd  ymm0,ymm0,ymmword ptr [r11+r8*8+70h]
14.10%    0x000002d1a946f56d: vmovdqu ymmword ptr [r11+r8*8+70h],ymm0
```

This trend can be seen to generally improve as 1024 is approached from below, and do bear in mind that this is a noisy measure. Interpret the numbers below as probabilities: were you to stop the execution of daxpy at random, at offset zero, you would have a 94% chance of finding yourself within the main vectorised loop. You would have a 50% chance of observing a store, and only 31% chance of observing a multiply or add. As we get further from 1024, the stores dominate the main loop, and the main loop comes to dominate the method. Again, this is approximate. When the arrays aren’t well aligned, we spend less time loading, less time multiplying and adding, and much more time storing.

classification offset = 0 offset = 6 offset = 12 offset = 18 offset = 24
add 22.79 21.46 15.41 7.77 8.03
load 12.19 11.95 15.55 21.9 21.19
multiply 8.61 7.7 9.54 13.15 8.33
store 50.29 51.3 49.16 42.34 44.56
main loop 93.88 92.41 89.66 85.16 82.11

The effect observed here is also a contributing factor to fluctuations in throughput observed in JDK-8150730.

### Garbage Collection

Is it necessary to make sure all arrays are of a size equal to a power of two and aligned with pages? In this microbenchmark, it’s easy to arrange that, for typical developers this probably isn’t feasible (which isn’t to say there aren’t people out there who do this). Fortunately, this isn’t necessary for most use cases. True to the title, this post has something to do with garbage collection. The arrays were allocated in order, and no garbage would be produced during the benchmarks, so the second array will be split across pages. Let’s put some code into the initialisation of the benchmark bound to trigger garbage collection:

``````  String acc = "";

@Setup(Level.Trial)
public void init() {
a = createDoubleArray(size);
b = createDoubleArray(size);
// don't do this in production
for (int i = 0; i < 10000; ++i) {
acc += UUID.randomUUID().toString();
}
}
``````

A miracle occurs: the code speeds up!

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
daxpy thrpt 1 10 6.854161 0.261247 ops/us 1000
daxpy thrpt 1 10 6.328602 0.163391 ops/us 1024

Why? G1 has rearranged the heap and that second array is no longer right next to the first array, and the bad luck of the initial placement has been undone. This makes the cost of garbage collection difficult to quantify, because if it takes resources with one hand it gives them back with another.

The benchmark code is available at github.

# Multiplying Matrices, Fast and Slow

I recently read a very interesting blog post about exposing Intel SIMD intrinsics via a fork of the Scala compiler (scala-virtualized), which reports multiplicative improvements in throughput over HotSpot JIT compiled code. The academic paper (SIMD Intrinsics on Managed Language Runtimes), which has been accepted at CGO 2018, proposes a powerful alternative to the traditional JVM approach of pairing dumb programmers with a (hopefully) smart JIT compiler. Lightweight Modular Staging (LMS) allows the generation of an executable binary from a high level representation: handcrafted representations of vectorised algorithms, written in a dialect of Scala, can be compiled natively and later invoked with a single JNI call. This approach bypasses C2 without incurring excessive JNI costs. The freely available benchmarks can be easily run to reproduce the results in the paper, which is an achievement in itself, but some of the Java implementations used as baselines look less efficient than they could be. This post is about improving the efficiency of the Java matrix multiplication the LMS generated code is benchmarked against. Despite finding edge cases where autovectorisation fails, I find it is possible to get performance comparable to LMS with plain Java (and a JDK upgrade).

Two implementations of Java matrix multiplication are provided in the NGen benchmarks: `JMMM.baseline` – a naive but cache unfriendly matrix multiplication – and `JMMM.blocked` which is supplied as an improvement. `JMMM.blocked` is something of a local maximum because it does manual loop unrolling: this actually removes the trigger for autovectorisation analysis. I provide a simple and cache-efficient Java implementation (with the same asymptotic complexity, the improvement is just technical) and benchmark these implementations using JDK8 and the soon to be released JDK10 separately.

``````public void fast(float[] a, float[] b, float[] c, int n) {
int in = 0;
for (int i = 0; i < n; ++i) {
int kn = 0;
for (int k = 0; k < n; ++k) {
float aik = a[in + k];
for (int j = 0; j < n; ++j) {
c[in + j] += aik * b[kn + j];
}
kn += n;
}
in += n;
}
}
``````

With JDK 1.8.0_131, the “fast” implementation is only 2x faster than the blocked algorithm; this is nowhere near fast enough to match LMS. In fact, LMS does a lot better than 5x blocked (6x-8x) on my Skylake laptop at 2.6GHz, and performs between 2x and 4x better than the improved implementation. Flops / Cycle is calculated as `size ^ 3 * 2 / CPU frequency Hz`.

```====================================================
Benchmarking MMM.jMMM.fast (JVM implementation)
----------------------------------------------------
Size (N) | Flops / Cycle
----------------------------------------------------
8 | 0.4994459272
32 | 1.0666533335
64 | 0.9429120397
128 | 0.9692385519
192 | 0.9796619688
256 | 1.0141446247
320 | 0.9894415771
384 | 1.0046245750
448 | 1.0221353392
512 | 0.9943527764
576 | 0.9952093603
640 | 0.9854689714
704 | 0.9947153752
768 | 1.0197765248
832 | 1.0479691069
896 | 1.0060121097
960 | 0.9937347412
1024 | 0.9056494897
====================================================

====================================================
Benchmarking MMM.nMMM.blocked (LMS generated)
----------------------------------------------------
Size (N) | Flops / Cycle
----------------------------------------------------
8 | 0.2500390686
32 | 3.9999921875
64 | 4.1626523901
128 | 4.4618695374
192 | 3.9598982956
256 | 4.3737341517
320 | 4.2412225389
384 | 3.9640163416
448 | 4.0957167537
512 | 3.3801071278
576 | 4.1869326167
640 | 3.8225244883
704 | 3.8648224140
768 | 3.5240611589
832 | 3.7941562681
896 | 3.1735179981
960 | 2.5856903789
1024 | 1.7817152313
====================================================

====================================================
Benchmarking MMM.jMMM.blocked (JVM implementation)
----------------------------------------------------
Size (N) | Flops / Cycle
----------------------------------------------------
8 | 0.3333854248
32 | 0.6336670915
64 | 0.5733484649
128 | 0.5987433798
192 | 0.5819900921
256 | 0.5473562109
320 | 0.5623263520
384 | 0.5583823292
448 | 0.5657882256
512 | 0.5430879470
576 | 0.5269635678
640 | 0.5595204791
704 | 0.5297557807
768 | 0.5493631388
832 | 0.5471832673
896 | 0.4769554752
960 | 0.4985080443
1024 | 0.4014589400
====================================================
```

JDK10 is about to be released so it’s worth looking at the effect of recent improvements to C2, including better use of AVX2 and support for vectorised FMA. Since LMS depends on scala-virtualized, which currently only supports Scala 2.11, the LMS implementation cannot be run with a more recent JDK so its performance running in JDK10 could only be extrapolated. Since its raison d’Ãªtre is to bypass C2, it could be reasonably assumed it is insulated from JVM performance improvements (or regressions). Measurements of floating point operations per cycle provide a sensible comparison, in any case.

Moving away from ScalaMeter, I created a JMH benchmark to see how matrix multiplication behaves in JDK10.

``````@OutputTimeUnit(TimeUnit.SECONDS)
@State(Scope.Benchmark)
public class MMM {

@Param({"8", "32", "64", "128", "192", "256", "320", "384", "448", "512" , "576", "640", "704", "768", "832", "896", "960", "1024"})
int size;

private float[] a;
private float[] b;
private float[] c;

@Setup(Level.Trial)
public void init() {
a = DataUtil.createFloatArray(size * size);
b = DataUtil.createFloatArray(size * size);
c = new float[size * size];
}

@Benchmark
public void fast(Blackhole bh) {
fast(a, b, c, size);
bh.consume(c);
}

@Benchmark
public void baseline(Blackhole bh) {
baseline(a, b, c, size);
bh.consume(c);
}

@Benchmark
public void blocked(Blackhole bh) {
blocked(a, b, c, size);
bh.consume(c);
}

//
// Baseline implementation of a Matrix-Matrix-Multiplication
//
public void baseline (float[] a, float[] b, float[] c, int n){
for (int i = 0; i < n; i += 1) {
for (int j = 0; j < n; j += 1) {
float sum = 0.0f;
for (int k = 0; k < n; k += 1) {
sum += a[i * n + k] * b[k * n + j];
}
c[i * n + j] = sum;
}
}
}

//
// Blocked version of MMM, reference implementation available at:
// http://csapp.cs.cmu.edu/2e/waside/waside-blocking.pdf
//
public void blocked(float[] a, float[] b, float[] c, int n) {
int BLOCK_SIZE = 8;
for (int kk = 0; kk < n; kk += BLOCK_SIZE) {
for (int jj = 0; jj < n; jj += BLOCK_SIZE) {
for (int i = 0; i < n; i++) {
for (int j = jj; j < jj + BLOCK_SIZE; ++j) {
float sum = c[i * n + j];
for (int k = kk; k < kk + BLOCK_SIZE; ++k) {
sum += a[i * n + k] * b[k * n + j];
}
c[i * n + j] = sum;
}
}
}
}
}

public void fast(float[] a, float[] b, float[] c, int n) {
int in = 0;
for (int i = 0; i < n; ++i) {
int kn = 0;
for (int k = 0; k < n; ++k) {
float aik = a[in + k];
for (int j = 0; j < n; ++j) {
c[in + j] = Math.fma(aik,  b[kn + j], c[in + j]);
}
kn += n;
}
in += n;
}
}
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size Ratio to blocked Flops/Cycle
baseline thrpt 1 10 1228544.82 38793.17392 ops/s 8 1.061598336 0.483857652
baseline thrpt 1 10 22973.03402 1012.043446 ops/s 32 1.302266947 0.57906183
baseline thrpt 1 10 2943.088879 221.57475 ops/s 64 1.301414733 0.593471609
baseline thrpt 1 10 358.010135 9.342801 ops/s 128 1.292889618 0.577539747
baseline thrpt 1 10 105.758366 4.275503 ops/s 192 1.246415143 0.575804515
baseline thrpt 1 10 41.465557 1.112753 ops/s 256 1.430003946 0.535135851
baseline thrpt 1 10 20.479081 0.462547 ops/s 320 1.154267894 0.516198866
baseline thrpt 1 10 11.686685 0.263476 ops/s 384 1.186535349 0.509027985
baseline thrpt 1 10 7.344184 0.269656 ops/s 448 1.166421127 0.507965526
baseline thrpt 1 10 3.545153 0.108086 ops/s 512 0.81796657 0.366017216
baseline thrpt 1 10 3.789384 0.130934 ops/s 576 1.327168294 0.557048123
baseline thrpt 1 10 1.981957 0.040136 ops/s 640 1.020965271 0.399660104
baseline thrpt 1 10 1.76672 0.036386 ops/s 704 1.168272442 0.474179037
baseline thrpt 1 10 1.01026 0.049853 ops/s 768 0.845514112 0.352024966
baseline thrpt 1 10 1.115814 0.03803 ops/s 832 1.148752171 0.494331667
baseline thrpt 1 10 0.703561 0.110626 ops/s 896 0.938435436 0.389298235
baseline thrpt 1 10 0.629896 0.052448 ops/s 960 1.081741651 0.428685898
baseline thrpt 1 10 0.407772 0.019079 ops/s 1024 1.025356561 0.336801424
blocked thrpt 1 10 1157259.558 49097.48711 ops/s 8 1 0.455782226
blocked thrpt 1 10 17640.8025 1226.401298 ops/s 32 1 0.444656782
blocked thrpt 1 10 2261.453481 98.937035 ops/s 64 1 0.456020355
blocked thrpt 1 10 276.906961 22.851857 ops/s 128 1 0.446704605
blocked thrpt 1 10 84.850033 4.441454 ops/s 192 1 0.461968485
blocked thrpt 1 10 28.996813 7.585551 ops/s 256 1 0.374219842
blocked thrpt 1 10 17.742052 0.627629 ops/s 320 1 0.447208892
blocked thrpt 1 10 9.84942 0.367603 ops/s 384 1 0.429003641
blocked thrpt 1 10 6.29634 0.402846 ops/s 448 1 0.435490676
blocked thrpt 1 10 4.334105 0.384849 ops/s 512 1 0.447472097
blocked thrpt 1 10 2.85524 0.199102 ops/s 576 1 0.419726816
blocked thrpt 1 10 1.941258 0.10915 ops/s 640 1 0.391453182
blocked thrpt 1 10 1.51225 0.076621 ops/s 704 1 0.40588053
blocked thrpt 1 10 1.194847 0.063147 ops/s 768 1 0.416344283
blocked thrpt 1 10 0.971327 0.040421 ops/s 832 1 0.430320551
blocked thrpt 1 10 0.749717 0.042997 ops/s 896 1 0.414837526
blocked thrpt 1 10 0.582298 0.016725 ops/s 960 1 0.39629231
blocked thrpt 1 10 0.397688 0.043639 ops/s 1024 1 0.328472491
fast thrpt 1 10 1869676.345 76416.50848 ops/s 8 1.615606743 0.736364837
fast thrpt 1 10 48485.47216 1301.926828 ops/s 32 2.748484496 1.222132271
fast thrpt 1 10 6431.341657 153.905413 ops/s 64 2.843897392 1.296875098
fast thrpt 1 10 840.601821 45.998723 ops/s 128 3.035683242 1.356053685
fast thrpt 1 10 260.386996 13.022418 ops/s 192 3.068790745 1.417684611
fast thrpt 1 10 107.895708 6.584674 ops/s 256 3.720950575 1.392453537
fast thrpt 1 10 56.245336 2.729061 ops/s 320 3.170170846 1.417728592
fast thrpt 1 10 32.917996 2.196624 ops/s 384 3.342125323 1.433783932
fast thrpt 1 10 20.960189 2.077684 ops/s 448 3.328948087 1.449725854
fast thrpt 1 10 14.005186 0.7839 ops/s 512 3.231390564 1.445957112
fast thrpt 1 10 8.827584 0.883654 ops/s 576 3.091713481 1.297675056
fast thrpt 1 10 7.455607 0.442882 ops/s 640 3.840605937 1.503417416
fast thrpt 1 10 5.322894 0.464362 ops/s 704 3.519850554 1.428638807
fast thrpt 1 10 4.308522 0.153846 ops/s 768 3.605919419 1.501303934
fast thrpt 1 10 3.375274 0.106715 ops/s 832 3.474910097 1.495325228
fast thrpt 1 10 2.320152 0.367881 ops/s 896 3.094703735 1.28379924
fast thrpt 1 10 2.057478 0.150198 ops/s 960 3.533376381 1.400249889
fast thrpt 1 10 1.66255 0.181116 ops/s 1024 4.180538513 1.3731919

Interestingly, the blocked algorithm is now the worst native JVM implementation. The code generated by C2 got a lot faster, but peaks at 1.5 flops/cycle, which still doesn’t compete with LMS. Why? Taking a look at the assembly, it’s clear that the autovectoriser choked on the array offsets and produced scalar SSE2 code, just like the implementations in the paper. I wasn’t expecting this.

```vmovss  xmm5,dword ptr [rdi+rcx*4+10h]
vmovss  dword ptr [rdi+rcx*4+10h],xmm5
```

Is this the end of the story? No, with some hacks and the cost of array allocation and a copy or two, autovectorisation can be tricked into working again to generate faster code:

``````
public void fast(float[] a, float[] b, float[] c, int n) {
float[] bBuffer = new float[n];
float[] cBuffer = new float[n];
int in = 0;
for (int i = 0; i < n; ++i) {
int kn = 0;
for (int k = 0; k < n; ++k) {
float aik = a[in + k];
System.arraycopy(b, kn, bBuffer, 0, n);
saxpy(n, aik, bBuffer, cBuffer);
kn += n;
}
System.arraycopy(cBuffer, 0, c, in, n);
Arrays.fill(cBuffer, 0f);
in += n;
}
}

private void saxpy(int n, float aik, float[] b, float[] c) {
for (int i = 0; i < n; ++i) {
c[i] += aik * b[i];
}
}
``````

Adding this hack into the NGen benchmark (back in JDK 1.8.0_131) I get closer to the LMS generated code, and beat it beyond L3 cache residency (6MB). LMS is still faster when both matrices fit in L3 concurrently, but by percentage points rather than a multiple. The cost of the hacky array buffers gives the game up for small matrices.

```====================================================
Benchmarking MMM.jMMM.fast (JVM implementation)
----------------------------------------------------
Size (N) | Flops / Cycle
----------------------------------------------------
8 | 0.2500390686
32 | 0.7710872405
64 | 1.1302489072
128 | 2.5113453810
192 | 2.9525859816
256 | 3.1180920385
320 | 3.1081563593
384 | 3.1458423577
448 | 3.0493148252
512 | 3.0551158263
576 | 3.1430376938
640 | 3.2169923048
704 | 3.1026513283
768 | 2.4190053777
832 | 3.3358586705
896 | 3.0755689237
960 | 2.9996690697
1024 | 2.2935654309
====================================================

====================================================
Benchmarking MMM.nMMM.blocked (LMS generated)
----------------------------------------------------
Size (N) | Flops / Cycle
----------------------------------------------------
8 | 1.0001562744
32 | 5.3330416826
64 | 5.8180867784
128 | 5.1717318641
192 | 5.1639907462
256 | 4.3418618628
320 | 5.2536572701
384 | 4.0801359215
448 | 4.1337007093
512 | 3.2678160754
576 | 3.7973028890
640 | 3.3557513664
704 | 4.0103133240
768 | 3.4188362575
832 | 3.2189488327
896 | 3.2316685219
960 | 2.9985655539
1024 | 1.7750946796
====================================================
```

With the benchmark below I calculate flops/cycle with improved JDK10 autovectorisation.

``````
@Benchmark
public void fastBuffered(Blackhole bh) {
fastBuffered(a, b, c, size);
bh.consume(c);
}

public void fastBuffered(float[] a, float[] b, float[] c, int n) {
float[] bBuffer = new float[n];
float[] cBuffer = new float[n];
int in = 0;
for (int i = 0; i < n; ++i) {
int kn = 0;
for (int k = 0; k < n; ++k) {
float aik = a[in + k];
System.arraycopy(b, kn, bBuffer, 0, n);
saxpy(n, aik, bBuffer, cBuffer);
kn += n;
}
System.arraycopy(cBuffer, 0, c, in, n);
Arrays.fill(cBuffer, 0f);
in += n;
}
}

private void saxpy(int n, float aik, float[] b, float[] c) {
for (int i = 0; i < n; ++i) {
c[i] = Math.fma(aik, b[i], c[i]);
}
}
``````

Just as in the modified NGen benchmark, this starts paying off once the matrices have 64 rows and columns. Finally, and it took an upgrade and a hack, I breached 4 Flops per cycle:

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size Flops / Cycle
fastBuffered thrpt 1 10 1047184.034 63532.95095 ops/s 8 0.412429404
fastBuffered thrpt 1 10 58373.56367 3239.615866 ops/s 32 1.471373026
fastBuffered thrpt 1 10 12099.41654 497.33988 ops/s 64 2.439838038
fastBuffered thrpt 1 10 2136.50264 105.038006 ops/s 128 3.446592911
fastBuffered thrpt 1 10 673.470622 102.577237 ops/s 192 3.666730488
fastBuffered thrpt 1 10 305.541519 25.959163 ops/s 256 3.943181586
fastBuffered thrpt 1 10 158.437372 6.708384 ops/s 320 3.993596774
fastBuffered thrpt 1 10 88.283718 7.58883 ops/s 384 3.845306266
fastBuffered thrpt 1 10 58.574507 4.248521 ops/s 448 4.051345968
fastBuffered thrpt 1 10 37.183635 4.360319 ops/s 512 3.839002314
fastBuffered thrpt 1 10 29.949884 0.63346 ops/s 576 4.40270151
fastBuffered thrpt 1 10 20.715833 4.175897 ops/s 640 4.177331789
fastBuffered thrpt 1 10 10.824837 0.902983 ops/s 704 2.905333492
fastBuffered thrpt 1 10 8.285254 1.438701 ops/s 768 2.886995686
fastBuffered thrpt 1 10 6.17029 0.746537 ops/s 832 2.733582608
fastBuffered thrpt 1 10 4.828872 1.316901 ops/s 896 2.671937962
fastBuffered thrpt 1 10 3.6343 1.293923 ops/s 960 2.473381573
fastBuffered thrpt 1 10 2.458296 0.171224 ops/s 1024 2.030442485

The code generated for the core of the loop looks better now:

```vmovdqu ymm1,ymmword ptr [r13+r11*4+10h]
vfmadd231ps ymm1,ymm3,ymmword ptr [r14+r11*4+10h]
vmovdqu ymmword ptr [r13+r11*4+10h],ymm1
```

These benchmark results can be compared on a line chart.

Given this improvement, it would be exciting to see how LMS can profit from JDK9 or JDK10 – does LMS provide the impetus to resume maintenance of scala-virtualized? L3 cache, which the LMS generated code seems to depend on for throughput, is typically shared between cores: a single thread rarely enjoys exclusive access. I would like to see benchmarks for the LMS generated code in the presence of concurrency.

# Autovectorised FMA in JDK10

Fused-multiply-add (FMA) allows floating point expressions of the form `a * x + b` to be evaluated in a single instruction, which is useful for numerical linear algebra. Despite the obvious appeal of FMA, JVM implementors are rather constrained when it comes to floating point arithmetic because Java programs are expected to be reproducible across versions and target architectures. FMA does not produce precisely the same result as the equivalent multiplication and addition instructions (this is caused by the compounding effect of rounding) so its use is a change in semantics rather than an optimisation; the user must opt in. To the best of my knowledge, support for FMA was first proposed in 2000, along with reorderable floating point operations, which would have been activated by a `fastfp` keyword, but this proposal was withdrawn. In Java 9, the intrinsic `Math.fma` was introduced to provide access to FMA for the first time.

### DAXPY Benchmark

A good use case to evaluate `Math.fma` is DAXPY from the Basic Linear Algebra Subroutine library. The code below will compile with JDK9+

``````@OutputTimeUnit(TimeUnit.MILLISECONDS)
public class DAXPY {

double s;

@Setup(Level.Invocation)
public void init() {
}

@Benchmark
public void daxpyFMA(DoubleData state, Blackhole bh) {
double[] a = state.data1;
double[] b = state.data2;
for (int i = 0; i < a.length; ++i) {
a[i] = Math.fma(s, b[i], a[i]);
}
bh.consume(a);
}

@Benchmark
public void daxpy(DoubleData state, Blackhole bh) {
double[] a = state.data1;
double[] b = state.data2;
for (int i = 0; i < a.length; ++i) {
a[i] += s * b[i];
}
bh.consume(a);
}
}
``````

Running this benchmark with Java 9, you may wonder why you bothered because the code is actually slower.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
daxpy thrpt 1 10 25.011242 2.259007 ops/ms 100000
daxpy thrpt 1 10 0.706180 0.046146 ops/ms 1000000
daxpyFMA thrpt 1 10 15.334652 0.271946 ops/ms 100000
daxpyFMA thrpt 1 10 0.623838 0.018041 ops/ms 1000000

This is because using `Math.fma` disables autovectorisation. Taking a look at `PrintAssembly` you can see that the naive `daxpy` routine exploits AVX2, whereas `daxpyFMA` reverts to scalar usage of SSE2.

```// daxpy routine, code taken from main vectorised loop
vmovdqu ymm1,ymmword ptr [r10+rdx*8+10h]
vmulpd  ymm1,ymm1,ymm2
vaddpd  ymm1,ymm1,ymmword ptr [r8+rdx*8+10h]
vmovdqu ymmword ptr [r8+rdx*8+10h],ymm1

// daxpyFMA routine
vmovsd  xmm2,qword ptr [rcx+r13*8+10h]
vmovsd  qword ptr [rcx+r13*8+10h],xmm2
```

Not to worry, this seems to have been fixed in JDK 10. Since Java 10’s release is around the corner, there are early access builds available for all platforms. Rerunning this benchmark, FMA no longer incurs costs, and it doesn’t bring the performance boost some people might expect. The benefit is that there is less floating point error because the total number of roundings is halved.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
daxpy thrpt 1 10 2582.363228 116.637400 ops/ms 1000
daxpy thrpt 1 10 405.904377 32.364782 ops/ms 10000
daxpy thrpt 1 10 25.210111 1.671794 ops/ms 100000
daxpy thrpt 1 10 0.608660 0.112512 ops/ms 1000000
daxpyFMA thrpt 1 10 2650.264580 211.342407 ops/ms 1000
daxpyFMA thrpt 1 10 389.274693 43.567450 ops/ms 10000
daxpyFMA thrpt 1 10 24.941172 2.393358 ops/ms 100000
daxpyFMA thrpt 1 10 0.671310 0.158470 ops/ms 1000000

```// vectorised daxpyFMA routine, code taken from main loop (you can still see the old code in pre/post loops)
vmovdqu ymm0,ymmword ptr [r9+r13*8+10h]
vfmadd231pd ymm0,ymm1,ymmword ptr [rbx+r13*8+10h]
vmovdqu ymmword ptr [r9+r13*8+10h],ymm0
```

# Spliterator Characteristics and Performance

The streams API has been around for a while now, and I’m a big fan of it. It allows for a clean declarative programming style, which permits various optimisations to occur, and keeps the pastafarians at bay. I also think the `Stream` is the perfect abstraction for data interchange across API boundaries. This is partly because a `Stream` is lazy, meaning you don’t need to pay for consumption until you actually need to, and partly because a `Stream` can only be used once and there can be no ambiguity about ownership. If you supply a `Stream` to an API, you must expect that it has been used and so must discard it. This almost entirely eradicates defensive copies and can mean that no intermediate data structures need ever exist. Despite my enthusiasm for this abstraction, there’s some weirdness in this API when you scratch beneath surface.

I wanted to find a way to quickly run length encode an `IntStream` and found it difficult to make this code as fast as I’d like it to be. The code is too slow because it’s necessary to inspect each `int`, even when there is enough context available to potentially apply optimisations such as skipping over ranges. It’s likely that I am experiencing the friction of treating `Stream` as a data interchange format, which wasn’t one of its design goals, but this led me to investigate spliterator characteristics (there is no contiguous characteristic, which could speed up RLE greatly) and their relationship with performance.

### Spliterator Characteristics

Streams have spliterators, which control iteration and splitting behaviour. If you want to process a stream in parallel, it is the spliterator which dictates how to split the stream, if possible. There’s more to a spliterator than parallel execution though, and each single threaded execution can be optimised based on the characteristics bit mask. The different characteristics are as follows:

• `ORDERED` promises that there is an order. For instance, `trySplit` is guaranteed to give a prefix of elements.
• `DISTINCT` a promise that each element in the stream is unique.
• `SORTED` a promise that the stream is already sorted.
• `SIZED` promises the size of the stream is known. This is not true when a call to `iterate` generates the stream.
• `NONNULL` promises that no elements in the stream are null.
• `IMMUTABLE` promises the underlying data will not change.
• `CONCURRENT` promises that the underlying data can be modified concurrently. Must not also be `IMMUTABLE`.
• `SUBSIZED` promises that the sizes of splits are known, must also be `SIZED`.

There’s javadoc for all of these flags, which should be your point of reference, and you need to read it because you wouldn’t guess based on relative performance. For instance, `IntStream.range(inclusive, exclusive)` creates an `RangeIntSpliterator` with the characteristics `ORDERED | SIZED | SUBSIZED | IMMUTABLE | NONNULL | DISTINCT | SORTED`. This means that this stream has no duplicates, no nulls, is already sorted in natural order, the size is known, and it will be chunked deterministically. The data and the iteration order never change, and if we split it, we will always get the same first chunk. So these code snippets should have virtually the same performance:

``````
@Benchmark
public long countRange() {
return IntStream.range(0, size).count();
}

@Benchmark
public long countRangeDistinct() {
return IntStream.range(0, size).distinct().count();
}
``````

This is completely at odds with observations. Even though the elements are already distinct, and metadata exists to support this, requesting the distinct elements decimates performance.

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
countRange thrpt 1 10 49.465729 1.804123 ops/us 262144
countRangeDistinct thrpt 1 10 0.000395 0.000002 ops/us 262144

It turns out this is because `IntStream.distinct` has a one-size-fits-all implementation which completely ignores the `DISTINCT` characteristic, and goes ahead and boxes the entire range.

``````    // from java.util.stream.IntPipeline
@Override
public final IntStream distinct() {
// While functional and quick to implement, this approach is not very efficient.
// An efficient version requires an int-specific map/set implementation.
return boxed().distinct().mapToInt(i -> i);
}
``````

There is even more observable weirdness. If we wanted to calculate the sum of the first 1000 natural numbers, these two snippets should have the same performance. Requesting what should be a redundant sort doubles the throughput.

``````
@Benchmark
public long headSum() {
return IntStream.range(0, size).limit(1000).sum();
}

@Benchmark
public long sortedHeadSum() {
return IntStream.range(0, size).sorted().limit(1000).sum();
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
headSum thrpt 1 10 0.209763 0.002478 ops/us 262144
sortedHeadSum thrpt 1 10 0.584227 0.006004 ops/us 262144

In fact, you would have a hard time finding a relationship between Spliterator characteristics and performance, but you can see cases of characteristics driving optimisations if you look hard enough, such as in `IntStream.count`, where the `SIZED` characteristic is used.

``````    // see java.util.stream.ReduceOps.makeIntCounting
@Override
public <P_IN> Long evaluateSequential(PipelineHelper<Integer> helper, Spliterator<P_IN> spliterator) {
if (StreamOpFlag.SIZED.isKnown(helper.getStreamAndOpFlags()))
return spliterator.getExactSizeIfKnown();
return super.evaluateSequential(helper, spliterator);
}
``````

This is a measurably worthwhile optimisation, when benchmarked against the unsized spliterator created by `IntStream.iterate`:

``````
@Benchmark
public long countIterator() {
return IntStream.iterate(0, i -> i < size, i -> i + 1).count();
}

@Benchmark
public long countRange() {
return IntStream.range(0, size).count();
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
countIterator thrpt 1 10 0.001198 0.001629 ops/us 262144
countRange thrpt 1 10 43.166065 4.628715 ops/us 262144

What about `limit`, that’s supposed to be useful for speeding up streams by limiting the amount of work done? Unfortunately not. It actually makes things potentially much worse. In `SliceOps.flags`, we see it will actually disable `SIZED` operations:

``````    //see java.util.stream.SliceOps
private static int flags(long limit) {
return StreamOpFlag.NOT_SIZED | ((limit != -1) ? StreamOpFlag.IS_SHORT_CIRCUIT : 0);
}
``````

This has a significant effect on performance, as can be seen in the following benchmark:

``````
@Benchmark
public long countRange() {
return IntStream.range(0, size).count();
}

@Benchmark
public long countHalfRange() {
return IntStream.range(0, size).limit(size / 2).count();
}
``````

Benchmark Mode Threads Samples Score Score Error (99.9%) Unit Param: size
countHalfRange thrpt 1 10 0.003632 0.003363 ops/us 262144
countRange thrpt 1 10 44.859998 6.191411 ops/us 262144

It’s almost as if there were grand plans involving characteristic based optimisation, and perhaps time ran out (`IntStream.distinct` has a very apologetic comment) or others were better on paper than in reality. In any case, it looks like they aren’t as influential as you might expect. Given that the relationship between the characteristics which exist and performance is flaky at best, it’s unlikely that a new one would get implemented, but I think the characteristic `CONTIGUOUS` is missing.