| 1 | #include "cpu_types.hpp" |
| 2 | |
| 3 | #include <torch/library.h> |
| 4 | |
| 5 | #include <cmath> |
| 6 | #include <cstdint> |
| 7 | #include <limits> |
| 8 | #include <mutex> |
| 9 | #include <random> |
| 10 | |
| 11 | namespace { |
| 12 | |
| 13 | constexpr int GUMBEL_TABLE_SIZE = 1 << 20; |
| 14 | constexpr int GUMBEL_TABLE_MASK = GUMBEL_TABLE_SIZE - 1; |
| 15 | |
| 16 | static float g_gumbel_table[GUMBEL_TABLE_SIZE]; |
| 17 | static std::once_flag g_gumbel_init_flag; |
| 18 | |
| 19 | static void init_gumbel_table() { |
| 20 | std::mt19937 rng(42); |
| 21 | std::uniform_real_distribution<float> uniform(1.0e-7f, 1.0f - 1.0e-7f); |
| 22 | for (int i = 0; i < GUMBEL_TABLE_SIZE; ++i) { |
| 23 | float u = uniform(rng); |
| 24 | g_gumbel_table[i] = -std::log(-std::log(u)); |
| 25 | } |
| 26 | } |
| 27 | |
| 28 | static inline void ensure_gumbel_table() { |
| 29 | std::call_once(g_gumbel_init_flag, init_gumbel_table); |
| 30 | } |
| 31 | |
| 32 | // Fused Gumbel-max kernel |
| 33 | |
| 34 | static void fused_gumbel_argmax_kernel(int64_t* __restrict__ output, |
| 35 | const float* __restrict__ logits, |
| 36 | const int64_t* __restrict__ seeds, |
| 37 | const int64_t batch_size, |
| 38 | const int64_t vocab_size) { |
| 39 | ensure_gumbel_table(); |
| 40 | constexpr int VEC_ELEM_NUM = vec_op::FP32Vec8::VEC_ELEM_NUM; |
| 41 | const int64_t vec_end = vocab_size - (vocab_size % VEC_ELEM_NUM); |
| 42 | |
| 43 | #pragma omp parallel for schedule(static) |
| 44 | for (int64_t b = 0; b < batch_size; ++b) { |
| 45 | const float* row = logits + b * vocab_size; |
| 46 | const uint32_t seed = static_cast<uint32_t>(seeds[b]); |
| 47 | |
| 48 | float best_score = -std::numeric_limits<float>::infinity(); |
| 49 | int64_t best_idx = 0; |
| 50 | |
| 51 | for (int64_t i = 0; i < vec_end; i += VEC_ELEM_NUM) { |
| 52 | vec_op::FP32Vec8 logit_vec(row + i); |
| 53 | |
| 54 | float gumbel_buf[VEC_ELEM_NUM]; |
| 55 | for (int k = 0; k < VEC_ELEM_NUM; ++k) |
| 56 | gumbel_buf[k] = g_gumbel_table[(seed + i + k) & GUMBEL_TABLE_MASK]; |
| 57 | vec_op::FP32Vec8 gumbel_vec(gumbel_buf); |
| 58 | |
| 59 | vec_op::FP32Vec8 score = logit_vec + gumbel_vec; |
| 60 | |
| 61 | float score_buf[VEC_ELEM_NUM]; |
| 62 | score.save(score_buf); |
| 63 | for (int k = 0; k < VEC_ELEM_NUM; ++k) { |
| 64 | if (score_buf[k] > best_score) { |
| 65 | best_score = score_buf[k]; |
| 66 | best_idx = i + k; |
| 67 | } |
| 68 | } |
| 69 | } |
| 70 | |
| 71 | for (int64_t i = vec_end; i < vocab_size; ++i) { |
| 72 | float s = row[i] + g_gumbel_table[(seed + i) & GUMBEL_TABLE_MASK]; |
| 73 | if (s > best_score) { |
| 74 | best_score = s; |
| 75 | best_idx = i; |
| 76 | } |
| 77 | } |
| 78 | |
| 79 | output[b] = best_idx; |
| 80 | } |
| 81 | } |
| 82 | |
| 83 | static void greedy_argmax_kernel(int64_t* __restrict__ output, |
| 84 | const float* __restrict__ logits, |
| 85 | const int64_t batch_size, |
| 86 | const int64_t vocab_size) { |
| 87 | constexpr int VEC_ELEM_NUM = vec_op::FP32Vec16::VEC_ELEM_NUM; |
| 88 | const int64_t vec_end = vocab_size - (vocab_size % VEC_ELEM_NUM); |
| 89 | |
| 90 | #pragma omp parallel for schedule(static) |
| 91 | for (int64_t b = 0; b < batch_size; ++b) { |
| 92 | const float* row = logits + b * vocab_size; |
| 93 | |
| 94 | vec_op::FP32Vec16 vmax(-std::numeric_limits<float>::infinity()); |
| 95 | for (int64_t i = 0; i < vec_end; i += VEC_ELEM_NUM) { |
| 96 | vmax = vmax.max(vec_op::FP32Vec16(row + i)); |
| 97 | } |
| 98 | float best_val = vmax.reduce_max(); |
| 99 | for (int64_t i = vec_end; i < vocab_size; ++i) { |
| 100 | if (row[i] > best_val) { |
| 101 | best_val = row[i]; |
| 102 | } |
| 103 | } |
| 104 | |
| 105 | int64_t best_idx = 0; |
| 106 | for (int64_t i = 0; i < vocab_size; ++i) { |
| 107 | if (row[i] == best_val) { |
| 108 | best_idx = i; |
| 109 | break; |
| 110 | } |
| 111 | } |
| 112 | output[b] = best_idx; |
| 113 | } |
| 114 | } |
| 115 | |
| 116 | } // namespace |
| 117 | |
| 118 | torch::Tensor fused_gumbel_argmax(const torch::Tensor& logits, |
| 119 | const torch::Tensor& seeds) { |
| 120 | TORCH_CHECK(logits.dim() == 2, "logits must be 2-D [batch, vocab]"); |
| 121 | TORCH_CHECK(logits.scalar_type() == torch::kFloat32, |
| 122 | "logits must be float32"); |
| 123 | TORCH_CHECK(seeds.dim() == 1 && seeds.size(0) == logits.size(0), |
| 124 | "seeds must be 1-D with batch_size elements"); |
| 125 | |
| 126 | auto logits_contig = logits.contiguous(); |
| 127 | auto output = torch::empty({logits_contig.size(0)}, torch::kInt64); |
| 128 | fused_gumbel_argmax_kernel( |
| 129 | output.data_ptr<int64_t>(), logits_contig.data_ptr<float>(), |
| 130 | seeds.data_ptr<int64_t>(), logits_contig.size(0), logits_contig.size(1)); |
| 131 | return output; |
| 132 | } |
| 133 | |
| 134 | torch::Tensor greedy_argmax(const torch::Tensor& logits) { |
| 135 | TORCH_CHECK(logits.dim() == 2, "logits must be 2-D [batch, vocab]"); |
| 136 | TORCH_CHECK(logits.scalar_type() == torch::kFloat32, |
| 137 | "logits must be float32"); |
| 138 | |
| 139 | // A single row is faster in ATen: OpenMP fork/join dominates the scan |
| 140 | // (~1 ms vs ~50 µs on s390x at vocab=32k). |
| 141 | if (logits.size(0) <= 1) { |
| 142 | return logits.argmax(/*dim=*/-1); |
| 143 | } |
| 144 | |
| 145 | auto logits_contig = logits.contiguous(); |
| 146 | auto output = torch::empty({logits_contig.size(0)}, torch::kInt64); |
| 147 | greedy_argmax_kernel(output.data_ptr<int64_t>(), |
| 148 | logits_contig.data_ptr<float>(), logits_contig.size(0), |
| 149 | logits_contig.size(1)); |
| 150 | return output; |
| 151 | } |