Matrix multiplication SGEMM: More work per thread

"WebGL 2 compute" NxN matrix multiplication C = A x B (SGEMM) v.3 demo.
See Kernel 3: More work per thread 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 3 benchmark.

Compute Shader 3

See also the page source
#version 310 es
#define TS 32u
#define WPT 8u
#define RTS 4u  // TS/WPT
layout (local_size_x = TS, local_size_y = RTS, local_size_z = 1) in;
layout (std430, binding = 0) readonly buffer ssbA {
  float A[];
};
layout (std430, binding = 1) readonly buffer ssbB {
  float B[];
};
layout (std430, binding = 2) writeonly buffer ssbC {
  float C[];
};
  uniform uvec3 MNK;
  shared float Asub[TS][TS];  // Local memory to fit a tile of
  shared float Bsub[TS][TS];  // TS*TS elements of A and B
void main() {
    uint M = MNK.x, N = MNK.y, K = MNK.z;

    // Thread identifiers
    uint row = gl_LocalInvocationID.x; // Local row ID (max: TS)
    uint col = gl_LocalInvocationID.y; // Local col ID (max: TS/WPT == RTS)
    uint globalRow = TS*gl_WorkGroupID.x + row; // Row ID of C (0..M)
    uint globalCol = TS*gl_WorkGroupID.y + col; // Col ID of C (0..N)

    // Initialise the accumulation registers
    float acc[WPT];
    for (uint w=0u; w < WPT; w++) {
        acc[w] = 0.0;
    }

    // Loop over all tiles
    uint numTiles = K/TS;
    for (uint t=0u; t < numTiles; t++) {

        // Load one tile of A and B into local memory
        for (uint w=0u; w < WPT; w++) {
            uint tiledRow = TS*t + row;
            uint tiledCol = TS*t + col;
            Asub[col + w*RTS][row] = A[(tiledCol + w*RTS)*M + globalRow];
            Bsub[col + w*RTS][row] = B[(globalCol + w*RTS)*K + tiledRow];
        }

        // Synchronise to make sure the tile is loaded
        memoryBarrierShared();
        barrier();

        // Perform the computation for a single tile
        for (uint k=0u; k < TS; k++) {
            for (uint w=0u; w < WPT; w++) {
                acc[w] += Asub[k][row] * Bsub[col + w*RTS][k];
            }
        }

        // Synchronise before loading the next tile
        barrier();
    }
    // Store the final result in C
    for (uint w=0u; w < WPT; w++) {
        C[(globalCol + w*RTS)*M + globalRow] = acc[w];
    }
}

Comments:

Thanks to Kentaro Kawakatsu for help with shared memory synchronization. Compute: Control barrier and shared memory within the local work group
SGEMM in WebGL2-compute     updated 1 Mar 2019