mention SIMD pospopcount really hard to comprehend

[libreriscv.git] / openpower / sv / cookbook / pospopcnt.mdwn

# Positional popcount SVP64

**Links**

* <https://bugs.libre-soc.org/show_bug.cgi?id=672>
* <https://github.com/clausecker/pospop/blob/master/countsse2_amd64.s>
* RISC-V Bitmanip Extension Document Version 0.94-draft Editor: Claire Wolf Symbiotic GmbH
<https://raw.githubusercontent.com/riscv/riscv-bitmanip/master/bitmanip-draft.pdf>

**Funding**

This project is funded through the [NGI Assure Fund](https://nlnet.nl/assure), a fund established by [NLnet](https://nlnet.nl) with financial support from the European Commission's [Next Generation Internet](https://ngi.eu) program. Learn more at the [NLnet project page](https://nlnet.nl/project/Libre-SOC-OpenPOWER-ISA).

[<img src="https://nlnet.nl/logo/banner.png" alt="NLnet foundation logo" width="20%" />](https://nlnet.nl)
[<img src="https://nlnet.nl/image/logos/NGIAssure_tag.svg" alt="NGI Assure Logo" width="20%" />](https://nlnet.nl/assure)

**Introduction**

Positional popcount in optimised assembler is typically done on SIMD
ISAs in around 500 lines.  The implementations are extraordinarily
complex and very hard to understand. Power ISA thanks to `bpermd` can be much
more efficient: with SVP64 even more so.  The reference implementation
showing the concept is below.

```
// Copyright (c) 2020 Robert Clausecker <fuz@fuz.su>
// count8 reference implementation for tests.  Do not alter.
func count8safe(counts *[8]int, buf []uint8) {
        for i := range buf {
                for j := 0; j < 8; j++ {
                        counts[j] += int(buf[i] >> j & 1)
                }
        }
}
```

A simple but still hardware-paralleliseable SVP64 assembler for
8-bit input values (`count8safe`) is as follows:

```
mtspr 9, 3                # move r3 to CTR
# VL = MIN(CTR,MAXVL=8), Rc=1 (CR0 set if CTR ends)
setvl 3,0,8,0,1,1         # set MVL=8, VL=MIN(MVL,CTR)
# load VL bytes (update r4 addr) but compressed (dw=8)
addi 6, 0, 0               # initialise all 64-bits of r6 to zero
sv.lbzu/pi/dw=8 *6, 1(4)   # should be /lf here as well
# gather performs the transpose (which gets us to positional..)
gbbd 8,6
# now those bits have been turned around, popcount and sum them
setvl 0,0,8,0,1,1          # set MVL=VL=8
sv.popcntd/sw=8 *24,*8     # do the (now transposed) popcount
sv.add *16,*16,*24         # and accumulate in results
# branch back if CTR still non-zero. works even though VL=8
sv.bc/all 16, *0, -0x28   # reduce CTR by VL and stop if -ve
```


Array popcount is just standard popcount function
([[!wikipedia Hamming weight]]) on an array of values, horizontally,
however positional popcount is different (vertical). Refer to Fig.1


<img src="/openpower/sv/cookbook/ArrayPopcnt.drawio.svg"
     alt="pospopcnt" width="70%" />

Positional popcount adds up the totals of each bit set to 1 in each
bit-position, of an array of input values. Refer to Fig.2

<img src="/openpower/sv/cookbook/PositionalPopcnt.drawio.svg"
     alt="pospopcnt" width="60%" />


# Visual representation of the pospopcount algorithm

In order to perform positional popcount we need to go through series of
steps shown below in figures 3, 4, 5 & 6.  The limitation to overcome is
that the CPU can only work on registers (horizontally) but the requirement
of pospopcount is to add *vertically*.  Part of the challenge is therefore
to perform an appropriate transpose of the data, in blocks that suit
the processor and the ISA capacity.

Fig.3 depicts how the data is divided into blocks of 8*8-bit.
 
<img src="/openpower/sv/cookbook/BlockDivision.drawio.svg"
     alt="pospopcnt" width="70%" />

Fig.4 shows that a transpose is needed on each block.  The gbbd
instruction is used for implementing the transpose function, preparing
the data for using the standard popcount instruction.

<img src="/openpower/sv/cookbook/Transpose.drawio.svg"
     alt="pospopcnt" width="60%" />

Now on each 8*8 block, 8 popcount instructions can be run each of
which is independent and therefore can be parallelised even by In-Order
multi-issue hardware(Fig.5).

<img src="/openpower/sv/cookbook/PopcountBlocks.drawio.svg"
     alt="pospopcnt" width="70%" />

Fig.6 depicts how each of the intermediate results are accumulated. It
is worth noting that each intermediate result is independent of the
other intermediate results and also parallel reduction can be applied
to all of them individually. This gives *two* opportunities for hardware
parallelism rather than one.

<img src="/openpower/sv/cookbook/ParallelAccumulate.drawio.svg"
     alt="pospopcnt" width="100%" />

In short this algorithm is very straightforward to implement thanks to the
two crucial instructions, `gbbd` and `popcntd`.  Below is a walkthrough
of the assembler, keeping it very simple, and exploiting only one of the
opportunities for parallelism (by not including the Parallel Reduction
opportunity mentioned above).

# Walkthrough of the assembler

Firstly the CTR (Counter) SPR is set up, and is key to looping
as outlined further, below

```
mtspr 9, 3                # move r3 to CTR
```

The Vector Length, which is limited to 8 (MVL - Maximum Vector Length)
is set up. A special "CTR" Mode is used which automatically uses the CTR
SPR rather than register RA. (*Note that RA is set to zero to indicate
this, because there is limited encoding space. See [[openpower/sv/setvl]]
instruction specification for details)*.

The result of this instruction is that if CTR is greater than 8, VL is set
to 8. If however CTR is less than or equal to 8, then VL is set to CTR.
Additionally, a copy of VL is placed into RT (r3 in this case), which
again is necessary as part of the limited encoding space but in some cases
(not here) this is desirable, and avoids a `mfspr` instruction to take
a copy of VL into a GPR.

```
# VL = r3 = MIN(CTR,MAXVL=8)
setvl 3,0,8,0,1,1         # set MVL=8, VL=MIN(MVL,CTR)
```

Next the data must be loaded in batches (of up to QTY 8of 8-bit
values). We must also not over-run the end of the array (which a normal
SIMD ISA would do, potentially causing a pagefault) and this is achieved
by when CTR <= 8: VL (Vector Length) in such circumstances being set to
the remaining elements.

The next instruction is `sv.lbzu/pi` - a Post-Increment variant of `lbzu`
which updates the Effective Address (in RA) **after** it has been used,
not before. Also of note is the element-width override of `dw=8` which
writes to *individual bytes* of r6, not to r6/7/8...13.

Of note here however is the precursor `addi` instruction, which clears
out the contents of r6 to zero before beginning the Vector/Looped Load
sequence.  This is *important* because SV does **not** zero-out parts
of a register when element-width overrides are in use.

```
# load VL bytes (update r4 addr) but compressed (dw=8)
addi 6, 0, 0               # initialise all 64-bits of r6 to zero
sv.lbzu/pi/dw=8 *6, 1(4)   # should be /lf here as well
```

The next instruction is quite straightforward, and has the effect of
turning (transposing) 8 values.  Each bit zero of each input byte is
placed into the first byte; each bit one of each input byte is placed
into the second byte and so on.

(*This instruction in VSX is called `vgbbd`, and it is present in many
other ISAs. RISC-V bitmanip-draft-0.94 calls it `bmatflip`, x86 calls it
GF2P7AFFINEQB. Power ISA, Cray, and many others all have this extremely
powerful and useful instruction*)

```
# gather performs the transpose (which gets us to positional..)
gbbd 8,6
```

Now the bits have been transposed in bytes, it is plain sailing to
perform QTY8 parallel popcounts.  However there is a trick going on:
we have set VL=MVL=8. To explain this: it covers the last case where
CTR may be between 1 and 8.

Remember what happened back at the Vector-Load, where r6 was cleared to
zero before-hand? This filled out the 8x8 transposed grid (`gbbd`) fully
with zeros prior to the actual transpose.  Now when we do the popcount,
there will be upper-numbered columns that are *not part of the result*
that contain *zero* and *consequently have no impact on the result*.

This elegant trick extends even to the accumulation of the
results. However before we get there it is necessary to describe why
`sw=8` has been added to `sv.popcntd`. What this is doing is treating
each **byte** of its input (starting at the first byte of r8) as an
independent quantity to be popcounted, but the result is *zero-extended*
to 64-bit and stored in r24, r25, r26 ... r31.

Therefore:

* r24 contains the popcount of the first byte of r8
* r25 contains the popcount of the second byte of r8
* ...
* r31 contains the popcount of the last (7th) byte of r8

```
# now those bits have been turned around, popcount and sum them
setvl 0,0,8,0,1,1          # set MVL=VL=8
sv.popcntd/sw=8 *24,*8     # do the (now transposed) popcount
```

With VL being hard-set to 8, we can now Vector-Parallel sum the partial
results r24-31 into the array of accumulators r16-r23. There was no need
to set VL=8 again because it was set already by the `setvl` above.

```
sv.add *16,*16,*24         # and accumulate in results
```

Finally, `sv.bc` loops back, subtracting VL from CTR.  Being based on
Power ISA Scalar Branch-Conditional, if CTR goes negative (which it will
because VL was set hard-coded to 8, above) it does not matter, the loop
will still be correctly exited. In other words we complete the correct
number of *blocks* (of length 8).

```
# branch back if CTR still non-zero. works even though VL=8
sv.bc/all 16, *0, -0x28   # reduce CTR by VL and stop if -ve
```

# Improvements

There exist many opportunities for parallelism that simpler hardware would
need to have in order to maximise performance.  On Out-of-Order hardware
the above extremely simple and very clear algorithm will achieve extreme
high levels of performance simply by exploiting standard Multi-Issue
Register Hazard Management.

However simpler hardware - in-order - will need a little bit of
assistance, and that can come in the form of expanding to QTY4 or QTY8
64-bit blocks (so that sv.popcntd uses MVL=VL=32 or MVL=VL=64), `gbbd`
becomes an `sv.gbbd` but VL being set to the block count (QTY4 or QTY8),
and the SV REMAP Parallel Reduction Schedule being applied to each
intermediary result rather than using an array of straight accumulators
`r16-r23`.  However this starts to push the boundaries of the number of
registers needed, so as an exercise is left for another time.

# Conclusion

Where a normal SIMD ISA requires explicit hand-crafted optimisation in
order to achieve full utilisation of the underlying hardware, Simple-V
instead can rely to a large extent on standard Multi-Issue hardware
to achieve similar performance, whilst crucially keeping the algorithm
implementation down to a shockingly-simple degree that makes it easy to
understand an easy to review.  Again also as with many other algorithms
when implemented in Simple-V SVP54, by keeping to a LOAD-COMPUTE-STORE
paradigm the L1 Data Cache usage is minimised, and in this case just
as with chacha20 the entire algorithm, being only 9 lines of assembler
fitting into 13 4-byte words it can fit into a single L1 I-Cache Line
without triggering Virtual Memory TLB misses.

Further performance improvements are achievable by using REMAP Parallel
Reduction, still fitting into a single L1 Cache line, but beginning to
approach the limit of the 128-long register file.

[[!tag svp64_cookbook ]]