aboutsummaryrefslogtreecommitdiff
path: root/src/core
diff options
context:
space:
mode:
Diffstat (limited to 'src/core')
-rw-r--r--src/core/kernel.c90
-rw-r--r--src/core/kernel.rs621
-rw-r--r--src/core/runtime.c2
-rw-r--r--src/core/runtime.rs585
4 files changed, 67 insertions, 1231 deletions
diff --git a/src/core/kernel.c b/src/core/kernel.c
index 377982f..063d47f 100644
--- a/src/core/kernel.c
+++ b/src/core/kernel.c
@@ -6,6 +6,8 @@
#include <stdlib.h>
#include <string.h>
+#include "parallel.h"
+
static char* dup_string(const char* s)
{
size_t n = strlen(s) + 1;
@@ -159,8 +161,49 @@ bool psi_fuse_kernels(struct PsiKernel a, struct PsiKernel b, struct PsiKernel*
return true;
}
-static struct PsiComplex* apply_kernel(const struct PsiComplex* state, struct PsiKernel kernel,
- size_t num_qubits)
+struct PsiKernelApply
+{
+ const struct PsiComplex* state;
+ struct PsiComplex* new_state;
+ struct PsiMatrix matrix;
+ const size_t* target_bits;
+ size_t g;
+ size_t gate_dim;
+ size_t non_target_mask;
+};
+
+static void kernel_apply_range(size_t start, size_t end, void* vctx)
+{
+ struct PsiKernelApply* c = vctx;
+
+ for (size_t i = start; i < end; i++)
+ {
+ size_t target_idx = 0;
+ for (size_t k = 0; k < c->g; k++)
+ if ((i >> c->target_bits[k]) & 1)
+ target_idx |= (size_t)1 << (c->g - 1 - k);
+
+ struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+ for (size_t j = 0; j < c->gate_dim; j++)
+ {
+ struct PsiComplex gate_elem = c->matrix.data[target_idx * c->gate_dim + j];
+ if (fabs(gate_elem.real) < 1e-15 && fabs(gate_elem.imaginary) < 1e-15)
+ continue;
+
+ size_t source_idx = i & c->non_target_mask;
+ for (size_t k = 0; k < c->g; k++)
+ if ((j >> (c->g - 1 - k)) & 1)
+ source_idx |= (size_t)1 << c->target_bits[k];
+
+ sum = psi_add_complex(sum, psi_mul_complex(gate_elem, c->state[source_idx]));
+ }
+
+ c->new_state[i] = sum;
+ }
+}
+
+static struct PsiComplex* run_kernel(const struct PsiComplex* state, struct PsiKernel kernel,
+ size_t num_qubits, bool parallel)
{
size_t dim = (size_t)1 << num_qubits;
size_t g = kernel.target_count;
@@ -178,35 +221,25 @@ static struct PsiComplex* apply_kernel(const struct PsiComplex* state, struct Ps
struct PsiComplex* new_state = malloc(dim * sizeof(struct PsiComplex));
assert(new_state != NULL);
- for (size_t i = 0; i < dim; i++)
- {
- size_t target_idx = 0;
- for (size_t k = 0; k < g; k++)
- if ((i >> target_bits[k]) & 1)
- target_idx |= (size_t)1 << (g - 1 - k);
+ struct PsiKernelApply ctx = {
+ state, new_state, kernel.matrix, target_bits, g, gate_dim, non_target_mask,
+ };
- struct PsiComplex sum = psi_new_complex(0.0, 0.0);
- for (size_t j = 0; j < gate_dim; j++)
- {
- struct PsiComplex gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
- if (fabs(gate_elem.real) < 1e-15 && fabs(gate_elem.imaginary) < 1e-15)
- continue;
-
- size_t source_idx = i & non_target_mask;
- for (size_t k = 0; k < g; k++)
- if ((j >> (g - 1 - k)) & 1)
- source_idx |= (size_t)1 << target_bits[k];
-
- sum = psi_add_complex(sum, psi_mul_complex(gate_elem, state[source_idx]));
- }
-
- new_state[i] = sum;
- }
+ if (parallel)
+ psi_parallel_for(dim, kernel_apply_range, &ctx);
+ else
+ kernel_apply_range(0, dim, &ctx);
free(target_bits);
return new_state;
}
+static struct PsiComplex* apply_kernel(const struct PsiComplex* state, struct PsiKernel kernel,
+ size_t num_qubits)
+{
+ return run_kernel(state, kernel, num_qubits, false);
+}
+
struct PsiKernelBatch psi_new_kernel_batch(size_t num_qubits)
{
struct PsiKernelBatch batch;
@@ -298,6 +331,13 @@ void psi_apply_kernel(struct PsiVector* state, struct PsiKernel kernel, size_t n
state->data = next;
}
+void psi_apply_kernel_parallel(struct PsiVector* state, struct PsiKernel kernel, size_t num_qubits)
+{
+ struct PsiComplex* next = run_kernel(state->data, kernel, num_qubits, true);
+ free(state->data);
+ state->data = next;
+}
+
static void push_kernel(struct PsiKernel** kernels, size_t* count, size_t* capacity,
struct PsiKernel kernel)
{
diff --git a/src/core/kernel.rs b/src/core/kernel.rs
deleted file mode 100644
index 7f73458..0000000
--- a/src/core/kernel.rs
+++ /dev/null
@@ -1,621 +0,0 @@
-use crate::maths::simd::{
- apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel, SimdCapability,
-};
-use crate::{complex, Complex, Matrix};
-use rayon::prelude::*;
-use std::collections::HashSet;
-
-#[derive(Clone, Copy, Debug, PartialEq, Eq)]
-pub enum GateType {
- Diagonal,
- NonDiagonal,
- Controlled,
-}
-
-#[derive(Clone)]
-pub struct Kernel {
- pub matrix: Matrix<Complex<f64>>,
- pub targets: Vec<usize>,
- pub name: String,
- pub gate_type: GateType,
-}
-
-impl Kernel {
- pub fn new(name: &str, matrix: Matrix<Complex<f64>>, targets: Vec<usize>) -> Self {
- let gate_type = Self::detect_gate_type(name, &matrix);
- Self {
- matrix,
- targets,
- name: name.to_string(),
- gate_type,
- }
- }
-
- fn detect_gate_type(name: &str, matrix: &Matrix<Complex<f64>>) -> GateType {
- let diagonal_gates = [
- "Z", "S", "T", "Sdg", "Tdg", "Rz", "P", "U1", "CZ", "CP", "CRz",
- ];
- if diagonal_gates.iter().any(|&g| name.starts_with(g)) {
- return GateType::Diagonal;
- }
-
- let controlled_gates = [
- "CNOT", "CZ", "SWAP", "CRx", "CRy", "CRz", "CP", "CCNOT", "CSWAP",
- ];
- if controlled_gates.iter().any(|&g| name.starts_with(g)) {
- return GateType::Controlled;
- }
-
- if matrix.rows == 2 && matrix.cols == 2 {
- let is_diag = matrix.data[1].real.abs() < 1e-10
- && matrix.data[1].imaginary.abs() < 1e-10
- && matrix.data[2].real.abs() < 1e-10
- && matrix.data[2].imaginary.abs() < 1e-10;
- if is_diag {
- return GateType::Diagonal;
- }
- }
-
- GateType::NonDiagonal
- }
-
- pub fn num_qubits(&self) -> usize {
- self.targets.len()
- }
-
- pub fn target_set(&self) -> HashSet<usize> {
- self.targets.iter().cloned().collect()
- }
-
- pub fn shares_qubits(&self, other: &Kernel) -> bool {
- self.targets.iter().any(|t| other.targets.contains(t))
- }
-
- pub fn commutes_with(&self, other: &Kernel) -> bool {
- if !self.shares_qubits(other) {
- return true;
- }
-
- if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal
- && self.targets == other.targets {
- return true;
- }
-
- false
- }
-
- pub fn can_fuse_with(&self, other: &Kernel) -> bool {
- if self.targets.len() != 1 || other.targets.len() != 1 {
- return false;
- }
- self.targets[0] == other.targets[0]
- }
-
- pub fn fuse(&self, other: &Kernel) -> Option<Kernel> {
- if !self.can_fuse_with(other) {
- return None;
- }
- let fused_matrix = other.matrix.dot(&self.matrix)?;
- let new_type =
- if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal {
- GateType::Diagonal
- } else {
- GateType::NonDiagonal
- };
- Some(Kernel {
- matrix: fused_matrix,
- targets: self.targets.clone(),
- name: format!("{}+{}", self.name, other.name),
- gate_type: new_type,
- })
- }
-}
-
-pub struct KernelBatch {
- kernels: Vec<Kernel>,
- num_qubits: usize,
-}
-
-impl KernelBatch {
- pub fn new(num_qubits: usize) -> Self {
- Self {
- kernels: Vec::new(),
- num_qubits,
- }
- }
-
- pub fn add(&mut self, kernel: Kernel) {
- self.kernels.push(kernel);
- }
-
- pub fn len(&self) -> usize {
- self.kernels.len()
- }
-
- pub fn is_empty(&self) -> bool {
- self.kernels.is_empty()
- }
-
- pub fn kernels(&self) -> &[Kernel] {
- &self.kernels
- }
-
- pub fn optimize(&mut self) {
- if self.kernels.len() < 2 {
- return;
- }
-
- let mut optimized: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
- let mut i = 0;
-
- while i < self.kernels.len() {
- let current = &self.kernels[i];
-
- if i + 1 < self.kernels.len() {
- let next = &self.kernels[i + 1];
- if let Some(fused) = current.fuse(next) {
- optimized.push(fused);
- i += 2;
- continue;
- }
- }
-
- optimized.push(current.clone());
- i += 1;
- }
-
- self.kernels = optimized;
- }
-
- pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
- for kernel in &self.kernels {
- *state = apply_kernel(state, kernel, self.num_qubits);
- }
- }
-
- pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
- for kernel in &self.kernels {
- *state = apply_kernel_parallel(state, kernel, self.num_qubits);
- }
- }
-
- pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
- 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<Complex<f64>>) {
- 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<f64>>) -> [[Complex<f64>; 2]; 2] {
- [
- [matrix.data[0], matrix.data[1]],
- [matrix.data[2], matrix.data[3]],
- ]
-}
-
-fn apply_kernel(state: &[Complex<f64>], kernel: &Kernel, num_qubits: usize) -> Vec<Complex<f64>> {
- let dim = 1 << num_qubits;
- let g = kernel.targets.len();
- let gate_dim = 1 << g;
-
- let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
- let mut non_target_mask: usize = (1 << num_qubits) - 1;
- for &pos in &target_bits {
- non_target_mask &= !(1 << pos);
- }
-
- let mut new_state = vec![complex!(0.0, 0.0); dim];
-
- for (i, new_val) in new_state.iter_mut().enumerate() {
- let mut target_idx = 0usize;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (i >> pos) & 1 == 1 {
- target_idx |= 1 << (g - 1 - k);
- }
- }
-
- let mut sum = complex!(0.0, 0.0);
-
- for j in 0..gate_dim {
- let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
- if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
- continue;
- }
-
- let mut source_idx = i & non_target_mask;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (j >> (g - 1 - k)) & 1 == 1 {
- source_idx |= 1 << pos;
- }
- }
-
- sum += gate_elem * state[source_idx];
- }
-
- *new_val = sum;
- }
-
- new_state
-}
-
-fn apply_kernel_parallel(
- state: &[Complex<f64>],
- kernel: &Kernel,
- num_qubits: usize,
-) -> Vec<Complex<f64>> {
- let dim = 1 << num_qubits;
- let g = kernel.targets.len();
- let gate_dim = 1 << g;
-
- let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
- let mut non_target_mask: usize = (1 << num_qubits) - 1;
- for &pos in &target_bits {
- non_target_mask &= !(1 << pos);
- }
-
- (0..dim)
- .into_par_iter()
- .map(|i| {
- let mut target_idx = 0usize;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (i >> pos) & 1 == 1 {
- target_idx |= 1 << (g - 1 - k);
- }
- }
-
- let mut sum = complex!(0.0, 0.0);
-
- for j in 0..gate_dim {
- let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
- if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
- continue;
- }
-
- let mut source_idx = i & non_target_mask;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (j >> (g - 1 - k)) & 1 == 1 {
- source_idx |= 1 << pos;
- }
- }
-
- sum += gate_elem * state[source_idx];
- }
-
- sum
- })
- .collect()
-}
-
-pub struct KernelBuilder {
- num_qubits: usize,
-}
-
-impl KernelBuilder {
- pub fn new(num_qubits: usize) -> Self {
- Self { num_qubits }
- }
-
- pub fn num_qubits(&self) -> usize {
- self.num_qubits
- }
-}
-
-#[derive(Clone)]
-pub struct ExecutionLayer {
- pub kernels: Vec<Kernel>,
-}
-
-impl ExecutionLayer {
- pub fn new() -> Self {
- Self {
- kernels: Vec::new(),
- }
- }
-
- pub fn can_add(&self, kernel: &Kernel) -> bool {
- !self.kernels.iter().any(|k| k.shares_qubits(kernel))
- }
-
- pub fn add(&mut self, kernel: Kernel) {
- self.kernels.push(kernel);
- }
-
- pub fn affected_qubits(&self) -> HashSet<usize> {
- self.kernels
- .iter()
- .flat_map(|k| k.targets.iter().cloned())
- .collect()
- }
-}
-
-impl Default for ExecutionLayer {
- fn default() -> Self {
- Self::new()
- }
-}
-
-pub struct StructureAwareKernelBatch {
- kernels: Vec<Kernel>,
- layers: Vec<ExecutionLayer>,
- num_qubits: usize,
- optimised: bool,
-}
-
-impl StructureAwareKernelBatch {
- pub fn new(num_qubits: usize) -> Self {
- Self {
- kernels: Vec::new(),
- layers: Vec::new(),
- num_qubits,
- optimised: false,
- }
- }
-
- pub fn add(&mut self, kernel: Kernel) {
- self.kernels.push(kernel);
- self.optimised = false;
- }
-
- pub fn len(&self) -> usize {
- self.kernels.len()
- }
-
- pub fn is_empty(&self) -> bool {
- self.kernels.is_empty()
- }
-
- pub fn kernels(&self) -> &[Kernel] {
- &self.kernels
- }
-
- pub fn layers(&self) -> &[ExecutionLayer] {
- &self.layers
- }
-
- pub fn num_layers(&self) -> usize {
- self.layers.len()
- }
-
- pub fn optimise(&mut self) {
- if self.optimised || self.kernels.len() < 2 {
- return;
- }
-
- self.reorder_commuting_gates();
- self.multi_pass_fusion();
- self.build_execution_layers();
- self.optimised = true;
- }
-
- fn reorder_commuting_gates(&mut self) {
- let mut changed = true;
- let mut iterations = 0;
- const MAX_ITERATIONS: usize = 100;
-
- while changed && iterations < MAX_ITERATIONS {
- changed = false;
- iterations += 1;
-
- for i in 0..self.kernels.len().saturating_sub(1) {
- let current = &self.kernels[i];
- let next = &self.kernels[i + 1];
-
- if current.targets.len() == 1
- && next.targets.len() == 1
- && current.targets[0] != next.targets[0]
- && current.commutes_with(next)
- {
- for j in (i + 2)..self.kernels.len() {
- let candidate = &self.kernels[j];
-
- if candidate.targets.len() == 1
- && candidate.targets[0] == current.targets[0]
- {
- let can_move = (i + 1..j).all(|k| {
- let between = &self.kernels[k];
- !between.shares_qubits(current) || current.commutes_with(between)
- });
-
- if can_move && current.can_fuse_with(candidate) {
- let kernel_to_move = self.kernels.remove(j);
- self.kernels.insert(i + 1, kernel_to_move);
- changed = true;
- break;
- }
- }
- }
- }
- }
- }
- }
-
- fn multi_pass_fusion(&mut self) {
- let mut changed = true;
- let mut iterations = 0;
- const MAX_ITERATIONS: usize = 50;
-
- while changed && iterations < MAX_ITERATIONS {
- changed = false;
- iterations += 1;
-
- let mut new_kernels: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
- let mut i = 0;
-
- while i < self.kernels.len() {
- if i + 1 < self.kernels.len() {
- let current = &self.kernels[i];
- let next = &self.kernels[i + 1];
-
- if let Some(fused) = current.fuse(next) {
- new_kernels.push(fused);
- i += 2;
- changed = true;
- continue;
- }
- }
-
- new_kernels.push(self.kernels[i].clone());
- i += 1;
- }
-
- self.kernels = new_kernels;
- }
- }
-
- fn build_execution_layers(&mut self) {
- self.layers.clear();
-
- for kernel in &self.kernels {
- let mut placed = false;
-
- for layer in &mut self.layers {
- if layer.can_add(kernel) {
- layer.add(kernel.clone());
- placed = true;
- break;
- }
- }
-
- if !placed {
- let mut new_layer = ExecutionLayer::new();
- new_layer.add(kernel.clone());
- self.layers.push(new_layer);
- }
- }
- }
-
- pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
- for kernel in &self.kernels {
- *state = apply_kernel(state, kernel, self.num_qubits);
- }
- }
-
- pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
- for kernel in &self.kernels {
- *state = apply_kernel_parallel(state, kernel, self.num_qubits);
- }
- }
-
- pub fn execute_layered(&self, state: &mut Vec<Complex<f64>>) {
- for layer in &self.layers {
- for kernel in &layer.kernels {
- *state = apply_kernel(state, kernel, self.num_qubits);
- }
- }
- }
-
- pub fn execute_layered_parallel(&self, state: &mut Vec<Complex<f64>>) {
- for layer in &self.layers {
- for kernel in &layer.kernels {
- *state = apply_kernel_parallel(state, kernel, self.num_qubits);
- }
- }
- }
-
- pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
- 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<Complex<f64>>) {
- 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 stats(&self) -> KernelStats {
- let single_qubit = self.kernels.iter().filter(|k| k.targets.len() == 1).count();
- let two_qubit = self.kernels.iter().filter(|k| k.targets.len() == 2).count();
- let multi_qubit = self.kernels.iter().filter(|k| k.targets.len() > 2).count();
- let diagonal = self
- .kernels
- .iter()
- .filter(|k| k.gate_type == GateType::Diagonal)
- .count();
-
- KernelStats {
- total_kernels: self.kernels.len(),
- single_qubit,
- two_qubit,
- multi_qubit,
- diagonal,
- execution_layers: self.layers.len(),
- }
- }
-}
-
-#[derive(Debug, Clone)]
-pub struct KernelStats {
- pub total_kernels: usize,
- pub single_qubit: usize,
- pub two_qubit: usize,
- pub multi_qubit: usize,
- pub diagonal: usize,
- pub execution_layers: usize,
-}
-
-impl std::fmt::Display for KernelStats {
- fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
- write!(
- f,
- "Kernels: {} (1q: {}, 2q: {}, 3q+: {}, diag: {}), Layers: {}",
- self.total_kernels,
- self.single_qubit,
- self.two_qubit,
- self.multi_qubit,
- self.diagonal,
- self.execution_layers
- )
- }
-}
diff --git a/src/core/runtime.c b/src/core/runtime.c
index 9c95214..599ff3f 100644
--- a/src/core/runtime.c
+++ b/src/core/runtime.c
@@ -220,6 +220,8 @@ static void execute_kernels(struct PsiVector* state, const struct PsiKernel* ker
else
psi_apply_single_qubit_gate_simd(state->data, gate, kernel.targets[0], num_qubits);
}
+ else if (use_parallel)
+ psi_apply_kernel_parallel(state, kernel, num_qubits);
else
psi_apply_kernel(state, kernel, num_qubits);
}
diff --git a/src/core/runtime.rs b/src/core/runtime.rs
deleted file mode 100644
index 6ea7b76..0000000
--- a/src/core/runtime.rs
+++ /dev/null
@@ -1,585 +0,0 @@
-use super::{
- GateOp, Kernel, KernelBatch, QuantumGate, QuantumRegister, QuantumState,
- StructureAwareKernelBatch,
-};
-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,
-};
-use crate::maths::simd::{apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel};
-use crate::maths::vector::Vector;
-use crate::{complex, Complex, Matrix};
-use rayon::prelude::*;
-
-const PARALLEL_THRESHOLD: usize = 8;
-
-#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
-pub struct RuntimeConfig {
- pub parallel: bool,
- pub simd: bool,
- pub batched: bool,
- pub structure_aware: bool,
- pub parallel_threshold: usize,
-}
-
-impl RuntimeConfig {
- pub fn new() -> Self {
- Self {
- parallel: false,
- simd: false,
- batched: false,
- structure_aware: false,
- parallel_threshold: PARALLEL_THRESHOLD,
- }
- }
-
- pub fn parallel(mut self) -> Self {
- self.parallel = true;
- self
- }
-
- pub fn simd(mut self) -> Self {
- self.simd = true;
- self
- }
-
- pub fn batched(mut self) -> Self {
- self.batched = true;
- self
- }
-
- pub fn structure_aware(mut self) -> Self {
- self.structure_aware = true;
- self
- }
-
- pub fn with_threshold(mut self, threshold: usize) -> Self {
- self.parallel_threshold = threshold;
- self
- }
-
- pub fn optimal() -> Self {
- Self::new().structure_aware().simd().parallel()
- }
-
- pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
- let dim = 1 << num_qubits;
- let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
- state[0] = complex!(1.0, 0.0);
-
- let use_parallel = self.parallel && num_qubits >= self.parallel_threshold;
-
- if self.structure_aware {
- let mut batch = Runtime::build_structure_aware_batch(num_qubits, operations);
- batch.optimise();
- self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
- } else if self.batched {
- let mut batch = Runtime::build_kernel_batch(num_qubits, operations);
- batch.optimize();
- self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
- } else {
- let batch = Runtime::build_kernel_batch(num_qubits, operations);
- self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
- }
-
- QuantumState::new(state)
- }
-
- fn execute_kernels(
- &self,
- state: &mut Vec<Complex<f64>>,
- kernels: &[Kernel],
- num_qubits: usize,
- use_parallel: bool,
- ) {
- for kernel in kernels {
- if self.simd && kernel.targets.len() == 1 {
- let gate = matrix_to_2x2(&kernel.matrix);
- if use_parallel {
- apply_single_qubit_gate_simd_parallel(
- state,
- &gate,
- kernel.targets[0],
- num_qubits,
- );
- } else {
- apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], num_qubits);
- }
- } else if use_parallel {
- *state = apply_gate_parallel(state, &kernel.matrix, &kernel.targets, num_qubits);
- } else {
- *state = apply_kernel_direct(state, kernel, num_qubits);
- }
- }
- }
-}
-
-impl std::fmt::Display for RuntimeConfig {
- fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
- let mut features = Vec::new();
- if self.structure_aware {
- features.push("structure-aware");
- }
- if self.batched && !self.structure_aware {
- features.push("batched");
- }
- if self.simd {
- features.push("SIMD");
- }
- if self.parallel {
- features.push("parallel");
- }
- if features.is_empty() {
- features.push("basic");
- }
- write!(f, "Runtime[{}]", features.join("+"))
- }
-}
-
-#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
-pub enum Runtime {
- #[default]
- BasicRT,
- BasicRTMT,
- BatchedRT,
- BatchedRTMT,
- SimdRT,
- SimdRTMT,
- StructureAwareRT,
- StructureAwareMT,
- WFEvolution,
- WFEvolutionMT,
- GPUAccelerated,
- Custom(RuntimeConfig),
-}
-
-impl Runtime {
- pub fn custom() -> RuntimeConfig {
- RuntimeConfig::new()
- }
-
- pub fn optimal() -> RuntimeConfig {
- RuntimeConfig::optimal()
- }
-
- pub fn to_config(&self) -> RuntimeConfig {
- match self {
- Runtime::BasicRT => RuntimeConfig::new(),
- Runtime::BasicRTMT => RuntimeConfig::new().parallel(),
- Runtime::BatchedRT => RuntimeConfig::new().batched(),
- Runtime::BatchedRTMT => RuntimeConfig::new().batched().parallel(),
- Runtime::SimdRT => RuntimeConfig::new().batched().simd(),
- Runtime::SimdRTMT => RuntimeConfig::new().batched().simd().parallel(),
- Runtime::StructureAwareRT => RuntimeConfig::new().structure_aware().simd(),
- Runtime::StructureAwareMT => RuntimeConfig::new().structure_aware().simd().parallel(),
- Runtime::Custom(config) => *config,
- _ => RuntimeConfig::new(),
- }
- }
-
- pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
- match self {
- Runtime::BasicRT => Self::compute_basic(num_qubits, operations),
- Runtime::BasicRTMT => Self::compute_basic_mt(num_qubits, operations),
- Runtime::Custom(config) => config.compute(num_qubits, operations),
- Runtime::WFEvolution => {
- unimplemented!("WFEvolution (Schrödinger equation) runtime not yet implemented")
- }
- Runtime::WFEvolutionMT => {
- unimplemented!(
- "WFEvolutionMT (multi-threaded Schrödinger) runtime not yet implemented"
- )
- }
- Runtime::GPUAccelerated => {
- unimplemented!("GPUAccelerated runtime not yet implemented")
- }
- _ => self.to_config().compute(num_qubits, operations),
- }
- }
-
- pub fn build_kernel_batch(num_qubits: usize, operations: &[GateOp]) -> KernelBatch {
- let mut batch = KernelBatch::new(num_qubits);
-
- for op in operations {
- if let Some(kernel) = Self::op_to_kernel(op) {
- batch.add(kernel);
- }
- }
-
- batch
- }
-
- fn op_to_kernel(op: &GateOp) -> Option<Kernel> {
- let (matrix, targets, name): (Matrix<Complex<f64>>, Vec<usize>, &str) = match op {
- GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t], "H"),
- GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t], "X"),
- GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t], "Y"),
- GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t], "Z"),
- GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t], "S"),
- GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t], "T"),
- GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t], "Sdg"),
- GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t], "Tdg"),
- GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t], "Sx"),
- GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t], "Sxdg"),
- GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t], "Rx"),
- GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t], "Ry"),
- GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t], "Rz"),
- GateOp::P(t, theta) => (p_matrix(*theta), vec![*t], "P"),
- GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t], "U1"),
- GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t], "U2"),
- GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t], "U3"),
- GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t], "CNOT"),
- GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t], "CZ"),
- GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b], "SWAP"),
- GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t], "CRx"),
- GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t], "CRy"),
- GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t], "CRz"),
- GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t], "CP"),
- GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t], "CCNOT"),
- GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2], "CSWAP"),
- GateOp::Measure(_, _) => return None,
- GateOp::Custom(gate, tgts) => {
- let qg = gate.to_quantum_gate();
- (qg.matrix, tgts.clone(), "Custom")
- }
- };
-
- Some(Kernel::new(name, matrix, targets))
- }
-
- pub fn build_structure_aware_batch(
- num_qubits: usize,
- operations: &[GateOp],
- ) -> StructureAwareKernelBatch {
- let mut batch = StructureAwareKernelBatch::new(num_qubits);
-
- for op in operations {
- if let Some(kernel) = Self::op_to_kernel(op) {
- batch.add(kernel);
- }
- }
-
- batch
- }
-
- fn compute_basic(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
- let names: Vec<String> = (0..num_qubits).map(|i| format!("q{}", i)).collect();
- let leaked_names: &'static [String] = Box::leak(names.into_boxed_slice());
- let name_refs: Vec<&'static str> = leaked_names.iter().map(|s| s.as_str()).collect();
-
- let mut register = QuantumRegister::new(
- Box::leak(Box::new("circuit".to_string())).as_str(),
- &name_refs,
- );
-
- for op in operations {
- match op {
- // Clifford gates
- GateOp::H(t) => register.apply_gate(&HADAMARD, &[*t]),
- GateOp::X(t) => register.apply_gate(&PAULI_X, &[*t]),
- GateOp::Y(t) => register.apply_gate(&PAULI_Y, &[*t]),
- GateOp::Z(t) => register.apply_gate(&PAULI_Z, &[*t]),
- GateOp::S(t) => register.apply_gate(&S_GATE, &[*t]),
- GateOp::CNOT(c, t) => register.apply_gate(&CNOT, &[*c, *t]),
- GateOp::CZ(c, t) => register.apply_gate(&CZ, &[*c, *t]),
- 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 {
- name: "Rx",
- matrix: rx_matrix(*theta),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::Ry(t, theta) => {
- let gate = QuantumGate {
- name: "Ry",
- matrix: ry_matrix(*theta),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::Rz(t, theta) => {
- let gate = QuantumGate {
- name: "Rz",
- matrix: rz_matrix(*theta),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::P(t, theta) => {
- let gate = QuantumGate {
- name: "P",
- matrix: p_matrix(*theta),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::U1(t, lambda) => {
- let gate = QuantumGate {
- name: "U1",
- matrix: u1_matrix(*lambda),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::U2(t, phi, lambda) => {
- let gate = QuantumGate {
- name: "U2",
- matrix: u2_matrix(*phi, *lambda),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
- GateOp::U3(t, theta, phi, lambda) => {
- let gate = QuantumGate {
- name: "U3",
- matrix: u3_matrix(*theta, *phi, *lambda),
- num_qubits: 1,
- };
- register.apply_gate(&gate, &[*t]);
- }
-
- // Controlled parametric gates
- GateOp::CRx(c, t, theta) => {
- let gate = QuantumGate {
- name: "CRx",
- matrix: crx_matrix(*theta),
- num_qubits: 2,
- };
- register.apply_gate(&gate, &[*c, *t]);
- }
- GateOp::CRy(c, t, theta) => {
- let gate = QuantumGate {
- name: "CRy",
- matrix: cry_matrix(*theta),
- num_qubits: 2,
- };
- register.apply_gate(&gate, &[*c, *t]);
- }
- GateOp::CRz(c, t, theta) => {
- let gate = QuantumGate {
- name: "CRz",
- matrix: crz_matrix(*theta),
- num_qubits: 2,
- };
- register.apply_gate(&gate, &[*c, *t]);
- }
- GateOp::CP(c, t, theta) => {
- let gate = QuantumGate {
- name: "CP",
- matrix: cp_matrix(*theta),
- num_qubits: 2,
- };
- register.apply_gate(&gate, &[*c, *t]);
- }
-
- // Measurement and custom gates
- GateOp::Measure(_, _) => {}
- GateOp::Custom(gate, targets) => {
- let quantum_gate = gate.to_quantum_gate();
- register.apply_gate(&quantum_gate, targets);
- }
- }
- }
-
- register.get_state()
- }
-
- fn compute_basic_mt(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
- // For small circuits, fall back to single-threaded (overhead not worth it)
- if num_qubits < PARALLEL_THRESHOLD {
- return Self::compute_basic(num_qubits, operations);
- }
-
- let dim = 1 << num_qubits;
-
- // Initialize state to |0...0⟩
- let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
- state[0] = complex!(1.0, 0.0);
-
- for op in operations {
- let (gate_matrix, targets): (Matrix<Complex<f64>>, Vec<usize>) = match op {
- // Clifford gates
- GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t]),
- GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t]),
- GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t]),
- GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t]),
- GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t]),
- GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t]),
- GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t]),
- 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]),
- GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t]),
- GateOp::P(t, theta) => (p_matrix(*theta), vec![*t]),
- 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) => {
- let quantum_gate = custom_gate.to_quantum_gate();
- state = apply_gate_parallel(&state, &quantum_gate.matrix, tgts, num_qubits);
- continue;
- }
- };
-
- state = apply_gate_parallel(&state, &gate_matrix, &targets, num_qubits);
- }
-
- QuantumState::new(state)
- }
-}
-
-/// Apply a gate to the state vector in parallel using sparse application
-/// This is O(2^n * 2^g) instead of O(2^2n) for full matrix multiplication
-fn apply_gate_parallel(
- state: &[Complex<f64>],
- gate_matrix: &Matrix<Complex<f64>>,
- targets: &[usize],
- num_qubits: usize,
-) -> Vec<Complex<f64>> {
- let dim = 1 << num_qubits;
- let g = targets.len();
- let gate_dim = 1 << g;
-
- // Convert target qubit indices to bit positions (from MSB)
- let target_bits: Vec<usize> = targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
- // Create a mask for non-target qubits
- let mut non_target_mask: usize = (1 << num_qubits) - 1;
- for &pos in &target_bits {
- non_target_mask &= !(1 << pos);
- }
-
- // Parallel computation of new state
- let new_state: Vec<Complex<f64>> = (0..dim)
- .into_par_iter()
- .map(|i| {
- // Extract the target qubit bits from index i
- let mut target_idx = 0usize;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (i >> pos) & 1 == 1 {
- target_idx |= 1 << (g - 1 - k);
- }
- }
-
- // Compute the contribution to state[i]
- let mut sum = complex!(0.0, 0.0);
-
- // For each possible input state that could contribute
- for j in 0..gate_dim {
- // Get the gate matrix element
- let gate_elem = gate_matrix.data[target_idx * gate_dim + j];
-
- // Skip if zero (sparse optimization)
- if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
- continue;
- }
-
- // Compute the source index by replacing target bits in i with bits from j
- let mut source_idx = i & non_target_mask;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (j >> (g - 1 - k)) & 1 == 1 {
- source_idx |= 1 << pos;
- }
- }
-
- sum += gate_elem * state[source_idx];
- }
-
- sum
- })
- .collect();
-
- new_state
-}
-
-fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] {
- [
- [matrix.data[0], matrix.data[1]],
- [matrix.data[2], matrix.data[3]],
- ]
-}
-
-fn apply_kernel_direct(
- state: &[Complex<f64>],
- kernel: &Kernel,
- num_qubits: usize,
-) -> Vec<Complex<f64>> {
- let dim = 1 << num_qubits;
- let g = kernel.targets.len();
- let gate_dim = 1 << g;
-
- let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
- let mut non_target_mask: usize = (1 << num_qubits) - 1;
- for &pos in &target_bits {
- non_target_mask &= !(1 << pos);
- }
-
- let mut new_state = vec![complex!(0.0, 0.0); dim];
-
- for (i, new_val) in new_state.iter_mut().enumerate() {
- let mut target_idx = 0usize;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (i >> pos) & 1 == 1 {
- target_idx |= 1 << (g - 1 - k);
- }
- }
-
- let mut sum = complex!(0.0, 0.0);
-
- for j in 0..gate_dim {
- let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
- if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
- continue;
- }
-
- let mut source_idx = i & non_target_mask;
- for (k, &pos) in target_bits.iter().enumerate() {
- if (j >> (g - 1 - k)) & 1 == 1 {
- source_idx |= 1 << pos;
- }
- }
-
- sum += gate_elem * state[source_idx];
- }
-
- *new_val = sum;
- }
-
- new_state
-}