[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>

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.
4. Move to the next row of the first matrix, and next column of the second matrix.

This for-loop uses the indices as shown above

```
    for i in range(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
```


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
```

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.

# 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 (OpenPOWER 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 SHAPE Special
Purpose Registers (SPRs), which were added alongside the SVP64 looping
system for complex element indexing. Without having SHAPE SPRs, only the most
basic, consecuting indexing of register elements (0,1,2,3...) would
be possible.

### SHAPE Remapping SPRs

* See [[sv/remap]] for the full break down of SPRs SHAPE0-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.

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 SHAPE0 SPR)
- 

## SVREMAP


## Appendix