From 8066223d25be51627416fb7a05c5234d8c499e55 Mon Sep 17 00:00:00 2001 From: hachem Date: Wed, 10 Dec 2025 09:26:56 +0100 Subject: [add]: simd --- README.md | 151 ++++++++++-- libpsi-core/Cargo.toml | 3 +- libpsi-core/src/core/kernel.rs | 44 ++++ libpsi-core/src/core/runtime.rs | 44 +++- libpsi-core/src/lib.rs | 1 + libpsi-core/src/maths/mod.rs | 2 + libpsi-core/src/maths/simd.rs | 510 ++++++++++++++++++++++++++++++++++++++++ libpsi-qasm/Cargo.toml | 1 + libpsi-visualizer/Cargo.toml | 2 +- tester/Cargo.toml | 1 + tester/src/main.rs | 8 + tester/src/simd.rs | 210 +++++++++++++++++ 12 files changed, 946 insertions(+), 31 deletions(-) create mode 100644 libpsi-core/src/maths/simd.rs create mode 100644 tester/src/simd.rs diff --git a/README.md b/README.md index 352afdd..91dd69d 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,137 @@ # $\psi$: A Quantum Computational Toolkit -$\psi$ is a powerful quantum computing toolkit designed for simulating quantum circuits on classical computers. - -> :warning: **Warning: Work in Progress** - -## About -- **`libpsi-core`**: A core library for writing and designing quantum circuits, used across all $\psi$ sub-projects. - - **`libpsi-core:runtime`**: A set of four runtimes for executing quantum circuits: - - **`BasicRuntime`**: A single-threaded runtime that executes a quantum circuit statistically, running it $n$ times. - - **`BasicRuntimeMT`**: A multi-threaded version of `BasicRuntime`, enabling parallel execution. - - **`WFEvolution`**: A runtime that applies quantum gates by evolving the quantum wave function over time steps $\Delta t$. - - **`WFEvolutionMT`**: A multi-threaded version of `WFEvolution` for faster execution. - - **`GPUAccelerated`**: A GPU-accelerated runtime for parallel execution of quantum circuits using [NVIDIA CUDA](https://developer.nvidia.com/cuda-toolkit). - - **`libpsi-core:maths`**: A comprehensive mathematics library featuring 32/64-bit complex numbers, vectors (both row and column), matrices, and more. - - **`libpsi-core:core`**: Contains all core quantum components, including quantum gates, classical/quantum bits, and quantum circuits. -- **`libpsi-visualizer`**: An extension of `libpsi-core` that offers visual representations of quantum circuits. It supports both text-based (ASCII) output in the terminal and graphical output using APIs like OpenGL and Vulkan. -- **`libpsi-qasmc`**: An [OpenQASM](https://openqasm.com/) compiler that enables you to write quantum programs, which are then compiled into native executables for classical computers using an [LLVM](https://llvm.org/) backend. + +$\psi$ is a powerful quantum computing toolkit designed for simulating quantum circuits and wavefunction dynamics on classical hardware. + +## Features + +### Quantum Gates + +**Clifford Gates:** +- Single-qubit: H, X, Y, Z, S +- Two-qubit: CNOT, CZ, SWAP +- Three-qubit: CCNOT (Toffoli), CSWAP (Fredkin) + +**Non-Clifford Gates:** +- Fixed: T, $S^\dagger$, $T^\dagger$, $\sqrt{X}$, $\sqrt{X}^\dagger$ +- Parametric rotations: $R_x(\theta)$, $R_y(\theta)$, $R_z(\theta)$, $P(\theta)$ +- General unitaries: $U_1(\lambda)$, $U_2(\phi, \lambda)$, $U_3(\theta, \phi, \lambda)$ +- Controlled parametric: $CR_x(\theta)$, $CR_y(\theta)$, $CR_z(\theta)$, $CP(\theta)$ + +**Custom Gates:** +- Define gates from unitary matrices +- Build composite gates from sequences of operations + +### Runtimes + +| Runtime | Description | +|---------|-------------| +| `BasicRT` | Single-threaded state vector simulation | +| `BasicRTMT` | Multi-threaded parallel simulation | +| `BatchedRT` | Kernel batching with gate fusion optimisation | +| `BatchedRTMT` | Multi-threaded kernel batching | +| `SimdRT` | SIMD-accelerated simulation (AVX2/AVX-512/NEON) | +| `SimdRTMT` | Multi-threaded SIMD acceleration | +| `WFEvolution` | Wave function time evolution (planned) | +| `GPUAccelerated` | CUDA GPU acceleration (planned) | + +### SIMD Acceleration + +Automatic detection and use of platform-specific SIMD instructions: +- **AVX-512**: Modern Intel/AMD processors +- **AVX2+FMA**: Older x86_64 processors +- **NEON**: ARM processors (Apple Silicon, etc.) +- **Scalar fallback**: Universal compatibility + +### Kernel Batching + +Optimisation system that: +- Groups consecutive single-qubit gates on the same qubit +- Fuses gate matrices to reduce operations +- Typically achieves 30–50% kernel reduction and significantly reduces memory bandwidth pressure. + +## Project Structure + +- **`libpsi-core`**: Core quantum simulation library + - `core`: Quantum gates, circuits, registers, and runtimes + - `maths`: Complex numbers, vectors, matrices, SIMD operations +- **`libpsi-visualizer`**: Circuit visualisation (ASCII horizontal/vertical) +- **`tester`**: Comprehensive test suite and benchmarks + +## Quick Start + +```rust +use libpsi_core::{QuantumCircuit, Runtime}; +use std::f64::consts::PI; + +fn main() { + // Create a 3-qubit circuit + let mut circuit = QuantumCircuit::new(3); + + // Build a GHZ state + circuit + .h(0) + .cnot(0, 1) + .cnot(0, 2); + + // Execute with SIMD acceleration + circuit.compute_with(Runtime::SimdRT); + + // Print the quantum state + println!("{}", circuit.state()); +} +``` + +### Parametric Gates + +```rust +let mut circuit = QuantumCircuit::new(2); +circuit + .rx(0, PI / 4.0) // Rotate around X + .ry(0, PI / 3.0) // Rotate around Y + .rz(1, PI / 2.0) // Rotate around Z + .crz(0, 1, PI / 4.0); // Controlled-Rz +``` + +### Custom Gates + +```rust +use libpsi_core::{CustomGateBuilder, complex, matrix}; + +// Build from operations +let bell_gate = CustomGateBuilder::new("BELL", 2) + .h(0) + .cnot(0, 1) + .build(); + +// Or from a matrix +let sqrt_x_matrix = matrix!( + [complex!(0.5, 0.5), complex!(0.5, -0.5)]; + [complex!(0.5, -0.5), complex!(0.5, 0.5)] +); +let sqrt_x = CustomGate::from_matrix("√X", sqrt_x_matrix); +``` + +## Running Tests + +```bash +# Run all tests +cargo run --package tester --release + +# Run specific test modules +cargo run --package tester --release -- clifford +cargo run --package tester --release -- non-clifford +cargo run --package tester --release -- kernels +cargo run --package tester --release -- simd +cargo run --package tester --release -- bench + +# Show help +cargo run --package tester --release -- help +``` ## Disclaimer -This project is a large and ongoing effort, and I try my hardest to deliver the advertised feature, some may not arrive as planned or according to any scheduled timeline. The development process is subject to change based on technical challenges, research priorities and the simultaneous management of multiple on-going projects, spanning both computer science, physics and unrelated domains. + +This project is under active development. Features and APIs may change. Some planned features may not arrive as scheduled due to technical challenges and research priorities. + +## License + +This project is made available under the Apache License, Version 2.0, allowing free use, modification, and distribution with proper attribution. Community contributions, improvements, and research collaborations are encouraged. Full licensing terms can be found in [LICENSE](LICENSE). \ No newline at end of file diff --git a/libpsi-core/Cargo.toml b/libpsi-core/Cargo.toml index 5f17157..01b98a5 100644 --- a/libpsi-core/Cargo.toml +++ b/libpsi-core/Cargo.toml @@ -2,9 +2,10 @@ name = "libpsi-core" version = "0.1.0" edition = "2021" +authors = ["Hachem"] [dependencies] lazy_static = "1.5.0" libm = "0.2.8" -rand = "0.8.5" +rand = "0.9.2" rayon = "1.10" diff --git a/libpsi-core/src/core/kernel.rs b/libpsi-core/src/core/kernel.rs index 3d75425..9f65d09 100644 --- a/libpsi-core/src/core/kernel.rs +++ b/libpsi-core/src/core/kernel.rs @@ -1,3 +1,6 @@ +use crate::maths::simd::{ + apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel, SimdCapability, +}; use crate::{complex, Complex, Matrix}; use rayon::prelude::*; @@ -108,6 +111,47 @@ impl KernelBatch { *state = apply_kernel_parallel(state, kernel, self.num_qubits); } } + + pub fn execute_simd(&self, state: &mut Vec>) { + for kernel in &self.kernels { + if kernel.targets.len() == 1 { + let gate = matrix_to_2x2(&kernel.matrix); + apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); + } else { + *state = apply_kernel(state, kernel, self.num_qubits); + } + } + } + + pub fn execute_simd_parallel(&self, state: &mut Vec>) { + for kernel in &self.kernels { + if kernel.targets.len() == 1 && self.num_qubits >= 10 { + let gate = matrix_to_2x2(&kernel.matrix); + apply_single_qubit_gate_simd_parallel( + state, + &gate, + kernel.targets[0], + self.num_qubits, + ); + } else if kernel.targets.len() == 1 { + let gate = matrix_to_2x2(&kernel.matrix); + apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits); + } else { + *state = apply_kernel_parallel(state, kernel, self.num_qubits); + } + } + } + + pub fn simd_capability(&self) -> SimdCapability { + SimdCapability::detect() + } +} + +fn matrix_to_2x2(matrix: &Matrix>) -> [[Complex; 2]; 2] { + [ + [matrix.data[0], matrix.data[1]], + [matrix.data[2], matrix.data[3]], + ] } fn apply_kernel(state: &[Complex], kernel: &Kernel, num_qubits: usize) -> Vec> { diff --git a/libpsi-core/src/core/runtime.rs b/libpsi-core/src/core/runtime.rs index 69b9d01..a75435e 100644 --- a/libpsi-core/src/core/runtime.rs +++ b/libpsi-core/src/core/runtime.rs @@ -1,9 +1,8 @@ use super::{GateOp, Kernel, KernelBatch, QuantumGate, QuantumRegister, QuantumState}; use crate::gates::{ - cp_matrix, crx_matrix, cry_matrix, crz_matrix, p_matrix, rx_matrix, ry_matrix, - rz_matrix, u1_matrix, u2_matrix, u3_matrix, CNOT, CZ, FREDKIN, HADAMARD, - PAULI_X, PAULI_Y, PAULI_Z, SDG_GATE, SWAP, SXDG_GATE, SX_GATE, S_GATE, TDG_GATE, - TOFFOLI, T_GATE, + cp_matrix, crx_matrix, cry_matrix, crz_matrix, p_matrix, rx_matrix, ry_matrix, rz_matrix, + u1_matrix, u2_matrix, u3_matrix, CNOT, CZ, FREDKIN, HADAMARD, PAULI_X, PAULI_Y, PAULI_Z, + SDG_GATE, SWAP, SXDG_GATE, SX_GATE, S_GATE, TDG_GATE, TOFFOLI, T_GATE, }; use crate::maths::vector::Vector; use crate::{complex, Complex, Matrix}; @@ -18,6 +17,8 @@ pub enum Runtime { BasicRTMT, BatchedRT, BatchedRTMT, + SimdRT, + SimdRTMT, WFEvolution, WFEvolutionMT, GPUAccelerated, @@ -30,6 +31,8 @@ impl Runtime { Runtime::BasicRTMT => Self::compute_basic_mt(num_qubits, operations), Runtime::BatchedRT => Self::compute_batched(num_qubits, operations, false), Runtime::BatchedRTMT => Self::compute_batched(num_qubits, operations, true), + Runtime::SimdRT => Self::compute_simd(num_qubits, operations, false), + Runtime::SimdRTMT => Self::compute_simd(num_qubits, operations, true), Runtime::WFEvolution => { unimplemented!("WFEvolution (Schrödinger equation) runtime not yet implemented") } @@ -111,6 +114,23 @@ impl Runtime { QuantumState::new(state) } + fn compute_simd(num_qubits: usize, operations: &[GateOp], parallel: bool) -> QuantumState { + let dim = 1 << num_qubits; + let mut state: Vec> = vec![complex!(0.0, 0.0); dim]; + state[0] = complex!(1.0, 0.0); + + let mut batch = Self::build_kernel_batch(num_qubits, operations); + batch.optimize(); + + if parallel && num_qubits >= PARALLEL_THRESHOLD { + batch.execute_simd_parallel(&mut state); + } else { + batch.execute_simd(&mut state); + } + + QuantumState::new(state) + } + fn compute_basic(num_qubits: usize, operations: &[GateOp]) -> QuantumState { let names: Vec = (0..num_qubits).map(|i| format!("q{}", i)).collect(); let leaked_names: &'static [String] = Box::leak(names.into_boxed_slice()); @@ -134,14 +154,14 @@ impl Runtime { GateOp::SWAP(a, b) => register.apply_gate(&SWAP, &[*a, *b]), GateOp::CCNOT(c1, c2, t) => register.apply_gate(&TOFFOLI, &[*c1, *c2, *t]), GateOp::CSWAP(c, t1, t2) => register.apply_gate(&FREDKIN, &[*c, *t1, *t2]), - + // Non-Clifford fixed gates GateOp::T(t) => register.apply_gate(&T_GATE, &[*t]), GateOp::Sdg(t) => register.apply_gate(&SDG_GATE, &[*t]), GateOp::Tdg(t) => register.apply_gate(&TDG_GATE, &[*t]), GateOp::Sx(t) => register.apply_gate(&SX_GATE, &[*t]), GateOp::Sxdg(t) => register.apply_gate(&SXDG_GATE, &[*t]), - + // Parametric single-qubit gates (non-Clifford for most angles) GateOp::Rx(t, theta) => { let gate = QuantumGate { @@ -199,7 +219,7 @@ impl Runtime { }; register.apply_gate(&gate, &[*t]); } - + // Controlled parametric gates GateOp::CRx(c, t, theta) => { let gate = QuantumGate { @@ -233,7 +253,7 @@ impl Runtime { }; register.apply_gate(&gate, &[*c, *t]); } - + // Measurement and custom gates GateOp::Measure(_, _) => {} GateOp::Custom(gate, targets) => { @@ -271,14 +291,14 @@ impl Runtime { GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b]), GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t]), GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2]), - + // Non-Clifford fixed gates GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t]), GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t]), GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t]), GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t]), GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t]), - + // Parametric single-qubit gates GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t]), GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t]), @@ -287,13 +307,13 @@ impl Runtime { GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t]), GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t]), GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t]), - + // Controlled parametric gates GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t]), GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t]), GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t]), GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t]), - + // Measurement (skip) and custom gates GateOp::Measure(_, _) => continue, GateOp::Custom(custom_gate, tgts) => { diff --git a/libpsi-core/src/lib.rs b/libpsi-core/src/lib.rs index ff3253c..96c77a0 100644 --- a/libpsi-core/src/lib.rs +++ b/libpsi-core/src/lib.rs @@ -5,6 +5,7 @@ pub use maths::complex::*; pub use maths::format::*; pub use maths::matrix::*; pub use maths::numeric::*; +pub use maths::simd::*; pub use maths::vector::*; pub use core::circuit::*; diff --git a/libpsi-core/src/maths/mod.rs b/libpsi-core/src/maths/mod.rs index e716181..85f4872 100644 --- a/libpsi-core/src/maths/mod.rs +++ b/libpsi-core/src/maths/mod.rs @@ -2,6 +2,7 @@ pub mod complex; pub mod format; pub mod matrix; pub mod numeric; +pub mod simd; pub mod vector; pub mod vector_ops; @@ -9,4 +10,5 @@ pub use complex::*; pub use format::*; pub use matrix::*; pub use numeric::*; +pub use simd::*; pub use vector::*; diff --git a/libpsi-core/src/maths/simd.rs b/libpsi-core/src/maths/simd.rs new file mode 100644 index 0000000..de0370c --- /dev/null +++ b/libpsi-core/src/maths/simd.rs @@ -0,0 +1,510 @@ +use crate::{complex, Complex}; + +#[cfg(target_arch = "x86_64")] +use std::arch::x86_64::*; + +#[cfg(target_arch = "aarch64")] +use std::arch::aarch64::*; + +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum SimdCapability { + None, + #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] + Avx2, + #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] + Avx512, + #[cfg(target_arch = "aarch64")] + Neon, +} + +impl SimdCapability { + pub fn detect() -> Self { + #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] + { + if is_x86_feature_detected!("avx512f") && is_x86_feature_detected!("avx512dq") { + return SimdCapability::Avx512; + } + if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") { + return SimdCapability::Avx2; + } + } + + #[cfg(target_arch = "aarch64")] + { + return SimdCapability::Neon; + } + + #[allow(unreachable_code)] + SimdCapability::None + } + + pub fn name(&self) -> &'static str { + match self { + SimdCapability::None => "Scalar", + #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] + SimdCapability::Avx2 => "AVX2+FMA", + #[cfg(any(target_arch = "x86_64", target_arch = "x86"))] + SimdCapability::Avx512 => "AVX-512", + #[cfg(target_arch = "aarch64")] + SimdCapability::Neon => "NEON", + } + } +} + +pub fn apply_single_qubit_gate_simd( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + let capability = SimdCapability::detect(); + + match capability { + #[cfg(target_arch = "x86_64")] + SimdCapability::Avx2 => unsafe { + apply_single_qubit_avx2(state, gate, target, num_qubits); + }, + #[cfg(target_arch = "x86_64")] + SimdCapability::Avx512 => unsafe { + apply_single_qubit_avx512(state, gate, target, num_qubits); + }, + #[cfg(target_arch = "aarch64")] + SimdCapability::Neon => unsafe { + apply_single_qubit_neon(state, gate, target, num_qubits); + }, + _ => { + apply_single_qubit_scalar(state, gate, target, num_qubits); + } + } +} + +#[cfg(target_arch = "x86_64")] +#[target_feature(enable = "avx2", enable = "fma")] +unsafe fn apply_single_qubit_avx2( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + let target_bit = num_qubits - 1 - target; + let step = 1 << target_bit; + let dim = 1 << num_qubits; + + let g00 = gate[0][0]; + let g01 = gate[0][1]; + let g10 = gate[1][0]; + let g11 = gate[1][1]; + + let pairs: Vec<(usize, usize)> = (0..dim) + .filter(|&i| (i >> target_bit) & 1 == 0) + .map(|i| (i, i | step)) + .collect(); + + let chunks = pairs.len() / 2; + + for chunk_idx in 0..chunks { + let (i0, j0) = pairs[chunk_idx * 2]; + let (i1, j1) = pairs[chunk_idx * 2 + 1]; + + let s0_re = _mm256_set_pd( + state[j1].real, + state[i1].real, + state[j0].real, + state[i0].real, + ); + let s0_im = _mm256_set_pd( + state[j1].imaginary, + state[i1].imaginary, + state[j0].imaginary, + state[i0].imaginary, + ); + + let g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real); + let g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary); + let g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real); + let g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary); + + let prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0)); + let prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0)); + + let prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1)); + let prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1)); + + let mut res0_re = [0.0f64; 4]; + let mut res0_im = [0.0f64; 4]; + let mut res1_re = [0.0f64; 4]; + let mut res1_im = [0.0f64; 4]; + + _mm256_storeu_pd(res0_re.as_mut_ptr(), prod0_re); + _mm256_storeu_pd(res0_im.as_mut_ptr(), prod0_im); + _mm256_storeu_pd(res1_re.as_mut_ptr(), prod1_re); + _mm256_storeu_pd(res1_im.as_mut_ptr(), prod1_im); + + state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); + state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); + state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); + state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); + } + + for &(i, j) in pairs.iter().skip(chunks * 2) { + let s0 = state[i]; + let s1 = state[j]; + + let new0 = complex!( + s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real + - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + + s0.imaginary * g00.real + + s1.real * g01.imaginary + + s1.imaginary * g01.real + ); + + let new1 = complex!( + s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real + - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + + s0.imaginary * g10.real + + s1.real * g11.imaginary + + s1.imaginary * g11.real + ); + + state[i] = new0; + state[j] = new1; + } +} + +#[cfg(target_arch = "x86_64")] +#[target_feature(enable = "avx512f", enable = "avx512dq")] +unsafe fn apply_single_qubit_avx512( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + let target_bit = num_qubits - 1 - target; + let step = 1 << target_bit; + let dim = 1 << num_qubits; + + let g00 = gate[0][0]; + let g01 = gate[0][1]; + let g10 = gate[1][0]; + let g11 = gate[1][1]; + + let pairs: Vec<(usize, usize)> = (0..dim) + .filter(|&i| (i >> target_bit) & 1 == 0) + .map(|i| (i, i | step)) + .collect(); + + let chunks = pairs.len() / 4; + + for chunk_idx in 0..chunks { + let base = chunk_idx * 4; + let (i0, j0) = pairs[base]; + let (i1, j1) = pairs[base + 1]; + let (i2, j2) = pairs[base + 2]; + let (i3, j3) = pairs[base + 3]; + + let s0_re = _mm512_set_pd( + state[j3].real, + state[i3].real, + state[j2].real, + state[i2].real, + state[j1].real, + state[i1].real, + state[j0].real, + state[i0].real, + ); + let s0_im = _mm512_set_pd( + state[j3].imaginary, + state[i3].imaginary, + state[j2].imaginary, + state[i2].imaginary, + state[j1].imaginary, + state[i1].imaginary, + state[j0].imaginary, + state[i0].imaginary, + ); + + let g_re_0 = _mm512_set_pd( + g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, + ); + let g_im_0 = _mm512_set_pd( + g01.imaginary, + g00.imaginary, + g01.imaginary, + g00.imaginary, + g01.imaginary, + g00.imaginary, + g01.imaginary, + g00.imaginary, + ); + let g_re_1 = _mm512_set_pd( + g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, + ); + let g_im_1 = _mm512_set_pd( + g11.imaginary, + g10.imaginary, + g11.imaginary, + g10.imaginary, + g11.imaginary, + g10.imaginary, + g11.imaginary, + g10.imaginary, + ); + + let prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0)); + let prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0)); + let prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1)); + let prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1)); + + let mut res0_re = [0.0f64; 8]; + let mut res0_im = [0.0f64; 8]; + let mut res1_re = [0.0f64; 8]; + let mut res1_im = [0.0f64; 8]; + + _mm512_storeu_pd(res0_re.as_mut_ptr(), prod0_re); + _mm512_storeu_pd(res0_im.as_mut_ptr(), prod0_im); + _mm512_storeu_pd(res1_re.as_mut_ptr(), prod1_re); + _mm512_storeu_pd(res1_im.as_mut_ptr(), prod1_im); + + state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]); + state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]); + state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]); + state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]); + state[i2] = complex!(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]); + state[j2] = complex!(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]); + state[i3] = complex!(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]); + state[j3] = complex!(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]); + } + + for &(i, j) in pairs.iter().skip(chunks * 4) { + let s0 = state[i]; + let s1 = state[j]; + + let new0 = complex!( + s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real + - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + + s0.imaginary * g00.real + + s1.real * g01.imaginary + + s1.imaginary * g01.real + ); + + let new1 = complex!( + s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real + - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + + s0.imaginary * g10.real + + s1.real * g11.imaginary + + s1.imaginary * g11.real + ); + + state[i] = new0; + state[j] = new1; + } +} + +#[cfg(target_arch = "aarch64")] +unsafe fn apply_single_qubit_neon( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + let target_bit = num_qubits - 1 - target; + let step = 1 << target_bit; + let dim = 1 << num_qubits; + + let g00 = gate[0][0]; + let g01 = gate[0][1]; + let g10 = gate[1][0]; + let g11 = gate[1][1]; + + let pairs: Vec<(usize, usize)> = (0..dim) + .filter(|&i| (i >> target_bit) & 1 == 0) + .map(|i| (i, i | step)) + .collect(); + + let chunks = pairs.len() / 2; + + for chunk_idx in 0..chunks { + let (i0, j0) = pairs[chunk_idx * 2]; + let (i1, j1) = pairs[chunk_idx * 2 + 1]; + + let s0_0 = state[i0]; + let s1_0 = state[j0]; + let s0_1 = state[i1]; + let s1_1 = state[j1]; + + let s0_re = vld1q_f64([s0_0.real, s0_1.real].as_ptr()); + let s0_im = vld1q_f64([s0_0.imaginary, s0_1.imaginary].as_ptr()); + let s1_re = vld1q_f64([s1_0.real, s1_1.real].as_ptr()); + let s1_im = vld1q_f64([s1_0.imaginary, s1_1.imaginary].as_ptr()); + + let g00_re = vdupq_n_f64(g00.real); + let g00_im = vdupq_n_f64(g00.imaginary); + let g01_re = vdupq_n_f64(g01.real); + let g01_im = vdupq_n_f64(g01.imaginary); + let g10_re = vdupq_n_f64(g10.real); + let g10_im = vdupq_n_f64(g10.imaginary); + let g11_re = vdupq_n_f64(g11.real); + let g11_im = vdupq_n_f64(g11.imaginary); + + let new0_re = vaddq_f64( + vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im), + vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im), + ); + let new0_im = vaddq_f64( + vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re), + vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re), + ); + + let new1_re = vaddq_f64( + vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im), + vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im), + ); + let new1_im = vaddq_f64( + vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re), + vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re), + ); + + state[i0] = complex!(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0)); + state[j0] = complex!(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0)); + state[i1] = complex!(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1)); + state[j1] = complex!(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1)); + } + + for &(i, j) in pairs.iter().skip(chunks * 2) { + let s0 = state[i]; + let s1 = state[j]; + + let new0 = complex!( + s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real + - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + + s0.imaginary * g00.real + + s1.real * g01.imaginary + + s1.imaginary * g01.real + ); + + let new1 = complex!( + s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real + - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + + s0.imaginary * g10.real + + s1.real * g11.imaginary + + s1.imaginary * g11.real + ); + + state[i] = new0; + state[j] = new1; + } +} + +fn apply_single_qubit_scalar( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + let target_bit = num_qubits - 1 - target; + let step = 1 << target_bit; + let dim = 1 << num_qubits; + + let g00 = gate[0][0]; + let g01 = gate[0][1]; + let g10 = gate[1][0]; + let g11 = gate[1][1]; + + for i in 0..dim { + if (i >> target_bit) & 1 == 1 { + continue; + } + + let j = i | step; + let s0 = state[i]; + let s1 = state[j]; + + let new0 = complex!( + s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real + - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + + s0.imaginary * g00.real + + s1.real * g01.imaginary + + s1.imaginary * g01.real + ); + + let new1 = complex!( + s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real + - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + + s0.imaginary * g10.real + + s1.real * g11.imaginary + + s1.imaginary * g11.real + ); + + state[i] = new0; + state[j] = new1; + } +} + +pub fn apply_single_qubit_gate_simd_parallel( + state: &mut [Complex], + gate: &[[Complex; 2]; 2], + target: usize, + num_qubits: usize, +) { + use rayon::prelude::*; + + let target_bit = num_qubits - 1 - target; + let step = 1 << target_bit; + let dim = 1 << num_qubits; + + let g00 = gate[0][0]; + let g01 = gate[0][1]; + let g10 = gate[1][0]; + let g11 = gate[1][1]; + + let pairs: Vec<(usize, usize)> = (0..dim) + .filter(|&i| (i >> target_bit) & 1 == 0) + .map(|i| (i, i | step)) + .collect(); + + let results: Vec<(usize, usize, Complex, Complex)> = pairs + .par_iter() + .map(|&(i, j)| { + let s0 = state[i]; + let s1 = state[j]; + + let new0 = complex!( + s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real + - s1.imaginary * g01.imaginary, + s0.real * g00.imaginary + + s0.imaginary * g00.real + + s1.real * g01.imaginary + + s1.imaginary * g01.real + ); + + let new1 = complex!( + s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real + - s1.imaginary * g11.imaginary, + s0.real * g10.imaginary + + s0.imaginary * g10.real + + s1.real * g11.imaginary + + s1.imaginary * g11.real + ); + + (i, j, new0, new1) + }) + .collect(); + + for (i, j, new0, new1) in results { + state[i] = new0; + state[j] = new1; + } +} + +pub fn get_simd_info() -> String { + let cap = SimdCapability::detect(); + format!("SIMD: {}", cap.name()) +} diff --git a/libpsi-qasm/Cargo.toml b/libpsi-qasm/Cargo.toml index e36f63f..f0ecc4a 100644 --- a/libpsi-qasm/Cargo.toml +++ b/libpsi-qasm/Cargo.toml @@ -2,5 +2,6 @@ name = "libpsi-qasm" version = "0.1.0" edition = "2021" +authors = ["Hachem"] [dependencies] diff --git a/libpsi-visualizer/Cargo.toml b/libpsi-visualizer/Cargo.toml index aae9e45..c1c9944 100644 --- a/libpsi-visualizer/Cargo.toml +++ b/libpsi-visualizer/Cargo.toml @@ -2,7 +2,7 @@ name = "libpsi-visualizer" version = "0.1.0" edition = "2021" +authors = ["Hachem"] [dependencies] libpsi-core = { path = "../libpsi-core" } -term_size = "0.3.2" diff --git a/tester/Cargo.toml b/tester/Cargo.toml index c536f89..e912889 100644 --- a/tester/Cargo.toml +++ b/tester/Cargo.toml @@ -2,6 +2,7 @@ name = "tester" version = "0.1.0" edition = "2021" +authors = ["Hachem"] [dependencies] libpsi-core ={ path = "../libpsi-core"} diff --git a/tester/src/main.rs b/tester/src/main.rs index 1fb4d32..854472c 100644 --- a/tester/src/main.rs +++ b/tester/src/main.rs @@ -4,6 +4,7 @@ mod common; mod custom_gates; mod kernels; mod non_clifford; +mod simd; use common::{print_benchmark_table, print_summary, BenchmarkResult}; use std::env; @@ -23,6 +24,7 @@ fn print_usage() { println!(" non-clifford Run non-Clifford gate tests only"); println!(" custom Run custom gate tests only"); println!(" kernels Run kernel batching tests only"); + println!(" simd Run SIMD acceleration tests only"); println!(" bench Run benchmark tests only"); println!(" help Show this help message"); println!(); @@ -31,6 +33,7 @@ fn print_usage() { println!(" tester clifford # Run only Clifford gate tests"); println!(" tester non-clifford # Run only rotation/parametric gate tests"); println!(" tester kernels # Run only kernel batching tests"); + println!(" tester simd # Run only SIMD tests"); println!(" tester custom bench # Run custom gates and benchmarks"); } @@ -54,6 +57,7 @@ fn main() { let run_non_clifford = run_all || args.iter().any(|a| a == "non-clifford"); let run_custom = run_all || args.iter().any(|a| a == "custom"); let run_kernels = run_all || args.iter().any(|a| a == "kernels"); + let run_simd = run_all || args.iter().any(|a| a == "simd"); let run_bench = run_all || args.iter().any(|a| a == "bench"); if run_clifford { @@ -72,6 +76,10 @@ fn main() { kernels::run_all(&mut results); } + if run_simd { + simd::run_all(&mut results); + } + if run_bench { benchmarks::run_all(&mut results); } diff --git a/tester/src/simd.rs b/tester/src/simd.rs new file mode 100644 index 0000000..7ab71b6 --- /dev/null +++ b/tester/src/simd.rs @@ -0,0 +1,210 @@ +use crate::common::{print_section, states_equal, BenchmarkResult}; +use libpsi_core::{get_simd_info, QuantumCircuit, Runtime}; +use std::f64::consts::PI; +use std::time::Instant; + +pub fn run_all(results: &mut Vec) { + println!("═══════════════════════════════════════════════════════════════"); + println!(" SIMD ACCELERATION TESTS"); + println!("═══════════════════════════════════════════════════════════════\n"); + + println!("Detected: {}\n", get_simd_info()); + + test_simd_correctness(results); + test_simd_vs_batched(results); + test_simd_large_circuits(results); +} + +pub fn test_simd_correctness(results: &mut Vec) { + print_section("SIMD Correctness Verification"); + + let test_cases: Vec<(&str, Box QuantumCircuit>)> = vec![ + ( + "Bell State", + Box::new(|| { + let mut c = QuantumCircuit::new(2); + c.h(0).cnot(0, 1); + c + }), + ), + ( + "GHZ-3", + Box::new(|| { + let mut c = QuantumCircuit::new(3); + c.h(0).cnot(0, 1).cnot(0, 2); + c + }), + ), + ( + "Rotation Chain", + Box::new(|| { + let mut c = QuantumCircuit::new(3); + c.rx(0, PI / 4.0) + .ry(0, PI / 4.0) + .rz(0, PI / 4.0) + .rx(1, PI / 3.0) + .ry(1, PI / 3.0); + c + }), + ), + ( + "Mixed Single-Qubit", + Box::new(|| { + let mut c = QuantumCircuit::new(4); + c.h(0).t(0).s(0).x(0).h(1).y(1).z(1).h(2).t(2).h(3).s(3); + c + }), + ), + ]; + + for (name, builder) in test_cases { + let mut basic = builder(); + basic.compute_with(Runtime::BasicRT); + + let mut simd = builder(); + simd.compute_with(Runtime::SimdRT); + + let match_result = states_equal(basic.state(), simd.state()); + + println!( + "{}: {}", + name, + if match_result { + "✓ Match" + } else { + "✗ MISMATCH" + } + ); + + results.push(BenchmarkResult { + name: format!("SIMD verify: {}", name), + basic_time: std::time::Duration::from_micros(0), + mt_time: std::time::Duration::from_micros(0), + results_match: match_result, + }); + } + println!(); +} + +pub fn test_simd_vs_batched(results: &mut Vec) { + print_section("SIMD vs Batched Runtime Comparison"); + + let test_cases: Vec<(&str, Box QuantumCircuit>)> = vec![ + ( + "Single-Qubit Heavy (6q)", + Box::new(|| { + let mut c = QuantumCircuit::new(6); + for q in 0..6 { + c.h(q).t(q).s(q).x(q).y(q).z(q); + } + c + }), + ), + ( + "Rotation Circuit (5q)", + Box::new(|| { + let mut c = QuantumCircuit::new(5); + for q in 0..5 { + c.rx(q, PI / 4.0).ry(q, PI / 3.0).rz(q, PI / 6.0); + } + c + }), + ), + ( + "Deep Single-Qubit (4q)", + Box::new(|| { + let mut c = QuantumCircuit::new(4); + for _ in 0..10 { + for q in 0..4 { + c.h(q).t(q); + } + } + c + }), + ), + ]; + + for (name, builder) in test_cases { + let mut batched = builder(); + let start = Instant::now(); + batched.compute_with(Runtime::BatchedRT); + let batched_time = start.elapsed(); + + let mut simd = builder(); + let start = Instant::now(); + simd.compute_with(Runtime::SimdRT); + let simd_time = start.elapsed(); + + let match_result = states_equal(batched.state(), simd.state()); + + let speedup = batched_time.as_secs_f64() / simd_time.as_secs_f64(); + println!( + "{}: Batched={:.2}μs, SIMD={:.2}μs, Speedup={:.2}x, Match={}", + name, + batched_time.as_secs_f64() * 1_000_000.0, + simd_time.as_secs_f64() * 1_000_000.0, + speedup, + if match_result { "✓" } else { "✗" } + ); + + results.push(BenchmarkResult { + name: format!("SIMD: {}", name), + basic_time: batched_time, + mt_time: simd_time, + results_match: match_result, + }); + } + println!(); +} + +pub fn test_simd_large_circuits(results: &mut Vec) { + print_section("SIMD on Large Circuits (Multi-threaded)"); + + let sizes = [8, 10, 12]; + + for &n in &sizes { + let builder = || { + let mut circuit = QuantumCircuit::new(n); + for i in 0..n { + circuit.h(i); + } + for i in 0..(n - 1) { + circuit.cnot(i, i + 1); + } + for i in 0..n { + circuit.t(i).s(i); + } + circuit + }; + + let mut batched_mt = builder(); + let start = Instant::now(); + batched_mt.compute_with(Runtime::BatchedRTMT); + let batched_time = start.elapsed(); + + let mut simd_mt = builder(); + let start = Instant::now(); + simd_mt.compute_with(Runtime::SimdRTMT); + let simd_time = start.elapsed(); + + let match_result = states_equal(batched_mt.state(), simd_mt.state()); + + let speedup = batched_time.as_secs_f64() / simd_time.as_secs_f64(); + println!( + "{}-qubit: BatchedMT={:.3}ms, SIMD_MT={:.3}ms, Speedup={:.2}x, Match={}", + n, + batched_time.as_secs_f64() * 1000.0, + simd_time.as_secs_f64() * 1000.0, + speedup, + if match_result { "✓" } else { "✗" } + ); + + results.push(BenchmarkResult { + name: format!("{}-qubit SIMD", n), + basic_time: batched_time, + mt_time: simd_time, + results_match: match_result, + }); + } + println!(); +} -- cgit v1.3