@@ -167,6 +167,117 @@ def general_mm(a, b, c, M, N, K):
167167 return c
168168
169169
170+ @libentry ()
171+ @libtuner (
172+ configs = runtime .get_tuned_config ("mm_self_transpose" ),
173+ key = ["M" , "K" , "stride_am" , "stride_ak" ],
174+ strategy = ["align32" , "align32" , "align32" , "align32" ],
175+ warmup = 2 ,
176+ rep = 4 ,
177+ )
178+ @triton .jit
179+ def mm_kernel_syrk (
180+ A ,
181+ C ,
182+ M ,
183+ K ,
184+ stride_am ,
185+ stride_ak ,
186+ stride_cm ,
187+ stride_cn ,
188+ BLOCK_M : tl .constexpr ,
189+ BLOCK_K : tl .constexpr ,
190+ ):
191+ pid = tl .program_id (0 )
192+
193+ # Packed lower-triangular launch domain:
194+ # pid = row * (row + 1) / 2 + col, where 0 <= col <= row.
195+ #
196+ # Invert the triangular-number indexing by solving:
197+ # row^2 + row - 2 * pid = 0
198+ # => row = (-1 + sqrt(1 + 8 * pid)) / 2
199+ #
200+ # We take floor(...) as the candidate row, then apply an integer +/-1 correction
201+ # because fp32 sqrt can be off near triangular-number boundaries.
202+ pid_f = pid .to (tl .float32 )
203+ pid_m = tl .floor ((tl .sqrt (8.0 * pid_f + 1.0 ) - 1.0 ) / 2.0 ).to (tl .int32 )
204+ tri_start = pid_m * (pid_m + 1 ) // 2
205+ pid_m = tl .where (tri_start > pid , pid_m - 1 , pid_m )
206+ next_tri_start = (pid_m + 1 ) * (pid_m + 2 ) // 2
207+ pid_m = tl .where (next_tri_start <= pid , pid_m + 1 , pid_m )
208+ tri_start = pid_m * (pid_m + 1 ) // 2
209+ pid_n = pid - tri_start
210+
211+ rm = pid_m * BLOCK_M + tl .arange (0 , BLOCK_M )
212+ rn = pid_n * BLOCK_M + tl .arange (0 , BLOCK_M )
213+ ram = tl .max_contiguous (tl .multiple_of (rm % M , BLOCK_M ), BLOCK_M ).to (tl .int64 )
214+ ran = tl .max_contiguous (tl .multiple_of (rn % M , BLOCK_M ), BLOCK_M ).to (tl .int64 )
215+ rm = rm .to (tl .int64 )
216+ rn = rn .to (tl .int64 )
217+ acc = tl .zeros ((BLOCK_M , BLOCK_M ), dtype = tl .float32 )
218+
219+ for start_k in range (0 , K , BLOCK_K ):
220+ rk = (start_k + tl .arange (0 , BLOCK_K )).to (tl .int64 )
221+ mask_k = rk < K
222+ a = tl .load (
223+ A + (ram [:, None ] * stride_am + rk [None , :] * stride_ak ),
224+ mask = mask_k [None , :],
225+ other = 0.0 ,
226+ )
227+ b = tl .load (
228+ A + (rk [:, None ] * stride_ak + ran [None , :] * stride_am ),
229+ mask = mask_k [:, None ],
230+ other = 0.0 ,
231+ )
232+ acc += tl .dot (a , b , out_dtype = tl .float32 , allow_tf32 = False )
233+
234+ out = acc .to (C .dtype .element_ty )
235+ c_ptr = C + (rm [:, None ] * stride_cm + rn [None , :] * stride_cn )
236+ mask = (rm < M )[:, None ] & (rn < M )[None , :]
237+ tl .store (c_ptr , out , mask = mask )
238+
239+ if pid_m > pid_n :
240+ c_t_ptr = C + (rn [:, None ] * stride_cm + rm [None , :] * stride_cn )
241+ mask_t = (rn < M )[:, None ] & (rm < M )[None , :]
242+ tl .store (c_t_ptr , tl .trans (out ), mask = mask_t )
243+
244+
245+ def is_syrk_transpose_pair (a , b ):
246+ return (
247+ a .ndim == 2
248+ and b .ndim == 2
249+ and a .shape [0 ] == b .shape [1 ]
250+ and a .shape [1 ] == b .shape [0 ]
251+ and a .stride (0 ) == b .stride (1 )
252+ and a .stride (1 ) == b .stride (0 )
253+ and a .storage_offset () == b .storage_offset ()
254+ and a .data_ptr () == b .data_ptr ()
255+ )
256+
257+
258+ def syrk_mm (a , c , M , K ):
259+ grid = lambda META : (
260+ # Number of tile rows is tiles = ceil(M / BLOCK_M).
261+ # Packed lower triangle contains:
262+ # 1 + 2 + ... + tiles = tiles * (tiles + 1) / 2
263+ triton .cdiv (M , META ["BLOCK_M" ])
264+ * (triton .cdiv (M , META ["BLOCK_M" ]) + 1 )
265+ // 2 ,
266+ )
267+ with torch_device_fn .device (a .device ):
268+ mm_kernel_syrk [grid ](
269+ a ,
270+ c ,
271+ M ,
272+ K ,
273+ a .stride (0 ),
274+ a .stride (1 ),
275+ c .stride (0 ),
276+ c .stride (1 ),
277+ )
278+ return c
279+
280+
170281def streamk_scenario (a , b , M , N , K ):
171282 # TODO: this my change sometime according to the realbenchmark result
172283 # Currently, the best configuration for streamk has only been tested on A100(capability[0] == 8).
@@ -185,6 +296,10 @@ def streamk_scenario(a, b, M, N, K):
185296
186297def mm (a , b ):
187298 device = a .device
299+ if is_syrk_transpose_pair (a , b ):
300+ M , K = a .shape
301+ c = torch .empty ((M , M ), device = device , dtype = a .dtype )
302+ return syrk_mm (a , c , M , K )
188303 # handle non-contiguous inputs if necessary
189304 if a .stride (0 ) > 1 and a .stride (1 ) > 1 :
190305 a = a .contiguous ()
@@ -206,6 +321,9 @@ def mm(a, b):
206321
207322
208323def mm_out (a , b , * , out ):
324+ if is_syrk_transpose_pair (a , b ):
325+ M , K = a .shape
326+ return syrk_mm (a , out , M , K )
209327 # handle non-contiguous inputs if necessary
210328 if a .stride (0 ) > 1 and a .stride (1 ) > 1 :
211329 a = a .contiguous ()
0 commit comments