There are many examples for “matmul in cuda” (e.g. from NVidia but also others), you could write your own kernel. To get the derivative, I’d recommend writing a function that keeps the indices of the maxima (similar to torch.max with dimension) and use the corresponding lookup in the backward (you might be able to use “native” PyTorch for that.
Best regards
Thomas