[QDP] [feature] Pr5 implicit hadamard engine

qdp/qdp-core/src/gpu/encodings/iqp.rs

@@ -64,6 +64,114 @@ impl IqpEncoder {
             num_qubits
         }
     }
+
+    /// Encode multiple IQP samples using Matrix-Free Implicit Hadamard Tensor Core engine
+    #[cfg(target_os = "linux")]
+    pub fn encode_batch_tc(
+        &self,
+        device: &Arc<CudaDevice>,
+        batch_data: &[f64],
+        num_samples: usize,
+        sample_size: usize,
+        num_qubits: usize,
+    ) -> Result<GpuStateVector> {
+        validate_qubit_count(num_qubits)?;
+        let expected_len = self.expected_data_len(num_qubits);
+
+        if sample_size != expected_len {
+            return Err(MahoutError::InvalidInput(format!(
+                "Batch sample size {} does not match required parameters {}",
+                sample_size, expected_len
+            )));
+        }
+
+        if batch_data.len() != num_samples * sample_size {
+            return Err(MahoutError::InvalidInput(format!(
+                "Batch data length {} does not match num_samples {} * sample_size {}",
+                batch_data.len(),
+                num_samples,
+                sample_size
+            )));
+        }
+
+        for (i, &val) in batch_data.iter().enumerate() {
+            if !val.is_finite() {
+                let sample_idx = i / sample_size;
+                let param_idx = i % sample_size;
+                return Err(MahoutError::InvalidInput(format!(
+                    "Sample {} parameter {} must be finite, got {}",
+                    sample_idx, param_idx, val
+                )));
+            }
+        }
+
+        let state_len = 1 << num_qubits;
+
+        let batch_state_vector = {
+            crate::profile_scope!("GPU::AllocBatchTC");
+            GpuStateVector::new_batch(device, num_samples, num_qubits, Precision::Float64)?
+        };
+
+        let input_bytes = std::mem::size_of_val(batch_data);
+        let data_gpu = {
+            crate::profile_scope!("GPU::H2D_BatchIqpDataTC");
+            device.htod_sync_copy(batch_data).map_err(|e| {
+                map_allocation_error(input_bytes, "IQP TC batch upload", Some(num_qubits), e)
+            })?
+        };
+
+        let state_ptr = batch_state_vector.ptr_f64().ok_or_else(|| {
+            MahoutError::InvalidInput(
+                "Batch state vector precision mismatch (expected float64 buffer)".to_string(),
+            )
+        })?;
+
+        {
+            crate::profile_scope!("GPU::BatchKernelLaunchTC");
+            let ret = unsafe {
+                qdp_kernels::launch_iqp_encode_tc(
+                    *data_gpu.device_ptr() as *const f64,
+                    state_ptr as *mut c_void,
+                    num_samples,
+                    state_len,
+                    num_qubits as u32,
+                    if self.enable_zz { 1 } else { 0 },
+                    std::ptr::null_mut(),
+                )
+            };
+
+            if ret != 0 {
+                return Err(MahoutError::KernelLaunch(format!(
+                    "Batch IQP TC encoding kernel failed: {} ({})",
+                    ret,
+                    cuda_error_to_string(ret)
+                )));
+            }
+        }
+
+        {
+            crate::profile_scope!("GPU::SynchronizeTC");
+            device
+                .synchronize()
+                .map_err(|e| MahoutError::Cuda(format!("Sync failed: {:?}", e)))?;
+        }
+
+        Ok(batch_state_vector)
+    }
+
+    #[cfg(not(target_os = "linux"))]
+    pub fn encode_batch_tc(
+        &self,
+        _device: &Arc<CudaDevice>,
+        _batch_data: &[f64],
+        _num_samples: usize,
+        _sample_size: usize,
+        _num_qubits: usize,
+    ) -> Result<GpuStateVector> {
+        Err(MahoutError::Cuda(
+            "CUDA unavailable (non-Linux stub)".to_string(),
+        ))
+    }
 }
 
 impl QuantumEncoder for IqpEncoder {

qdp/qdp-core/src/lib.rs

@@ -212,6 +212,29 @@ impl QdpEngine {
         self.encode_batch_for_pipeline(batch_data, num_samples, sample_size, num_qubits, encoding)
     }
 
+    /// Encode a batch of IQP samples via the Ozaki Kronecker Tensor Core path.
+    #[cfg(target_os = "linux")]
+    pub fn encode_batch_tc(
+        &self,
+        batch_data: &[f64],
+        num_samples: usize,
+        sample_size: usize,
+        num_qubits: usize,
+    ) -> Result<*mut DLManagedTensor> {
+        crate::profile_scope!("Mahout::EncodeBatchTC");
+
+        let encoder = gpu::encodings::IqpEncoder::full();
+        let state_vector = encoder.encode_batch_tc(
+            &self.device,
+            batch_data,
+            num_samples,
+            sample_size,
+            num_qubits,
+        )?;
+
+        Ok(state_vector.to_dlpack())
+    }
+
     /// Same as [`encode_batch`](Self::encode_batch) with a resolved [`Encoding`] (no string parse).
     pub(crate) fn encode_batch_for_pipeline(
         &self,

qdp/qdp-kernels/build.rs

@@ -26,9 +26,9 @@
 use std::env;
 use std::process::Command;
 
-const DEFAULT_CUBIN_ARCHES: &[&str] = &["75", "80", "86", "89", "90", "100", "120"];
-const DEFAULT_PTX_CANDIDATES: &[&str] = &["120", "100", "90", "89", "86", "80", "75"];
-const LEGACY_FALLBACK_ARCHES: &[&str] = &["75", "80", "86"];
+const DEFAULT_CUBIN_ARCHES: &[&str] = &["80", "86", "89", "90", "100", "120"];
+const DEFAULT_PTX_CANDIDATES: &[&str] = &["120", "100", "90", "89", "86", "80"];
+const LEGACY_FALLBACK_ARCHES: &[&str] = &["80", "86"];
 
 fn add_sm_target(build: &mut cc::Build, arch: &str) {
     build.flag("-gencode");
@@ -164,6 +164,7 @@ fn main() {
     println!("cargo:rerun-if-changed=src/angle.cu");
     println!("cargo:rerun-if-changed=src/validation.cu");
     println!("cargo:rerun-if-changed=src/iqp.cu");
+    println!("cargo:rerun-if-changed=src/iqp_tc.cu");
     println!("cargo:rerun-if-changed=src/phase.cu");
     println!("cargo:rerun-if-env-changed=QDP_NO_CUDA");
     println!("cargo:rerun-if-env-changed=QDP_CUDA_ARCH_LIST");
@@ -224,6 +225,8 @@ fn main() {
         .file("src/angle.cu")
         .file("src/validation.cu")
         .file("src/iqp.cu")
+        .file("src/iqp_tc.cu")
+        .file("src/ImplicitHadamardOzaki.cu")
         .file("src/phase.cu")
         .compile("kernels");
 }

qdp/qdp-kernels/src/AdaptiveOzaki.h

@@ -0,0 +1,86 @@
+#pragma once
+
+#include <cuda_runtime.h>
+#include <vector>
+#include <memory>
+#include <cstdint>
+
+namespace ozaki {
+
+struct PreScanStats {
+    double dynamic_range;
+    double sparsity;
+    int effective_rank;
+};
+
+enum class ExecutionMode {
+    Phase22,
+    Phase23Hetero,
+    Phase24ExtremeMix,
+    Phase25Stacking,
+    Phase26HybridOzaki
+};
+
+struct OzakiConfig {
+    double target_epsilon = 1e-10;
+    size_t vram_limit_mb = 6144;
+    ExecutionMode mode = ExecutionMode::Phase22;
+    bool enable_fp64 = true;
+    bool enable_fp32 = false; // New for Phase 24
+    bool enable_fp16 = false; // New for Phase 24
+    bool enable_residual = true;
+    bool enable_partial_residual = false;
+    bool enable_kernel_fp64 = false;
+    int split_fp64_bits = 8;
+    int split_fp32_bits = 0;  // New for Phase 24
+    int split_fp16_bits = 0;  // New for Phase 24
+    int split_tc_bits = 32;
+    int split_residual_bits = 13;
+    int warp_fp64 = 2;
+    int warp_tc = 4;
+    int warp_residual = 2;
+    double tile_error_guard = 1e-9;
+};
+
+class AdaptiveOzakiEngine {
+public:
+    explicit AdaptiveOzakiEngine(const OzakiConfig& config);
+    ~AdaptiveOzakiEngine();
+
+    void execute(const double* d_A, const double* d_B, double* d_C, int m, int n, int k);
+
+    // Workspace management
+    void allocateWorkspace(int m, int n, int k);
+    void freeWorkspace();
+
+    static void bitShiftToINT8(const double* d_src, int8_t* d_dst_int8, int elem_count, int exponent_shift);
+    static void computeResidualGradedRing(const double* d_A, const double* d_B, const double* d_C, double* d_R, int m, int n, int k);
+    static void accumulateSlicedProduct(const int8_t* d_A_int8, const int8_t* d_B_int8, double* d_C, int m, int n, int k, int total_shift);
+    static void accumulateFusedSlicedProductWMMA(const double* d_A, const double* d_B, double* d_C, int m, int n, int k, int shift_A, int shift_B);
+
+protected:
+    PreScanStats analyzeMatrix(const double* d_A, const double* d_B, int m, int n, int k);
+    int calculateOptimalTile(int m, int n, int k);
+
+private:
+    void profileHardwareRatio();
+
+    OzakiConfig config_;
+    double ratio_fp64_tc = -1.0; // -1 means not profiled yet
+    double ratio_fp32_tc = -1.0;
+    double ratio_fp16_tc = -1.0;
+    double ratio_tf32_tc = -1.0;
+    double ratio_int32_tc = -1.0;
+
+    // Workspace pointers
+    int8_t *dA8_h = nullptr, *dA8_l = nullptr, *dB8_h = nullptr, *dB8_l = nullptr;
+    uint64_t *dmA_h = nullptr, *dmA_l = nullptr, *dmB_h = nullptr, *dmB_l = nullptr;
+    int* d_global_work_queue = nullptr;
+
+    // Low-precision buffers for cross-terms
+    float *dA_hi_f32 = nullptr, *dA_low_f32 = nullptr;
+    float *dB_hi_f32 = nullptr, *dB_low_f32 = nullptr;
+
+    bool workspace_allocated_ = false;
+};
+}