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

# SVP64 REMAP Worked Example: Matrix Multiply

Links

* [Online matrix calculator](https://matrix.reshish.com/multCalculation.php)
* [[sv/remap]]
* <https://bugs.libre-soc.org/show_bug.cgi?id=701>
* video <https://m.youtube.com/watch?v=BbhOA8ULKt4> slides
  <https://ftp.libre-soc.org/openpower_2021.pdf>
* setvl instruction: [[sv/setvl]]
* REMAP python function based on yield (shows how indices are generated):
[[sv/remapyield.py]]

TODO: Include screenshots

One of the most powerful and versatile modes of the REMAP engine (a part of
the SVP64 feature set) is the ability to perform matrix
multiplication with all elements within a scalar register file.

This is done by converting the index used to iterate over the operand and
result matrices to the actual index inside the scalar register file.

## Worked example - manual (conventional method)

The matrix multiply looks like this:

```
    mat_X * mat_Y = mat_Z
```

When multiplying non-square matrices (rows != columns), to determine the
dimension of the result when matrix X has `a` rows and `b` columns and matrix
Y has `b` rows and `c` columns:

```
    X_axb * Y_bxc = Z_axc
```

The result matrix will have number of rows of the first matrix, and number of
columns of the second matrix.


For this example, the following values will be used for the operand matrices X
and Y, result Z shown for completeness.

    X =| 1 2 3 |  Y =  |  6  7 |  Z = |  52  58 |
       | 3 4 5 |       |  8  9 |      | 100 112 |
                       | 10 11 |

Matrix X has 2 rows, 3 columns (2x3), and matrix Y has 3 rows, 2 columns.

To determine the final dimensions of the resultant matrix Z, take the number
of rows from matrix X (2) and number of columns from matrix Y (2).

The method usually taught in linear algebra course to students is the
following (outer product):

1. Start with the first row of the first matrix, and first column of the
second matrix.
2. Multiply each element in the row by each element in the column, and sum
with the current value of the result matrix element (multiply-add-accumulate).
Store result in the first row, first column entry.
3. Move to the next column of the second matrix, and next column of the result
matrix. If there are no more columns in the second matrix, go back to first
column (second matrix), and move to next row (first and result matrices).
If there are no more rows left, result matrix has been fully computed.
4. Repeat step 2.

This for-loop uses the indices as shown above

```
    for i in range(0, mat_X_num_rows):
        for k in range(0, mat_Y_num_cols):
            for j in range(0, mat_X_num_cols): # or mat_Y_num_rows
                mat_Z[i][k] += mat_X[i][j] * mat_Y[j][k]
```

Calculations:

```
    | 1 2 3 |   |  6  7 | = | (1*6 + 2*8 + 3*10) (1*7 + 2*9 3*11) |
    | 3 4 5 | * |  8  9 |   | (3*6 + 4*8 + 5*10) (3*7 + 4*9 5*11) |
                | 10 11 |
                
    | 1 2 3 |   |  6  7 | = | ( 6 + 16 + 30) ( 7 + 18 + 33) |
    | 3 4 5 | * |  8  9 |   | (18 + 32 + 50) (21 + 36 + 55) |
                | 10 11 |
                
    | 1 2 3 |   |  6  7 | = |  52  58 |
    | 3 4 5 | * |  8  9 |   | 100 112 |
                | 10 11 |
```

For the algorithm, assign indeces to matrices as follows:

```
    Index | 0 1 2 3 4 5 |
    Mat X | 1 2 3 3 4 5 |

    Index | 0 1 2 3  4  5 |
    Mat Y | 6 7 8 9 10 11 |

    Index |  0  1   2   3 |
    Mat Z | 52 58 100 112 |
```

(Start with the first row, then assign index left-to-right, top-to-bottom.)

Index list:

```
    | Mat X | Mat Y | Mat Z |
    |   0   |   0   |   0   |
    |   1   |   2   |   0   |
    |   2   |   4   |   0   |
    |   0   |   1   |   1   |
    |   1   |   3   |   1   |
    |   2   |   5   |   1   |
    |   3   |   0   |   2   |
    |   4   |   2   |   2   |
    |   5   |   4   |   2   |
    |   3   |   1   |   3   |
    |   4   |   3   |   3   |
    |   5   |   5   |   3   |
```

Worked example broken down into individual multiply-add accumulates:

[[!img outer_product_worked_example.jpg size="600x" ]]

The issue with this algorithm is that the result matrix element is the same
for three consecutive operations, and where each element is stored in CPU
registers, the same register will be written to three times and thus causing
consistent stalling.

## Inner Product

A slight modification to the order of the loops in the algorithm massively
reduces the chance of read-after-write hazards, as the result matrix
element (and thus register) changes with every multiply-add operation.

The code:

```
    for i in range(mat_X_num_rows):
        for j in range(0, mat_X_num_cols): # or mat_Y_num_rows
            for k in range(0, mat_Y_num_cols):
                mat_Z[i][k] += mat_X[i][j] * mat_Y[j][k]
```

Index list:

```
    | Mat X | Mat Y | Mat Z |
    |   0   |   0   |   0   |
    |   0   |   1   |   1   |
    |   3   |   0   |   2   |
    |   3   |   1   |   3   |
    |   1   |   2   |   0   |
    |   1   |   3   |   1   |
    |   4   |   2   |   2   |
    |   4   |   3   |   3   |
    |   2   |   4   |   0   |
    |   2   |   5   |   1   |
    |   5   |   4   |   2   |
    |   5   |   5   |   3   |
```

Worked example for inner product:

[[!img inner_product_worked_example.jpg size="600x" ]]

The index for the result matrix changes with every operation, and thus the
consecutive multiply-add instruction doesn't depend on the previous write
register.

Outer and inner product indices side-by-side:

```
    |   Outer Product       |   Inner Product       |
    | Mat X | Mat Y | Mat Z | Mat X | Mat Y | Mat Z |
    |   0   |   0   |   0   |   0   |   0   |   0   |
    |   1   |   2   |   0   |   0   |   1   |   1   |
    |   2   |   4   |   0   |   3   |   0   |   2   |
    |   0   |   1   |   1   |   3   |   1   |   3   |
    |   1   |   3   |   1   |   1   |   2   |   0   |
    |   2   |   5   |   1   |   1   |   3   |   1   |
    |   3   |   0   |   2   |   4   |   2   |   2   |
    |   4   |   2   |   2   |   4   |   3   |   3   |
    |   5   |   4   |   2   |   2   |   4   |   0   |
    |   3   |   1   |   3   |   2   |   5   |   1   |
    |   4   |   3   |   3   |   5   |   4   |   2   |
    |   5   |   5   |   3   |   5   |   5   |   3   |
```

# SVP64 instructions implementing matrix multiply

* SVP64 assembler example:
[unit test](https://git.libre-soc.org/?p=openpower-isa.git;a=blob;f=src/openpower/decoder/isa/test_caller_svp64_matrix.py;hb=30f2d8a8e92ad2939775f19e6a0f387499e9842b#l56)
* svremap and svshape instructions defined in:
[[sv/rfc/ls009]
* Multiple-Add Low Doubleword instruction pseudo-code (Power ISA 3.0C
Book I, section 3.3.9): [[openpower/isa/fixedarith]]

*(Need to check if first arg of svremap correct, then one shown works with
ISACaller)*

```
    svshape 2, 2, 3, 0, 0
    svremap 31, 1, 2, 3, 0, 0, 0
    sv.maddld *0, *16, *32, *0
```

## svshape

The `svshape` instruction is a convenient way to access the `SVSHAPE` Special
Purpose Registers (SPRs), which were added alongside the SVP64 looping
system for complex element indexing. Without having "Re-shaping" SPRs, only the most
basic, consecuting indexing of register elements (0,1,2,3...) would
be possible.

### SVSHAPE Remapping SPRs

* See [[sv/remap]] for the full break down of SPRs `SVSHAPE0-3`.

For Matrix Multiply, SHAPE0 SPR is used:

|0:5   |6:11  | 12:17   | 18:20   | 21:23   |24:27 |28:29  |30:31|
|----- |----- | ------- | ------- | ------  |------|------ |---- |
|xdimsz|ydimsz| zdimsz  | permute | invxyz  |offset|skip   |mode |


skip:

- 0b00 indicates no dimensions to be skipped
- 0b01 - skip '1st dim'
- 0b10 - skip '2nd dim'
- 0b11 - skip '3rd dim'

invxyz (3-bit; 1 for x, 1 for y, 1 for z):

- If corresponding dim bit is zero, start index from zero and increment
- If bit set, start from xdimsz-1 (x dimension size, or whichever dimension
bit is being looked at) and decrement down to zero.

offset is used to offset the result by `offset` elements (important for when
using element width overrides are used).

xdimsz, ydimsz, zdimsz are offset by 1, such that 0-0b111111 correspond to
1-64. A value of xdimsz=2 would indicate that in the first dimension there are
3 elements in the array.

With the example Matrix X (2 rows, 3 columns, or 2x3 matrix), xdimsz=1,
ydimsz=2, zdimsz=0.

permute setting:

| permute | order | array format             |
| ------- | ----- | ------------------------ |
| 000     | 0,1,2 | (xdim+1)(ydim+1)(zdim+1) |
| 001     | 0,2,1 | (xdim+1)(zdim+1)(ydim+1) |
| 010     | 1,0,2 | (ydim+1)(xdim+1)(zdim+1) |
| 011     | 1,2,0 | (ydim+1)(zdim+1)(xdim+1) |
| 100     | 2,0,1 | (zdim+1)(xdim+1)(ydim+1) |
| 101     | 2,1,0 | (zdim+1)(ydim+1)(xdim+1) |
| 110     | 0,1   | Indexed (xdim+1)(ydim+1) |
| 111     | 1,0   | Indexed (ydim+1)(xdim+1) |

Permute re-arranges the order of the nested for-loops used to iterate over the
three dimensions. This allows for in-place transpose, in-place rotate, matrix
multiply, convolutions, without the limitation of Power-of-Two matrices.

For normal matrix multiply, the permute setting is 0b010 (order 1,0,2,
or swap x and y loops).

(*NOTE:* This is done automatically by the Matrix-Multiply REMAP mode, `SVRM=0`.)

Limitations of Matrix REMAP are currently:

- Vector Length (VL) limited to 127, and up to 127 Multiply-Add Accumulates
(MAC), or other operations may be performed in total.
For matrix multiply, it means both operand matrices and result matrix can
have no more than 127 elements in total.
(Larger matrices can be split into tiles to circumvent this issue, out
of scope of this document).
- `svshape` instruction only provides part of the Matrix REMAP capability.
For rotation and mirroring, `SVSHAPE` SPRs must be programmed directly (thus
requiring more assembler instructions). Future revisions of SVP64 will
provide more comprehensive capacity, mitigating the need to write to `SVSHAPE`
SPRs directly.

Going back to the assembler instruction used to setup the shape for matrix
multiply:

```
    svshape 2, 2, 3, 0, 0
```

breakdown:

- SVxd=2, SVyd=2, SVzd=3
- SVRM=0 (Matrix mode, uses `SVSHAPE0` SPR)
- vf=0 (not using Vertical-First mode)

To determine the `SVxd`/`SVyd`/`SVzd` settings:

- `SVxd` is equal to the number of columns in the second operand matrix.
Matrix Y has 2 columns, so `SVxd=2`.
- `SVyd` is equal to the number of rows in the first operand matrix.
Matrix X has 2 rows, so `SVyd=2`.
- `SVzd` is equal to the number of columns in the first operand matrix,
or the number of rows in the second operand matrix.
Matrix X has 3 columns, matrix Y has 3 rows, so `SVzd=3`.

Table form

```
    SVxd |         mat_Y_num_cols
    SVyd |         mat_X_num_rows
    SVzd | mat_X_num_cols OR mat_Y_num_rows
```

## SVREMAP

```
    svremap 31, 1, 2, 3, 0, 0, 0
```

Assigns the configured SVSHAPEs to the relevant operand/result registers
of the consecutive instruction/s (depending on if REMAP is set to persistent).

## maddld - Multiply-Add Low Doubleword VA-form

```
    sv.maddld *0, *16, *32, *0
```

A standard instruction available since version 3.0 of PowerISA.

*Temporary note:* maddld (Multiply-Add Low Doubleword) in the 3.1b version
of the PowerISA spec is in the Linux Compliancy Subset, not SFS or SFFS.
See page 1477 of the document, or page 1503 of the pdf.

This instruction can be used as a multiply-add accumulate by setting the
third operand to be the same as the result register, which functions as
an accumulator.

## Appendix


[[!tag svp64_cookbook ]]