def METHOD_NAME(alpha: dc.float64, beta: dc.float64, C: dc.float64[NI, NJ], A: dc.float64[NI, NK], B: dc.float64[NK,
                                                                                                                NJ]):
    lA = np.empty((lNI, lNKa), dtype=A.dtype)
    lB = np.empty((lNKb, lNJ), dtype=B.dtype)
    lC = np.empty((lNI, lNJ), dtype=A.dtype)
    Av = np.reshape(A, (Px, lNI, Py, lNKa))
    A2 = np.transpose(Av, axes=(0, 2, 1, 3))
    Bv = np.reshape(B, (Px, lNKb, Py, lNJ))
    B2 = np.transpose(Bv, axes=(0, 2, 1, 3))
    Cv = np.reshape(C, (Px, lNI, Py, lNJ))
    C2 = np.transpose(Cv, axes=(0, 2, 1, 3))
    dc.comm.Scatter(A2, lA)
    dc.comm.Scatter(B2, lB)
    dc.comm.Scatter(C2, lC)
    tmp = distr.MatMult(lA, lB, (NI, NJ, NK))
    lC[:] = alpha * tmp + beta * lC
    dc.comm.Gather(lC, C2)
    C[:] = np.transpose(C2, (0, 2, 1, 3))