Added distributed simulation code

CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 3.9)
+
+project(tangle)
+
+find_package(MPI REQUIRED)
+
+# Add headers
+include_directories(include)
+
+file(GLOB SOURCES "src/*.c")
+
+add_executable(main ${SOURCES})
+
+target_link_libraries(main PRIVATE MPI::MPI_C)

include/communication.h

@@ -0,0 +1,11 @@
+#ifndef __COMMUNICATION_H
+#define __COMMUNICATION_H
+
+#include "tangle.h"
+
+void send_top(TangleState state, int node);
+void recieve_top(TangleState state, int node);
+void send_bottom(TangleState state, int node);
+void recieve_bottom(TangleState state, int node);
+
+#endif

include/tangle.h

@@ -0,0 +1,69 @@
+#ifndef __TANGLE_H
+#define __TANGLE_H
+
+#include <complex.h>
+
+// Common Types 
+typedef float _Complex cfloat;
+typedef long long int llint;
+
+// Helper Method
+#define cfloat(r, i) ((float)(r) + ((float)(i)) * I)
+
+// Tangle Types
+typedef struct {
+    int rank;
+    int nodes;
+} TangleEnvironment;
+
+typedef struct {
+    int num_qubits;
+    cfloat *amps;
+
+    int rank;
+    int nodes;
+    int k;
+    int m;
+
+    llint node_amps;
+    llint temp_amps;
+} TangleState;
+
+typedef struct {
+    cfloat A;
+    cfloat B;
+    cfloat C;
+    cfloat D;
+} TangleGate;
+
+// Environment/Distributed Functions
+TangleEnvironment initialize_environment();
+void destroy_environemt(TangleEnvironment env);
+
+// Quantum State Functions
+TangleState create_state(int num_qubits, TangleEnvironment env);
+void print_state(TangleState state);
+
+// Quantum Gates
+void apply_gate(TangleState state, int target, TangleGate u);
+void apply_diagonal_gate(TangleState state, int target, TangleGate u);
+void apply_antidiagonal_gate(TangleState state, int target, TangleGate u);
+void apply_controlled_gate(TangleState state, int control, int target,
+                           TangleGate u);
+void apply_controlled_gate(TangleState state, int control, int target,
+                           TangleGate u);
+void apply_controlled_diagonal_gate(TangleState state, int control, int target,
+                                    TangleGate u);
+void apply_controlled_antidiagonal_gate(TangleState state, int control,
+                                        int target, TangleGate u);
+
+void X(TangleState state, int target);
+void H(TangleState state, int target);
+void Z(TangleState state, int target);
+void CZ(TangleState state, int control, int target);
+void CX(TangleState state, int control, int target);
+
+// Measurmement
+llint measure(TangleState state);
+
+#endif

src/communication.c

@@ -0,0 +1,48 @@
+#include "mpi.h"
+#include "tangle.h"
+
+void communication_helper(TangleState state, int node, llint idx_1, llint idx_2,
+                          int tag_1, int tag_2) {
+  // There is a limit to the amount of data that can be transfered
+  llint message_size = 1LL << 29; // 28 for a double, 27 for a long quad
+  
+  // It will always be a state.temp_amps sized block transfered.
+  // This can be decomposed if needed
+  if (state.temp_amps < message_size) {
+    message_size = state.temp_amps;
+  }
+
+  // Message size is a power of two, so we can break up the message
+
+  for (llint i = 0; i < state.temp_amps / message_size; i++) {
+    MPI_Sendrecv(&state.amps[idx_1 + (i * message_size)], message_size, MPI_COMPLEX, node, tag_1,
+               &state.amps[idx_2 + (i * message_size)], message_size, MPI_COMPLEX, node, tag_2,
+               MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+  }
+
+  
+}
+
+// Send the top half of the state,
+// Recieve into the temporary storage
+void send_top(TangleState state, int node) {
+  communication_helper(state, node, 0, state.node_amps, 1, 2);
+}
+
+// Send the temporary storage,
+// recieve into the top half
+void recieve_top(TangleState state, int node) {
+  communication_helper(state, node, state.node_amps, 0, 3, 4);
+}
+
+// Send the bottom half of the state
+// Recieve into the temporary state
+void send_bottom(TangleState state, int node) {
+  communication_helper(state, node, state.temp_amps, state.node_amps, 2, 1);
+}
+
+// Send the temporary storage
+// Receive into the bottom half
+void recieve_bottom(TangleState state, int node) {
+  communication_helper(state, node, state.node_amps, state.temp_amps, 4, 3);
+}

src/distributed.c

@@ -0,0 +1,31 @@
+#include "tangle.h"
+
+#include "mpi.h"
+
+TangleEnvironment initialize_environment() {
+  TangleEnvironment env;
+
+  int rank, nodes, initialized;
+  MPI_Initialized(&initialized);
+
+  if (!initialized) {
+    MPI_Init(0, 0);
+  }
+
+  MPI_Comm_size(MPI_COMM_WORLD, &nodes);
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  env.rank = rank;
+  env.nodes = nodes;
+
+  return env;
+}
+
+void destroy_environment(TangleEnvironment env) {
+  int finalized;
+  MPI_Finalized(&finalized);
+
+  if (!finalized) {
+    MPI_Finalize();
+  }
+}

src/gate_application.c

@@ -0,0 +1,268 @@
+#include "mpi.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "communication.h"
+#include "tangle.h"
+
+llint is_bit_clear(llint n, llint pos) { return n & (1 << pos); }
+
+void apply_gate(TangleState state, int target, TangleGate u) {
+  if (target < state.m) {
+    // Communication isn't required
+    llint stride = 1LL << (target + 1LL);
+    for (llint i = 0; i < state.node_amps; i += stride) {
+      for (llint j = i; j < i + (1LL << target); j++) {
+        llint zero = j;
+        llint one = j + (1LL << target);
+
+        cfloat zero_amp = state.amps[zero];
+        cfloat one_amp = state.amps[one];
+
+        state.amps[zero] = (u.A * zero_amp) + (u.B * one_amp);
+        state.amps[one] = (u.C * zero_amp) + (u.D * one_amp);
+      }
+    }
+  } else {
+    // Communication is required
+    llint stride = 1LL << (target - state.m);
+    llint global_state = state.rank * state.node_amps;
+
+    llint node_i, node_j;
+    if (is_bit_clear(global_state, target) == 0) {
+      node_i = state.rank;
+      node_j = node_i + stride;
+    } else {
+      node_j = state.rank;
+      node_i = node_j - stride;
+    }
+
+    if (node_i == state.rank) {
+      // s1 --> dt
+      // d2 --> st
+      send_top(state, node_j);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i + state.temp_amps;
+        llint one = i;
+
+        cfloat zero_amp = (&state.amps[0])[zero];
+        cfloat one_amp = (&state.amps[state.node_amps])[one];
+
+        (&state.amps[0])[zero] = (u.A * zero_amp) + (u.B * one_amp);
+        (&state.amps[state.node_amps])[one] =
+            (u.C * zero_amp) + (u.D * one_amp);
+      }
+
+      // st --> d2
+      // dt --> s1
+      recieve_top(state, node_j);
+    } else {
+      // s2 --> dt
+      // d1 --> st
+      send_bottom(state, node_i);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i;
+        llint one = i;
+
+        cfloat zero_amp = (&state.amps[state.node_amps])[zero];
+        cfloat one_amp = (&state.amps[0])[one];
+
+        (&state.amps[state.node_amps])[zero] =
+            (u.A * zero_amp) + (u.B * one_amp);
+        (&state.amps[0])[one] = (u.C * zero_amp) + (u.D * one_amp);
+      }
+
+      // st --> d1
+      // dt --> s2
+      recieve_bottom(state, node_i);
+    }
+  }
+  MPI_Barrier(MPI_COMM_WORLD);
+}
+
+void apply_diagonal_gate(TangleState state, int target, TangleGate u) {
+  // Communication can always be avoided
+  for (llint i = 0; i < state.node_amps; i++) {
+    if (is_bit_clear(i, target) == 0) {
+      state.amps[i] = u.A * state.amps[i];
+    } else {
+      state.amps[i] = u.D * state.amps[i];
+    }
+  }
+
+  MPI_Barrier(MPI_COMM_WORLD);
+}
+
+void apply_antidiagonal_gate(TangleState state, int target, TangleGate u) {
+  if (target < state.m) {
+    // Communication isn't required
+    llint stride = 1LL << (target + 1LL);
+    for (llint i = 0; i < state.node_amps; i += stride) {
+      for (llint j = i; j < i + (1LL << target); j++) {
+        llint zero = j;
+        llint one = j + (1LL << target);
+
+        cfloat zero_amp = state.amps[zero];
+        cfloat one_amp = state.amps[one];
+
+        state.amps[zero] = u.B * one_amp;
+        state.amps[one] = u.C * zero_amp;
+      }
+    }
+  } else {
+    // Communication is required
+    llint stride = 1LL << (target - state.m);
+    llint global_state = state.rank * state.node_amps;
+
+    llint node_i, node_j;
+    if (is_bit_clear(global_state, target) == 0) {
+      node_i = state.rank;
+      node_j = node_i + stride;
+    } else {
+      node_j = state.rank;
+      node_i = node_j - stride;
+    }
+
+    if (node_i == state.rank) {
+      // s1 --> dt
+      // d2 --> st
+      send_top(state, node_j);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i + state.temp_amps;
+        llint one = i;
+
+        cfloat zero_amp = (&state.amps[0])[zero];
+        cfloat one_amp = (&state.amps[state.node_amps])[one];
+
+        (&state.amps[0])[zero] = u.B * one_amp;
+        (&state.amps[state.node_amps])[one] = u.C * zero_amp;
+      }
+
+      // st --> d2
+      // dt --> s1
+      recieve_top(state, node_j);
+    } else {
+      // s2 --> dt
+      // d1 --> st
+      send_bottom(state, node_i);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i;
+        llint one = i;
+
+        cfloat zero_amp = (&state.amps[state.node_amps])[zero];
+        cfloat one_amp = (&state.amps[0])[one];
+
+        (&state.amps[state.node_amps])[zero] = u.B * one_amp;
+        (&state.amps[0])[one] = u.C * zero_amp;
+      }
+
+      // st --> d1
+      // dt --> s2
+      recieve_bottom(state, node_i);
+    }
+  }
+  MPI_Barrier(MPI_COMM_WORLD);
+}
+
+void apply_controlled_gate(TangleState state, int control, int target,
+                           TangleGate u) {
+  if (target < state.m) {
+    // Communication isn't required
+    llint stride = 1LL << (target + 1LL);
+    // TODO: Better stride for controlled gates
+    for (llint i = 0; i < state.node_amps; i += stride) {
+      for (llint j = i; j < i + (1LL << target); j++) {
+        if (is_bit_clear(j, control) != 0) {
+          llint zero = j;
+          llint one = j + (1LL << target);
+
+          cfloat zero_amp = state.amps[zero];
+          cfloat one_amp = state.amps[one];
+
+          state.amps[zero] = (u.A * zero_amp) + (u.B * one_amp);
+          state.amps[one] = (u.C * zero_amp) + (u.D * one_amp);
+        }
+      }
+    }
+  } else {
+    // Communication is required
+    llint stride = 1LL << (target - state.m);
+    llint global_state = state.rank * state.node_amps;
+
+    llint node_i, node_j;
+
+    if (is_bit_clear(global_state, target) == 0) {
+      node_i = state.rank;
+      node_j = node_i + stride;
+    } else {
+      node_j = state.rank;
+      node_i = node_j - stride;
+    }
+
+    if (node_i == state.rank) {
+      // s1 --> dt
+      // d2 --> st
+      send_top(state, node_j);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i + state.temp_amps;
+        llint one = i;
+
+        if (is_bit_clear(zero, control) != 0) {
+          cfloat zero_amp = (&state.amps[0])[zero];
+          cfloat one_amp = (&state.amps[state.node_amps])[one];
+
+          (&state.amps[0])[zero] = (u.A * zero_amp) + (u.B * one_amp);
+          (&state.amps[state.node_amps])[one] =
+              (u.C * zero_amp) + (u.D * one_amp);
+        }
+      }
+
+      // st --> d2
+      // dt --> s1
+      recieve_top(state, node_j);
+    } else {
+      // s2 --> dt
+      // d1 --> st
+      send_bottom(state, node_i);
+
+      for (llint i = 0; i < state.temp_amps; i++) {
+        llint zero = i;
+        llint one = i;
+
+        if (is_bit_clear(zero, control) != 0) {
+
+          cfloat zero_amp = (&state.amps[state.node_amps])[zero];
+          cfloat one_amp = (&state.amps[0])[one];
+
+          (&state.amps[state.node_amps])[zero] =
+              (u.A * zero_amp) + (u.B * one_amp);
+          (&state.amps[0])[one] = (u.C * zero_amp) + (u.D * one_amp);
+        }
+      }
+      // st --> d1
+      // dt --> s2
+      recieve_bottom(state, node_i);
+    }
+  }
+  MPI_Barrier(MPI_COMM_WORLD);
+}
+
+// There is some faster way of doing this
+// TODO: implement
+void apply_controlled_diagonal_gate(TangleState state, int control, int target,
+                                    TangleGate u) {
+  apply_controlled_gate(state, control, target, u);
+}
+
+// There is some faster way of doing this
+// TODO: implement
+void apply_controlled_antidiagonal_gate(TangleState state, int control,
+                                        int target, TangleGate u) {
+  apply_controlled_gate(state, control, target, u);
+}

src/gates.c

@@ -0,0 +1,53 @@
+#include <math.h>
+
+#include "tangle.h"
+
+void X(TangleState state, int target) {
+  TangleGate x;
+  x.A = 0;
+  x.B = 1;
+  x.C = 1;
+  x.D = 0;
+
+  apply_antidiagonal_gate(state, target, x);
+}
+
+void CX(TangleState state, int control, int target) {
+  TangleGate x;
+  x.A = 0;
+  x.B = 1;
+  x.C = 1;
+  x.D = 0;
+
+  apply_controlled_antidiagonal_gate(state, control, target, x);
+}
+
+void H(TangleState state, int target) {
+  TangleGate u;
+  u.A = 0.70710678118654752440;
+  u.B = 0.70710678118654752440;
+  u.C = 0.70710678118654752440;
+  u.D = -0.70710678118654752440;
+
+  apply_gate(state, target, u);
+}
+
+void Z(TangleState state, int target) {
+  TangleGate u;
+  u.A = 1;
+  u.B = 0;
+  u.C = 0;
+  u.D = -1;
+
+  apply_diagonal_gate(state, target, u);
+}
+
+void CZ(TangleState state, int control, int target) {
+  TangleGate u;
+  u.A = 1;
+  u.B = 0;
+  u.C = 0;
+  u.D = -1;
+
+  apply_controlled_diagonal_gate(state, control, target, u);
+}

src/main.c

@@ -0,0 +1,27 @@
+#include <stdio.h>
+#include "mpi.h"
+#include "tangle.h"
+
+int main() {
+  // Create the environment
+  TangleEnvironment env = initialize_environment();
+
+  int num_qubits = 24;
+  TangleState state = create_state(num_qubits, env);
+
+  for (int i = 0; i < num_qubits; i++) {
+    X(state, i);
+  }
+
+  llint result = measure(state);
+  
+
+  if (state.rank == 0) {
+    printf("Result: %lld\n", result);
+  }
+
+  // Finish the computation
+  destroy_environment(env);
+
+  return 0;
+}

src/measurment.c

@@ -0,0 +1,68 @@
+#include "mpi.h"
+#include "tangle.h"
+#include <complex.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+llint measure(TangleState state) {
+  float sample_value = 1.0;
+  llint measured = -1;
+
+  if (state.rank == 0) {
+    // The master node selects a random value between 0 and 1
+    srand(time(0));
+    sample_value = (float)rand() / (float)RAND_MAX;
+
+    // Then, the master node goes through it's own amplitudes, to
+    // attempt to select a measurment outcome
+    llint i;
+    for (i = 0; i < state.node_amps && sample_value > 0; i++) {
+      float amp_len = cabs(state.amps[i]);
+      sample_value = sample_value - (amp_len * amp_len);
+    }
+
+    // We now have to send the reduced random value to the other nodes
+    if (state.nodes > 1) {
+      MPI_Send(&sample_value, 1, MPI_REAL, 1, 5, MPI_COMM_WORLD);
+    }
+
+    // If the sample value is <= 0, then we have a measurment outcome.
+    // No other nodes will send a measurment value if sample value <= 0
+    if (sample_value <= 0) {
+      measured = state.rank * (1 << state.m) + (i - 1);
+    } else {
+      // Otherwise, we have to look on other nodes
+      MPI_Recv(&measured, 1, MPI_LONG_LONG_INT, MPI_ANY_SOURCE, 6,
+               MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    }
+  } else {
+    // This is not a master node, so we need to wait for the previous reduced
+    // random value
+    MPI_Recv(&sample_value, 1, MPI_REAL, state.rank - 1, 5, MPI_COMM_WORLD,
+             MPI_STATUS_IGNORE);
+
+    // Attempt to reduce the random value if needed
+    if (sample_value > 0) {
+      llint i;
+      for (i = 0; i < state.node_amps && sample_value > 0; i++) {
+        float amp_len = cabs(state.amps[i]);
+        sample_value = sample_value - (amp_len * amp_len);
+      }
+
+      // If the sample_value is now <= 0, then we have a measurment value, which
+      // can be communicated
+      if (sample_value <= 0) {
+        measured = state.rank * (1 << state.m) + (i - 1);
+        MPI_Send(&measured, 1, MPI_LONG_LONG_INT, 0, 6, MPI_COMM_WORLD);
+      }
+    }
+
+    // Again, the reduced random value can be passed on
+    if ((state.rank + 1) < state.nodes) {
+      MPI_Send(&sample_value, 1, MPI_REAL, state.rank + 1, 5, MPI_COMM_WORLD);
+    }
+  }
+
+  return measured;
+}

src/state.c

@@ -0,0 +1,48 @@
+#include "mpi.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "tangle.h"
+
+TangleState create_state(int num_qubits, TangleEnvironment env) {
+  TangleState state;
+
+  state.num_qubits = num_qubits;
+  state.rank = env.rank;
+  state.nodes = env.nodes;
+  state.k = log2(env.nodes);
+  state.m = num_qubits - state.k;
+
+  state.node_amps = 1LL << state.m;
+  state.temp_amps = state.node_amps / 2;
+
+  // Allocate the memory
+  llint num_amplitudes = state.node_amps + state.temp_amps;
+  state.amps = malloc(num_amplitudes * sizeof(cfloat));
+
+  for (int i = 0; i < num_amplitudes; i++) {
+    state.amps[i] = cfloat(0.0, 0.0);
+  }
+
+  if (state.rank == 0) {
+    state.amps[0] = cfloat(1.0, 0.0);
+  }
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
+  return state;
+}
+
+void print_state(TangleState state) {
+  for (int id = 0; id < state.nodes; id++) {
+    if (state.rank == id) {
+      for (llint i = 0; i < (1LL << state.m); i++) {
+        cfloat amp = state.amps[i];
+        printf("%lld: %f + i%f\n", state.rank * (1 << state.m) + i, creal(amp),
+               cimag(amp));
+      }
+    }
+    MPI_Barrier(MPI_COMM_WORLD);
+  }
+}