Matrix multiplication SGEMM: Wider loads with register blocking

"WebGL2-compute" NxN matrix multiplication C = A x B (SGEMM) v.7 demo.
See Kernel 7: Wider loads with register blocking by Cedric Nugteren.
All A, B elements are random (0 - 1). Error "er1" is calculated as the sum of |CCPU - CGPU |/(N*N) for all matrix elements. er2 = max(|CCPU - CGPU |). See also Shader 7 benchmark.

Compute Shader 7

See also the page source
#version 310 es
#define WIDTH 4u                     // The vector-width (in number of floats)
#define TSM 128u                     // The tile-size in dimension M
#define TSN 128u                     // The tile-size in dimension N
#define TSK 16u                      // The tile-size in dimension K
#define WPTM 8u                      // The amount of work-per-thread in dimension M
#define WPTN 8u                      // The amount of work-per-thread in dimension N
#define LPTA ((TSK*WPTM*WPTN)/(TSN)) // The amount of loads-per-thread for A
#define LPTB ((TSK*WPTM*WPTN)/(TSM)) // The amount of loads-per-thread for B
#define RTSM 16u    // The reduced tile-size in dimension M (TSM/WPTM number of threads)
#define RTSN 16u    // The reduced tile-size in dimension N (TSN/WPTN number of threads)
#define MOD2(x,y) ((x) % (y))
#define DIV2(x,y) ((x) / (y))

layout (local_size_x = RTSM, local_size_y = RTSN, local_size_z = 1) in;
layout (std430, binding = 0) readonly buffer ssbA {
  vec4 A[];
};
layout (std430, binding = 1) readonly buffer ssbB {
  vec4 B[];
};
layout (std430, binding = 2) writeonly buffer ssbC {
  float C[];
};
  uniform uvec3 MNK;
  shared float Asub[TSK][TSM];    // Local memory to fit a tile of A and B
  shared float Bsub[TSK][TSN];
void main() {
    uint M = MNK.x, N = MNK.y, K = MNK.z;

    // Thread identifiers
    uint tidm = gl_LocalInvocationID.x; // Local row ID (max: TSM/WPTM == RTSM)
    uint tidn = gl_LocalInvocationID.y; // Local col ID (max: TSN/WPTN == RTSN)
    uint offsetM = TSM*gl_WorkGroupID.x; // Work-group offset
    uint offsetN = TSN*gl_WorkGroupID.y; // Work-group offset

    // Allocate register space
    float Areg;
    float Breg[WPTN];
    float acc[WPTM][WPTN];

    // Initialise the accumulation registers
    for (uint wm=0u; wm < WPTM; wm++) {
        for (uint wn=0u; wn < WPTN; wn++) {
            acc[wm][wn] = 0.0;
        }
    }
    // Loop over all tiles
    uint numTiles = K/TSK;
    uint t=0u;
    do {
        // Load one tile of A and B into local memory
        for (uint la=0u; la < LPTA/WIDTH; la++) {
            uint tid = tidn*RTSM + tidm;
            uint id = la*RTSN*RTSM + tid;
            uint row = MOD2(id,TSM/WIDTH);
            uint col = DIV2(id,TSM/WIDTH);

            // Load the values (wide vector load)
            uint tiledIndex = TSK*t + col;
            vec4 vecA = A[tiledIndex*(M/WIDTH) + offsetM/WIDTH + row];
            vec4 vecB = B[tiledIndex*(N/WIDTH) + offsetN/WIDTH + row];

            // Store the loaded vectors into local memory
            Asub[col][WIDTH*row + 0u] = vecA.x;
            Asub[col][WIDTH*row + 1u] = vecA.y;
            Asub[col][WIDTH*row + 2u] = vecA.z;
            Asub[col][WIDTH*row + 3u] = vecA.w;
            Bsub[col][WIDTH*row + 0u] = vecB.x;
            Bsub[col][WIDTH*row + 1u] = vecB.y;
            Bsub[col][WIDTH*row + 2u] = vecB.z;
            Bsub[col][WIDTH*row + 3u] = vecB.w;
        }
        // Synchronise to make sure the tile is loaded
        barrier();

        // Loop over the values of a single tile
        for (uint k=0u; k < TSK; k++) {

            // Cache the values of Bsub in registers
            for (uint wn=0u; wn < WPTN; wn++) {
                uint col = tidn + wn*RTSN;
                Breg[wn] = Bsub[k][col];
            }

            // Perform the computation
            for (uint wm=0u; wm < WPTM; wm++) {
                uint row = tidm + wm*RTSM;
                Areg = Asub[k][row];
                for (uint wn=0u; wn < WPTN; wn++) {
                    acc[wm][wn] += Areg * Breg[wn];
                }
            }
        }
        // Synchronise before loading the next tile
        barrier();

        // Next tile
        t++;
    } while (t < numTiles);

    // Store the final result in C
    for (uint wm=0u; wm < WPTM; wm++) {
        uint globalRow = offsetM + tidm + wm*RTSM;
        for (uint wn=0u; wn < WPTN; wn++) {
            uint globalCol = offsetN + tidn + wn*RTSN;
            C[globalCol*M + globalRow] = acc[wm][wn];
        }
    }
}
Simple transpose kernel for a PxQ matrix see at Shader v.5
SGEMM in WebGL2-compute     updated 7 Apr 2019