Added OpenCL Backend

Cargo.toml

@@ -4,3 +4,11 @@ version = "1.0.0"
 authors = ["Adam Kelly <adamkelly2201@gmail.com>"]
 
 [dependencies]
+num-complex = "0.1.0"
+ocl = "0.18.0"
+rand = "0.5.1"
+
+[[example]]
+name = "bell"
+doc = false
+path = "examples/bell-state.rs"

examples/bell-state.rs

@@ -0,0 +1,15 @@
+extern crate qcgpu;
+
+use qcgpu::Simulator;
+
+fn main() {
+    println!("Creating Bell State");
+
+    let mut sim = Simulator::new_opencl(2);
+
+    sim.h(0);
+    sim.cx(0, 1);
+
+    println!("Measurment Results:");
+    println!("{}", sim.measure());
+}

src/backends/mod.rs

@@ -1,3 +1,3 @@
 mod opencl;
 
-pub use self::opencl::OpenCL;
+pub use self::opencl::OpenCL;

src/backends/opencl/kernel.cl

@@ -0,0 +1,277 @@
+typedef float2 complex_f;
+
+/*
+ * Addition of two complex numbers:
+ *
+ * a + b = (Re(a) + Re(b)) + i(Im(a) + Im(b))
+ */
+static complex_f add(complex_f a, complex_f b)
+{
+    return (complex_f)(a.x + b.x, a.y + b.y);
+}
+
+/*
+ * Negation of a complex numbers:
+ *
+ * -a = -(Re(a) - i(Im(a))
+ */
+static complex_f neg(complex_f a)
+{
+    return (complex_f)(-a.x, -a.y);
+}
+
+/*
+ * Multiplication of two complex numbers:
+ *
+ * a * b =
+ *   ((Re(a) * Re(b)) - (Im(a) * Im(b)))
+ * + ((Im(a) * Re(b)) + (Re(a) * Im(b)))i
+ */
+static complex_f mul(complex_f a, complex_f b)
+{
+    return (complex_f)(
+        (a.x * b.x) - (a.y * b.y),
+        (a.y * b.x) + (a.x * b.y));
+}
+/**
+ * Absolute value of a complex number
+ *
+ * |a| = √(Re(a)^2 + Im(a)^2)
+ */
+static float complex_abs(complex_f a)
+{
+    return sqrt((a.x * a.x) + (a.y * a.y));
+}
+
+static complex_f cexp(float a) {
+    return (complex_f)(cos(a), sin(a));
+}
+/*
+ * Applies a single qubit gate to the register.
+ * The gate matrix must be given in the form:
+ *
+ *  A B
+ *  C D
+ */
+__kernel void apply_gate(
+    __global complex_f *const amplitudes,
+    __global complex_f *amps,
+    uint target,
+    complex_f A,
+    complex_f B,
+    complex_f C,
+    complex_f D)
+{
+    uint const state = get_global_id(0);
+    complex_f const amp = amplitudes[state];
+
+    uint const zero_state = state & (~(1 << target));
+    uint const one_state = state | (1 << target);
+
+    uint const bit_val = (((1 << target) & state) > 0) ? 1 : 0;
+
+    if (bit_val == 0)
+    {
+        // Bitval = 0
+
+        amps[state] = add(mul(A, amp), mul(B, amplitudes[one_state]));
+    }
+    else
+    {
+        amps[state] = add(mul(D, amp), mul(C, amplitudes[zero_state]));
+    }
+}
+
+/*
+ * Applies a controlled single qubit gate to the register.
+ */
+__kernel void apply_controlled_gate(
+    __global complex_f *const amplitudes,
+    __global complex_f *amps,
+    uint control,
+    uint target,
+    complex_f A,
+    complex_f B,
+    complex_f C,
+    complex_f D)
+{
+    uint const state = get_global_id(0);
+    complex_f const amp = amplitudes[state];
+
+    uint const zero_state = state & (~(1 << target));
+    uint const one_state = state | (1 << target);
+
+    uint const bit_val = (((1 << target) & state) > 0) ? 1 : 0;
+    uint const control_val = (((1 << control) & state) > 0) ? 1 : 0;
+
+    if (control_val == 0)
+    {
+        // Control is 0, don't apply gate
+        amps[state] = amp;
+    }
+    else
+    {
+        // control is 1, apply gate.
+        if (bit_val == 0)
+        {
+            // Bitval = 0
+            amps[state] = add(mul(A, amp), mul(B, amplitudes[one_state]));
+        }
+        else
+        {
+            amps[state] = add(mul(D, amp), mul(C, amplitudes[zero_state]));
+        }
+    }
+}
+
+/*
+ * Applies a controlled-controlled single qubit gate to the register.
+ */
+__kernel void apply_controlled_controlled_gate(
+    __global complex_f *const amplitudes,
+    __global complex_f *amps,
+    uint control1,
+    uint control2,
+    uint target,
+    complex_f A,
+    complex_f B,
+    complex_f C,
+    complex_f D)
+{
+    uint const state = get_global_id(0);
+    complex_f const amp = amplitudes[state];
+
+    uint const zero_state = state & (~(1 << target));
+    uint const one_state = state | (1 << target);
+
+    uint const bit_val = (((1 << target) & state) > 0) ? 1 : 0;
+    uint const control1_val = (((1 << control1) & state) > 0) ? 1 : 0;
+    uint const control2_val = (((1 << control2) & state) > 0) ? 1 : 0;
+
+    if (control1_val == 0 || control2_val == 0)
+    {
+        // Control is 0, don't apply gate
+        amps[state] = amp;
+    }
+    else
+    {
+        // control is 1, apply gate.
+        if (bit_val == 0)
+        {
+            // Bitval = 0
+            amps[state] = add(mul(A, amp), mul(B, amplitudes[one_state]));
+        }
+        else
+        {
+            amps[state] = add(mul(D, amp), mul(C, amplitudes[zero_state]));
+        }
+    }
+}
+
+/*
+ * Swaps the states of two qubits in the register
+ */
+__kernel void swap(
+    __global complex_f *const amplitudes,
+    __global complex_f *amps,
+    uint first_qubit,
+    uint second_qubit)
+{
+    uint const state = get_global_id(0);
+
+    uint const first_bit_mask = 1 << first_qubit;
+    uint const second_bit_mask = 1 << second_qubit;
+
+    uint const new_second_bit = ((state & first_bit_mask) >> first_qubit) << second_qubit;
+    uint const new_first_bit = ((state & second_bit_mask) >> second_qubit) << first_qubit;
+
+    uint const new_state = (state & !first_bit_mask & !second_bit_mask) | new_first_bit | new_second_bit;
+
+    amps[new_state] = amplitudes[state];
+}
+
+static uint pow_mod(uint x, uint y, uint n)
+{
+    uint r = 1;
+    while (y > 0)
+    {
+        if ((y & 1) == 1)
+        {
+            r = r * x % n;
+        }
+        y /= 2;
+        x = x * x % n;
+    }
+
+    return r;
+}
+
+/*
+ *  Calculates f(a) = x^a mod N
+ */
+__kernel void apply_pow_mod(
+    __global complex_f *const amplitudes,
+    __global complex_f *amps,
+    uint x,
+    uint n,
+    uint input_width,
+    uint output_width)
+{
+    uint input_bit_range_from = output_width;
+    uint input_bit_range_to = output_width + input_width - 1;
+    uint target_bit_range_from = 0;
+    uint target_bit_range_to = output_width - 1;
+
+    uint high_bit_mask = (1 << (input_bit_range_to + 1)) - 1;
+    uint target_bit_mask = ((1 << (1 + target_bit_range_to - target_bit_range_from)) - 1) << target_bit_range_from;
+
+    uint const state = get_global_id(0);
+
+    uint input = (state & high_bit_mask) >> input_bit_range_from;
+    uint result = (pow_mod(x, input, n) << target_bit_range_from) & target_bit_mask;
+    uint result_state = state ^ result;
+
+    if (result_state == state)
+    {
+        amps[state] = amplitudes[state];
+    }
+    else
+    {
+        amps[state] = amplitudes[result_state];
+        amps[result_state] = amplitudes[state];
+    }
+
+    amps[result_state] = amplitudes[state];
+}
+
+/**
+ * Calculates The Probabilities Of A State Vector
+ */
+__kernel void calculate_probabilities(
+    __global complex_f *const amplitudes,
+    __global float *probabilities)
+{
+    uint const state = get_global_id(0);
+    complex_f amp = amplitudes[state];
+
+    probabilities[state] = complex_abs(mul(amp, amp));
+}
+
+/**
+ * Initializes a register to the value 1|0..100...0>
+ *                                          ^ target
+ */
+__kernel void initialize_register(
+    __global complex_f *amplitudes,
+    uint const target)
+{
+    uint const state = get_global_id(0);
+    if (state == target)
+    {
+        amplitudes[state] = (complex_f)(1, 0);
+    }
+    else
+    {
+        amplitudes[state] = (complex_f)(0, 0);
+    }
+}

src/backends/opencl/mod.rs

@@ -1,26 +1,140 @@
 use traits::Backend;
 use gate::Gate;
 
-use std::collections::HashMap;
+use ocl::{Buffer, MemFlags, ProQue};
+use num_complex::{Complex, Complex32};
+use rand::random;
+
+// OpenCL Kernel
+pub static KERNEL: &'static str = include_str!("kernel.cl");
 
 pub struct OpenCL {
+    /// OpenCL Buffer for the state vector
+    pub buffer: Buffer<Complex<f32>>,
+    pro_que: ProQue,
+}
+
+impl OpenCL {
+    pub fn new(num_qubits: u8) -> OpenCL {
+        // How many amplitudes needed?
+        let num_amps = 2_usize.pow(u32::from(num_qubits)) as usize;
+
+        let ocl_pq = ProQue::builder()
+            .src(KERNEL)
+            .device(0)
+            .dims(num_amps)
+            .build()
+            .expect("Error Building ProQue");
+
+        let buffer: Buffer<Complex32> = Buffer::builder()
+            .queue(ocl_pq.queue().clone())
+            .flags(MemFlags::new().read_write())
+            .len(num_amps)
+            .build()
+            .expect("Source Buffer");
+
+        let apply = ocl_pq
+            .kernel_builder("initialize_register")
+            .arg(&buffer)
+            .arg(0)
+            .build()
+            .unwrap();
+
+        unsafe {
+            apply.enq().unwrap();
+        }
+
+        OpenCL {
+            pro_que: ocl_pq,
+            buffer,
+        }
+    }
+
+    fn get_probabilities(&self) -> Vec<f32> {
+        let result_buffer: Buffer<f32> = self.pro_que.create_buffer().unwrap();
+
+        let apply = self.pro_que
+            .kernel_builder("calculate_probabilities")
+            .arg(&self.buffer)
+            .arg(&result_buffer)
+            .build()
+            .unwrap();
 
+        unsafe {
+            apply.enq().unwrap();
+        }
+
+        let mut vec_result = vec![0.0f32; self.buffer.len()];
+        result_buffer.read(&mut vec_result).enq().unwrap();
+
+        vec_result
+    }
 }
 
 impl Backend for OpenCL {
-  fn apply_gate(&mut self, gate: Gate, target: u8) {
+    fn apply_gate(&mut self, gate: Gate, target: u8) {
+        // create a temporary vector with the source buffer
+        let result_buffer: Buffer<Complex32> = self.pro_que.create_buffer().unwrap();
+
+        let apply = self.pro_que
+            .kernel_builder("apply_gate")
+            .arg(&self.buffer)
+            .arg(&result_buffer)
+            .arg(i32::from(target))
+            .arg(gate.a)
+            .arg(gate.b)
+            .arg(gate.c)
+            .arg(gate.d)
+            .build()
+            .unwrap();
+
+        unsafe {
+            apply.enq().unwrap();
+        }
+
+        self.buffer = result_buffer;
+    }
 
-  }
+    fn apply_controlled_gate(&mut self, gate: Gate, control: u8, target: u8) {
+        let result_buffer: Buffer<Complex32> = self.pro_que.create_buffer().unwrap();
 
-  fn apply_controlled_gate(&mut self, gate: Gate, control: u8, target: u8) {
+        let apply = self.pro_que
+            .kernel_builder("apply_controlled_gate")
+            .arg(&self.buffer)
+            .arg(&result_buffer)
+            .arg(control)
+            .arg(target)
+            .arg(gate.a)
+            .arg(gate.b)
+            .arg(gate.c)
+            .arg(gate.d)
+            .build()
+            .unwrap();
 
-  } 
+        unsafe {
+            apply.enq().unwrap();
+        }
 
-  fn measure(&mut self) -> u8 {
-    1
-  }
+        self.buffer = result_buffer;
+    }
 
-  fn measure_many(&mut self, iters: u64) -> HashMap<u8, u64> {
-    HashMap::new()
-  }
-}
+    fn measure(&self) -> u8 {
+        let probabilities = self.get_probabilities();
+
+        let mut key = random::<f32>();
+        if key > 1.0 {
+            key %= 1.0;
+        }
+
+        let mut i = 0;
+        while i < probabilities.len() {
+            key -= probabilities[i];
+            if key <= 0.0 {
+                break;
+            }
+            i += 1;
+        }
+
+        i as u8
+    }
+}

src/gate.rs

@@ -1,2 +1,46 @@
-pub struct Gate { 
+use num_complex::Complex32;
+use std::fmt;
+use std::f32::consts::FRAC_1_SQRT_2;
+
+pub struct Gate {
+    pub a: Complex32,
+    pub b: Complex32,
+    pub c: Complex32,
+    pub d: Complex32,
+}
+
+impl fmt::Display for Gate {
+    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
+        write!(f, "[[{}, {}], [{}, {}]]", self.a, self.b, self.c, self.d)
+    }
+}
+
+/// Hadamard Gate
+///
+/// [0.70710678118, 0.70710678118]
+///
+/// [0.70710678118, -0.70710678118]
+#[inline]
+pub fn h() -> Gate {
+    Gate {
+        a: Complex32::new(FRAC_1_SQRT_2, 0.0),
+        b: Complex32::new(FRAC_1_SQRT_2, 0.0),
+        c: Complex32::new(FRAC_1_SQRT_2, 0.0),
+        d: Complex32::new(-FRAC_1_SQRT_2, 0.0),
+    }
+}
+
+/// Pauli X / NOT Gate
+///
+/// [0, 1]
+///
+/// [1, 0]
+#[inline]
+pub fn x() -> Gate {
+    Gate {
+        a: Complex32::new(0.0, 0.0),
+        b: Complex32::new(1.0, 0.0),
+        c: Complex32::new(1.0, 0.0),
+        d: Complex32::new(0.0, 0.0),
+    }
 }

src/lib.rs

@@ -1,7 +1,38 @@
+extern crate num_complex;
+extern crate ocl;
+extern crate rand;
+
 pub mod gate;
 pub mod traits;
 pub mod backends;
 
+use backends::OpenCL;
+use gate::{h, x};
+
 pub struct Simulator {
-    backend: Box<traits::Backend>
+    backend: Box<traits::Backend>,
 }
+
+impl Simulator {
+    pub fn new_opencl(num_qubits: u8) -> Simulator {
+        Simulator {
+            backend: Box::new(OpenCL::new(num_qubits))
+        }
+    }
+
+    pub fn x(&mut self, target: u8) {
+        self.backend.apply_gate(x(), target);
+    }
+
+    pub fn h(&mut self, target: u8) {
+        self.backend.apply_gate(h(), target);
+    }
+
+    pub fn cx(&mut self, control: u8, target: u8) {
+        self.backend.apply_controlled_gate(x(), control, target);
+    }
+
+    pub fn measure(&self) -> u8 {
+        self.backend.measure()
+    }
+}

src/traits.rs

@@ -1,9 +1,7 @@
 use gate::Gate;
-use std::collections::HashMap;
 
 pub trait Backend {
     fn apply_gate(&mut self, gate: Gate, target: u8);
     fn apply_controlled_gate(&mut self, gate: Gate, control: u8, target: u8);
-    fn measure(&mut self) -> u8;
-    fn measure_many(&mut self, iters: u64) -> HashMap<u8, u64>;
+    fn measure(&self) -> u8;
 }