chore: remove PR1 agent comments, trim kernel docs, add PR4 benchmark

Author: aloha1357

qdp/qdp-kernels/src/amplitude.cu

@@ -41,27 +41,18 @@ __global__ void amplitude_encode_kernel(
     double v1 = 0.0;
     double v2 = 0.0;
 
-    // Vectorized Load Optimization:
-    // If we are well within bounds, treat input as double2 to issue a single 128-bit load instruction.
-    // Use __ldg() to pull through the read-only cache; cudaMalloc aligns to 256 bytes so the
-    // reinterpret_cast<double2*> load is naturally aligned.
+    // double2 load via __ldg when aligned and in bounds.
     if (state_idx_base + 1 < input_len) {
-        // Reinterpret cast to load two doubles at once
         const double2 loaded = __ldg(reinterpret_cast<const double2*>(input) + idx);
         v1 = loaded.x;
         v2 = loaded.y;
     }
-    // Handle edge case: Odd input length
     else if (state_idx_base < input_len) {
         v1 = __ldg(input + state_idx_base);
-        // v2 remains 0.0
     }
 
-    // Write output:
-    // Apply pre-calculated reciprocal (multiplication is faster than division)
     state[state_idx_base]     = make_cuDoubleComplex(v1 * inv_norm, 0.0);
 
-    // Check boundary for the second element (state_len is usually power of 2, but good to be safe)
     if (state_idx_base + 1 < state_len) {
         state[state_idx_base + 1] = make_cuDoubleComplex(v2 * inv_norm, 0.0);
     }
@@ -82,7 +73,6 @@ __global__ void amplitude_encode_kernel_f32(
     float v2 = 0.0f;
 
     if (state_idx_base + 1 < input_len) {
-        // Mirror the double kernel: cached vectorized load for two floats
         const float2 loaded = __ldg(reinterpret_cast<const float2*>(input) + idx);
         v1 = loaded.x;
         v2 = loaded.y;
@@ -225,17 +215,7 @@ int launch_amplitude_encode_f32(
     return (int)cudaGetLastError();
 }
 
-/// Optimized batch amplitude encoding kernel
-///
-/// Memory Layout (row-major):
-/// - input_batch: [sample0_data | sample1_data | ... | sampleN_data]
-/// - state_batch: [sample0_state | sample1_state | ... | sampleN_state]
-///
-/// Optimizations:
-/// 1. Vectorized double2 loads for 128-bit memory transactions when aligned
-/// 2. Grid-stride loop for arbitrary batch sizes
-/// 3. Coalesced memory access within warps
-/// 4. Scalar fallback for misaligned sample bases and odd tails
+/// Batch amplitude encoding kernel (grid-stride, vectorized loads when aligned).
 __global__ void amplitude_encode_batch_kernel(
     const double* __restrict__ input_batch,
     cuDoubleComplex* __restrict__ state_batch,
@@ -244,25 +224,20 @@ __global__ void amplitude_encode_batch_kernel(
     size_t input_len,
     size_t state_len
 ) {
-    // Grid-stride loop pattern for flexibility
-    const size_t elements_per_sample = state_len / 2;  // Each thread handles 2 elements
+    const size_t elements_per_sample = state_len / 2;
     const size_t total_work = num_samples * elements_per_sample;
     const size_t stride = gridDim.x * blockDim.x;
 
     size_t global_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
-    // Process elements in grid-stride fashion
     for (size_t idx = global_idx; idx < total_work; idx += stride) {
-        // Decompose linear index into (sample, element_pair)
         const size_t sample_idx = idx / elements_per_sample;
         const size_t elem_pair = idx % elements_per_sample;
 
-        // Calculate base addresses (strength-reduced)
         const size_t input_base = sample_idx * input_len;
         const size_t state_base = sample_idx * state_len;
         const size_t elem_offset = elem_pair * 2;
 
-        // Load inverse norm (cached by L1)
         const double inv_norm = inv_norms[sample_idx];
 
         double v1, v2;
@@ -281,34 +256,20 @@ __global__ void amplitude_encode_batch_kernel(
                 ? __ldg(sample_input + elem_offset + 1)
                 : 0.0;
         } else {
-            // Padding region
             v1 = v2 = 0.0;
         }
 
-        // Normalize and write as complex numbers
-        // Compiler will optimize multiplications
         const cuDoubleComplex c1 = make_cuDoubleComplex(v1 * inv_norm, 0.0);
         const cuDoubleComplex c2 = make_cuDoubleComplex(v2 * inv_norm, 0.0);
 
-        // Write to global memory (coalesced within warp)
         state_batch[state_base + elem_offset] = c1;
         if (elem_offset + 1 < state_len) {
             state_batch[state_base + elem_offset + 1] = c2;
         }
     }
 }
 
-/// Optimized batch amplitude encoding kernel (float32)
-///
-/// Memory Layout (row-major):
-/// - input_batch: [sample0_data | sample1_data | ... | sampleN_data]
-/// - state_batch: [sample0_state | sample1_state | ... | sampleN_state]
-///
-/// Optimizations:
-/// 1. Vectorized float2 loads for 64-bit memory transactions
-/// 2. Grid-stride loop for arbitrary batch sizes
-/// 3. Coalesced memory access within warps
-/// 4. Minimized register pressure
+/// Batch amplitude encoding kernel (float32).
 __global__ void amplitude_encode_batch_kernel_f32(
     const float* __restrict__ input_batch,
     cuComplex* __restrict__ state_batch,
@@ -317,25 +278,20 @@ __global__ void amplitude_encode_batch_kernel_f32(
     size_t input_len,
     size_t state_len
 ) {
-    // Grid-stride loop pattern for flexibility
     const size_t elements_per_sample = state_len / 2;
     const size_t total_work = num_samples * elements_per_sample;
     const size_t stride = gridDim.x * blockDim.x;
 
     size_t global_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
-    // Process elements in grid-stride fashion
     for (size_t idx = global_idx; idx < total_work; idx += stride) {
-        // Decompose linear index into (sample, element_pair)
         const size_t sample_idx = idx / elements_per_sample;
         const size_t elem_pair = idx % elements_per_sample;
 
-        // Calculate base addresses (strength-reduced)
         const size_t input_base = sample_idx * input_len;
         const size_t state_base = sample_idx * state_len;
         const size_t elem_offset = elem_pair * 2;
 
-        // Load inverse norm (cached by L1)
         const float inv_norm = inv_norms[sample_idx];
 
         float v1, v2;
@@ -357,11 +313,9 @@ __global__ void amplitude_encode_batch_kernel_f32(
             v1 = v2 = 0.0f;
         }
 
-        // Normalize and write as complex numbers
         const cuComplex c1 = make_cuComplex(v1 * inv_norm, 0.0f);
         const cuComplex c2 = make_cuComplex(v2 * inv_norm, 0.0f);
 
-        // Write to global memory (coalesced within warp)
         state_batch[state_base + elem_offset] = c1;
         if (elem_offset + 1 < state_len) {
             state_batch[state_base + elem_offset + 1] = c2;
@@ -397,14 +351,9 @@ int launch_amplitude_encode_batch(
 
     cuDoubleComplex* state_complex_d = static_cast<cuDoubleComplex*>(state_batch_d);
 
-    // Optimal configuration for modern GPUs (SM 7.0+)
-    // - Block size: DEFAULT_BLOCK_SIZE threads (8 warps, good occupancy)
-    // - Grid size: Enough blocks to saturate GPU, but not excessive
     const int blockSize = DEFAULT_BLOCK_SIZE;
     const size_t total_work = num_samples * (state_len / 2);
 
-    // Calculate grid size: aim for high occupancy without too many blocks
-    // Limit to reasonable number of blocks to avoid scheduler overhead
     const size_t blocks_needed = (total_work + blockSize - 1) / blockSize;
     const size_t max_blocks = MAX_GRID_BLOCKS;
     const size_t gridSize = (blocks_needed < max_blocks) ? blocks_needed : max_blocks;
@@ -462,7 +411,6 @@ __global__ void l2_norm_kernel(
     size_t input_len,
     double* __restrict__ out_accum
 ) {
-    // Vectorized double2 loads for bandwidth and coalescing
     const size_t vec_idx = blockIdx.x * blockDim.x + threadIdx.x;
     const size_t stride = gridDim.x * blockDim.x;
 
@@ -497,7 +445,6 @@ __global__ void l2_norm_kernel_f32(
     size_t input_len,
     float* __restrict__ out_accum
 ) {
-    // Vectorized float2 loads for bandwidth and coalescing
     const size_t vec_idx = blockIdx.x * blockDim.x + threadIdx.x;
     const size_t stride = gridDim.x * blockDim.x;
 

qdp/qdp-kernels/src/iqp.cu

@@ -43,8 +43,6 @@ __device__ double compute_phase(
 
     // Single-qubit Z terms: sum_i x_i * data[i]
     for (unsigned int i = 0; i < num_qubits; ++i) {
-        // PR1 Optimization: Use arithmetic multiplication instead of conditional branch
-        // (if x_i != 0) to eliminate warp divergence across threads.
         phase += data[i] * (double)((x >> i) & 1U);
     }
 
@@ -53,7 +51,6 @@ __device__ double compute_phase(
         unsigned int pair_idx = num_qubits;
         for (unsigned int i = 0; i < num_qubits; ++i) {
             for (unsigned int j = i + 1; j < num_qubits; ++j) {
-                // PR1 Optimization: Boolean arithmetic avoids branching for ZZ terms.
                 phase += data[pair_idx] * (double)(((x >> i) & 1U) & ((x >> j) & 1U));
                 pair_idx++;
             }
@@ -89,10 +86,7 @@ __device__ cuDoubleComplex compute_amplitude_naive(
     return make_cuDoubleComplex(real_sum, imag_sum);
 }
 
-// ============================================================================
-// Naive Implementation: O(2^n) per amplitude, O(4^n) for the full state
-// (kept as fallback for small n and verification)
-// ============================================================================
+// Naive O(2^n) per amplitude; fallback when FWT overhead dominates.
 
 __global__ void iqp_encode_kernel_naive(
     const double* __restrict__ data,
@@ -112,9 +106,7 @@ __global__ void iqp_encode_kernel_naive(
 }
 
 
-// ============================================================================
-// FWT O(n * 2^n) Implementation
-// ============================================================================
+// FWT O(n * 2^n) path.
 
 // Step 1: Compute f[x] = exp(i*theta(x)) for all x.
 // Uses a grid-stride loop so large state vectors can reuse a fixed launch size.
@@ -262,9 +254,7 @@ __global__ void normalize_state_kernel(
     }
 }
 
-// ============================================================================
-// Naive O(4^n) Batch Implementation (kept as fallback)
-// ============================================================================
+// Naive batch fallback for small n.
 
 __global__ void iqp_encode_batch_kernel_naive(
     const double* __restrict__ data_batch,
@@ -278,7 +268,6 @@ __global__ void iqp_encode_batch_kernel_naive(
     const size_t total_elements = num_samples * state_len;
     const size_t stride = gridDim.x * blockDim.x;
     const size_t state_mask = state_len - 1;
-    // Normalize by 1/2^n (state_len = 2^n) - hoisted outside the loop
     const double norm = 1.0 / (double)state_len;
 
     for (size_t global_idx = blockIdx.x * blockDim.x + threadIdx.x;
@@ -295,9 +284,7 @@ __global__ void iqp_encode_batch_kernel_naive(
 }
 
 
-// ============================================================================
-// FWT O(n * 2^n) Batch Implementation
-// ============================================================================
+// FWT batch path.
 
 // Step 1: Compute the normalized phase vector for all samples in batch.
 __global__ void iqp_phase_batch_kernel(

qdp/qdp-kernels/src/iqp_tc.cu

@@ -1,3 +1,19 @@
+//
+// Licensed to the Apache Software Foundation (ASF) under one or more
+// contributor license agreements.  See the NOTICE file distributed with
+// this work for additional information regarding copyright ownership.
+// The ASF licenses this file to You under the Apache License, Version 2.0
+// (the "License"); you may not use this file except in compliance with
+// the License.  You may obtain a copy of the License at
+//
+//    http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
 // iqp_tc.cu
 #include <cuda_runtime.h>
 #include <cuComplex.h>
@@ -26,7 +42,7 @@ __device__ double compute_phase_tc(
     return phase;
 }
 
-// PR3: Shared-memory FWT path (Operator Fusion)
+// Shared-memory FWT path (Operator Fusion)
 // Fuses Phase computation, Fast Walsh-Hadamard Transform, and Normalization
 // entirely within Shared Memory. This completely avoids DRAM roundtrips for N <= 12.
 __global__ void iqp_phase_fwt_normalize_tc_kernel(
@@ -162,7 +178,7 @@ void iqp_tc_launch_transpose(const double* d_in, double* d_out, int B, int rows,
     iqp_tc_batch_transpose_kernel<<<grid, block, 0, stream>>>(d_in, d_out, B, rows, cols);
 }
 
-// PR4: Naive Implicit Hadamard GEMM (Fallback before PR5/6 Tensor Core integration)
+// Naive Implicit Hadamard GEMM (Fallback before PR5/6 Tensor Core integration)
 // Computes Y = X * H_K where H_K is a KxK Hadamard matrix generated on-the-fly.
 __global__ void naive_implicit_hadamard_gemm_kernel(const double* __restrict__ X, double* __restrict__ Y, int B, int M, int K, double norm) {
     int k = blockIdx.x * blockDim.x + threadIdx.x;
@@ -185,7 +201,7 @@ void launch_naive_implicit_hadamard(const double* d_in, double* d_out, int B, in
     naive_implicit_hadamard_gemm_kernel<<<grid, block, 0, stream>>>(d_in, d_out, B, M, K, norm);
 }
 
-// PR2: Recombine Real and Imaginary parts back into cuDoubleComplex
+// Recombine Real and Imaginary parts back into cuDoubleComplex
 // This restores the memory layout after Tensor Core matrix multiplications.
 __global__ void recombine_complex_kernel(
     const double* __restrict__ real_part,
@@ -209,7 +225,7 @@ extern "C" int launch_iqp_encode_tc(
     cudaStream_t  stream
 ) {
     if (num_qubits <= FWT_SHARED_MEM_THRESHOLD) {
-        // PR3: For N <= 12, use the fused Shared Memory FWT kernel
+        // For N <= 12, use the fused Shared Memory FWT kernel
         double norm_factor = 1.0 / (double)state_len;
         unsigned int data_len = num_qubits;
         // Request max dynamic shared memory for this kernel
@@ -218,7 +234,7 @@ extern "C" int launch_iqp_encode_tc(
             data_batch_d, static_cast<cuDoubleComplex*>(state_batch_d), num_samples, state_len, num_qubits, data_len, enable_zz, norm_factor
         );
     } else {
-        // PR4: Blocked TC-FWT (Kronecker Product Decomposition)
+        // Blocked TC-FWT (Kronecker Product Decomposition)
         size_t m_samples = num_samples;
         size_t total_elements = m_samples * state_len;
 
@@ -236,18 +252,18 @@ extern "C" int launch_iqp_encode_tc(
         cudaMalloc(&d_out_imag, total_elements * sizeof(double));
         cudaMalloc(&d_temp_real, total_elements * sizeof(double));
         cudaMalloc(&d_temp_imag, total_elements * sizeof(double));
-        
+
         // 1. Initialize Phase (Split Real/Imag)
         unsigned int data_len = num_qubits;
         const size_t blocks = (total_elements + DEFAULT_BLOCK_SIZE - 1) / DEFAULT_BLOCK_SIZE;
         iqp_phase_split_kernel<<<blocks, DEFAULT_BLOCK_SIZE, 0, stream>>>(
             data_batch_d, d_state_real, d_state_imag, num_samples, state_len, num_qubits, data_len, enable_zz
         );
-        
+
         double norm_factor = 1.0 / (double)state_len;
 
         // 3. TC-FWT Step 1: Z = X * H_{n2} (X shape: B*dim1 x dim2)
-        // PR4: Uses Naive GEMM Placeholder. PR5/6 will replace this with Ozaki Implicit Engine.
+        // Uses Naive GEMM Placeholder. PR5/6 will replace this with Ozaki Implicit Engine.
         launch_naive_implicit_hadamard(d_state_real, d_out_real, num_samples * dim1, dim2, dim2, 1.0, stream);
         launch_naive_implicit_hadamard(d_state_imag, d_out_imag, num_samples * dim1, dim2, dim2, 1.0, stream);
 
@@ -262,19 +278,19 @@ extern "C" int launch_iqp_encode_tc(
         // 6. TC-FWT Step 4: Transpose back (B, dim2, dim1) -> (B, dim1, dim2)
         iqp_tc_launch_transpose(d_out_real, d_temp_real, num_samples, dim2, dim1, stream);
         iqp_tc_launch_transpose(d_out_imag, d_temp_imag, num_samples, dim2, dim1, stream);
-        
+
         // 7. Recombine and Write back
         recombine_complex_kernel<<<blocks, DEFAULT_BLOCK_SIZE, 0, stream>>>(
             d_temp_real, d_temp_imag, static_cast<cuDoubleComplex*>(state_batch_d), total_elements
         );
-        
+
         cudaFree(d_state_real);
         cudaFree(d_state_imag);
         cudaFree(d_out_real);
         cudaFree(d_out_imag);
         cudaFree(d_temp_real);
         cudaFree(d_temp_imag);
     }
-    
+
     return (int)cudaSuccess;
 }

qdp/qdp-kernels/src/phase.cu

@@ -49,16 +49,12 @@ __global__ void phase_encode_kernel(
     // φ(idx) = Σ_k phases[k] * b_k,  b_k = (idx >> k) & 1
     double phi = 0.0;
     for (unsigned int bit = 0; bit < num_qubits; ++bit) {
-        // PR1 Optimization: Use cast & multiplication instead of `if ((idx >> bit) & 1U)`
-        // to avoid thread divergence in the GPU warp.
         phi += phases[bit] * (double)((idx >> bit) & 1U);
     }
 
     double re, im;
     sincos(phi, &im, &re);   // re = cos(φ), im = sin(φ)
 
-    // PR1 Optimization: norm_factor is pre-calculated on the host (CPU) and passed 
-    // down to save GPU cycles that would be spent calculating pow(1/sqrt(2), n) repeatedly.
     state[idx] = make_cuDoubleComplex(norm_factor * re, norm_factor * im);
 }
 
@@ -83,14 +79,12 @@ __global__ void phase_encode_batch_kernel(
 
         double phi = 0.0;
         for (unsigned int bit = 0; bit < num_qubits; ++bit) {
-            // PR1 Optimization: Use cast & multiplication to eliminate warp divergence
             phi += phases[bit] * (double)((element_idx >> bit) & 1U);
         }
 
         double re, im;
         sincos(phi, &im, &re);
 
-        // PR1 Optimization: norm_factor is pre-calculated on the host (CPU)
         state_batch[global_idx] = make_cuDoubleComplex(norm_factor * re, norm_factor * im);
     }
 }