def matmul_cupydense(left, right, scale=1, out=None):
    """
    Perform the operation
        ``out := scale * (left @ right) + out``
    where `left`, `right` and `out` are matrices.  `scale` is a complex scalar,
    defaulting to 1.
    """
    # This may be done  more naturally with gemm from cupy.cublas
    # but the GPU timings which are not very different.
    if out is None:
        return (left @ right) * scale
    else:
        out += (left @ right) * scale
        return out