|
| 1 | +import logging |
| 2 | + |
| 3 | +import torch |
| 4 | +import triton |
| 5 | +import triton.language as tl |
| 6 | + |
| 7 | +logger = logging.getLogger("flag_gems").getChild(__name__.lstrip(".")) |
| 8 | + |
| 9 | + |
| 10 | +@triton.jit |
| 11 | +def _digamma(x): |
| 12 | + PI = 3.141592653589793 |
| 13 | + |
| 14 | + needs_reflect = x <= 0.0 |
| 15 | + xp = tl.where(needs_reflect, 1.0 - x, x) |
| 16 | + |
| 17 | + t = tl.minimum(xp, 1.0e19) |
| 18 | + r = xp + 3.0 |
| 19 | + inv_sum = (3.0 * t * t + 6.0 * t + 2.0) / (t * (t + 1.0) * (t + 2.0)) |
| 20 | + |
| 21 | + rinv = 1.0 / r |
| 22 | + z = rinv * rinv |
| 23 | + poly = z * ( |
| 24 | + 1.0 / 12.0 + z * (-1.0 / 120.0 + z * (1.0 / 252.0 + z * (-1.0 / 240.0))) |
| 25 | + ) |
| 26 | + psi = tl.math.log(r) - 0.5 * rinv - poly - inv_sum |
| 27 | + |
| 28 | + cot = tl.math.cos(PI * x) / tl.math.sin(PI * x) |
| 29 | + psi = tl.where(needs_reflect, psi - PI * cot, psi) |
| 30 | + return psi |
| 31 | + |
| 32 | + |
| 33 | +@triton.jit |
| 34 | +def digamma_kernel(in_ptr, out_ptr, n_elements, BLOCK_SIZE: tl.constexpr): |
| 35 | + pid = tl.program_id(0) |
| 36 | + offsets = pid * BLOCK_SIZE + tl.arange(0, BLOCK_SIZE) |
| 37 | + mask = offsets < n_elements |
| 38 | + x = tl.load(in_ptr + offsets, mask=mask, other=1.0) |
| 39 | + y = _digamma(x.to(tl.float32)) |
| 40 | + tl.store(out_ptr + offsets, y, mask=mask) |
| 41 | + |
| 42 | + |
| 43 | +def digamma(input): |
| 44 | + logger.debug("GEMS_KUNLUNXIN DIGAMMA") |
| 45 | + output = torch.empty_like(input) |
| 46 | + n = input.numel() |
| 47 | + if n == 0: |
| 48 | + return output |
| 49 | + if n >= (1 << 22): |
| 50 | + if input.dtype == torch.bfloat16: |
| 51 | + BLOCK_SIZE, NUM_WARPS = 16384, 16 |
| 52 | + else: |
| 53 | + BLOCK_SIZE, NUM_WARPS = 8192, 8 |
| 54 | + else: |
| 55 | + BLOCK_SIZE, NUM_WARPS = 1024, 4 |
| 56 | + grid = (triton.cdiv(n, BLOCK_SIZE),) |
| 57 | + digamma_kernel[grid](input, output, n, BLOCK_SIZE=BLOCK_SIZE, num_warps=NUM_WARPS) |
| 58 | + return output |
0 commit comments