Merge pull request #459 from quantumlib/split-fusion-and-simulation

Split gate fusion and simulation.

Author: sergeisakov

docs/usage.md

@@ -22,12 +22,23 @@ Sample circuits are provided in
 |`-d maxtime` | maximum time |
 |`-t num_threads` | number of threads to use|
 |`-f max_fused_size` | maximum fused gate size|
-|`-v verbosity` | verbosity level (0,1,>1)|
+|`-v verbosity` | verbosity level (0,1,2,3,4,5)|
 |`-z` | set flush-to-zero and denormals-are-zeros MXCSR control flags|
 
 qsim_base computes all the amplitudes and just prints the first eight of them
 (or a smaller number for 1- or 2-qubit circuits).
 
+Verbosity levels are described in the following table.
+
+| Verbosity level | Description |
+|-----------------|-------------|
+| 0 | no additional information|
+| 1 | add total simulation runtime|
+| 2 | add initialization runtime and fuser runtime|
+| 3 | add basic fuser statistics|
+| 4 | add simulation runtime for each fused gate|
+| 5 | additional fuser information (qubit indices for each fused gate)|
+
 Example:
 ```
 ./qsim_base.x -c ../circuits/circuit_q24 -d 16 -t 8 -v 1
@@ -46,7 +57,7 @@ Example:
 |`-d maxtime` | maximum time |
 |`-t num_threads` | number of threads to use|
 |`-f max_fused_size` | maximum fused gate size|
-|`-v verbosity` | verbosity level (0,1,>1)|
+|`-v verbosity` | verbosity level (0,1,2,3,4,5)|
 |`-z` | set flush-to-zero and denormals-are-zeros MXCSR control flags|
 
 qsim_von_neumann computes all the amplitudes and calculates the von Neumann
@@ -77,7 +88,7 @@ Example:
 |`-o output_files` | comma-separated list of amplitude output files|
 |`-t num_threads` | number of threads to use|
 |`-f max_fused_size` | maximum fused gate size|
-|`-v verbosity` | verbosity level (0,1,>1)|
+|`-v verbosity` | verbosity level (0,1,2,3,4,5)|
 |`-z` | set flush-to-zero and denormals-are-zeros MXCSR control flags|
 
 qsim_amplitudes reads input files of bitstrings, computes the corresponding
@@ -112,7 +123,7 @@ Example:
 |`-t traj0` | starting trajectory |
 |`-n num_trajectories ` | number of trajectories to run starting with `traj0` |
 |`-f max_fused_size` | maximum fused gate size|
-|`-v verbosity` | verbosity level (0,1)|
+|`-v verbosity` | verbosity level (0,1,2,3,4,5)|
 
 qsim_qtrajectory_cuda runs on GPUs. qsim_qtrajectory_cuda performs quantum
 trajactory simulations with amplitude damping and phase damping noise channels.
@@ -145,7 +156,7 @@ Example:
 |`-p num_prefix_gates` | number of prefix gates|
 |`-r num_root_gates` | number of root gates|
 |`-t num_threads` | number of threads to use|
-|`-v verbosity` | verbosity level (0,>0)|
+|`-v verbosity` | verbosity level (0,1,4,5)|
 |`-z` | set flush-to-zero and denormals-are-zeros MXCSR control flags|
 
 qsimh_base just computes and just prints the first eight amplitudes. The hybrid
@@ -226,7 +237,7 @@ maximum "time".
 |`-i input_file` | bitstring input file|
 |`-o output_file` | amplitude output file|
 |`-t num_threads` | number of threads to use|
-|`-v verbosity` | verbosity level (0,>0)|
+|`-v verbosity` | verbosity level (0,1,4,5)|
 |`-z` | set flush-to-zero and denormals-are-zeros MXCSR control flags|
 
 qsimh_amplitudes reads the input file of bitstrings, computes the corresponding

lib/expect.h

@@ -123,8 +123,8 @@ std::complex<double> ExpectationValue(
         break;
       }
 
-      auto matrix = CalculateFusedMatrix<typename Simulator::fp_type>(fgate);
-      auto r = simulator.ExpectationValue(fgate.qubits, matrix.data(), state);
+      auto r = simulator.ExpectationValue(
+          fgate.qubits, fgate.matrix.data(), state);
       eval += str.weight * r;
     }
   }

lib/fuser.h

@@ -50,6 +50,10 @@ struct GateFused {
    * Ordered list of component gates.
    */
   std::vector<const Gate*> gates;
+  /**
+   * Fused gate matrix.
+   */
+  Matrix<typename Gate::fp_type> matrix;
 };
 
 /**
@@ -134,16 +138,14 @@ class Fuser {
 /**
  * Multiplies component gate matrices of a fused gate.
  * @param gate Fused gate.
- * @return Matrix product of component matrices.
  */
-template <typename fp_type, typename FusedGate>
-inline Matrix<fp_type> CalculateFusedMatrix(const FusedGate& gate) {
-  Matrix<fp_type> matrix;
-  MatrixIdentity(unsigned{1} << gate.qubits.size(), matrix);
+template <typename FusedGate>
+inline void CalculateFusedMatrix(FusedGate& gate) {
+  MatrixIdentity(unsigned{1} << gate.qubits.size(), gate.matrix);
 
   for (auto pgate : gate.gates) {
     if (gate.qubits.size() == pgate->qubits.size()) {
-      MatrixMultiply(gate.qubits.size(), pgate->matrix, matrix);
+      MatrixMultiply(gate.qubits.size(), pgate->matrix, gate.matrix);
     } else {
       unsigned mask = 0;
 
@@ -157,11 +159,31 @@ inline Matrix<fp_type> CalculateFusedMatrix(const FusedGate& gate) {
       }
 
       MatrixMultiply(mask, pgate->qubits.size(), pgate->matrix,
-                     gate.qubits.size(), matrix);
+                     gate.qubits.size(), gate.matrix);
     }
   }
+}
 
-  return matrix;
+/**
+ * Multiplies component gate matrices for a range of fused gates.
+ * @param gbeg, gend The iterator range [gbeg, gend) of fused gates.
+ */
+template <typename Iterator>
+inline void CalculateFusedMatrices(Iterator gbeg, Iterator gend) {
+  for (auto g = gbeg; g != gend; ++g) {
+    if (g->kind != gate::kMeasurement) {
+      CalculateFusedMatrix(*g);
+    }
+  }
+}
+
+/**
+ * Multiplies component gate matrices for a vector of fused gates.
+ * @param gates The vector of fused gates.
+ */
+template <typename FusedGate>
+inline void CalculateFusedMatrices(std::vector<FusedGate>& gates) {
+  CalculateFusedMatrices(gates.begin(), gates.end());
 }
 
 }  // namespace qsim

lib/fuser_basic.h

@@ -53,29 +53,29 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
 
   /**
    * Stores sets of gates that can be applied together. Only one- and
-   * two-qubit gates will get fused. Gates fused with this method are not
-   * multiplied together until ApplyFusedGate is called on the output.
-   * To respect specific time boundaries while fusing gates, use the other
-   * version of this method below.
+   * two-qubit gates will get fused. To respect specific time boundaries while
+   * fusing gates, use the other version of this method below.
    * @param param Options for gate fusion.
    * @param num_qubits The number of qubits acted on by 'gates'.
    * @param gates The gates (or pointers to the gates) to be fused.
    *   Gate times of the gates that act on the same qubits should be ordered.
    *   Gates that are out of time order should not cross the time boundaries
    *   set by measurement gates.
+   * @param fuse_matrix If true, multiply gate matrices together.
    * @return A vector of fused gate objects. Each element is a set of gates
    *   acting on a specific pair of qubits which can be applied as a group.
    */
   static std::vector<GateFused> FuseGates(const Parameter& param,
                                           unsigned num_qubits,
-                                          const std::vector<Gate>& gates) {
-    return FuseGates(param, num_qubits, gates.cbegin(), gates.cend(), {});
+                                          const std::vector<Gate>& gates,
+                                          bool fuse_matrix = true) {
+    return FuseGates(
+        param, num_qubits, gates.cbegin(), gates.cend(), {}, fuse_matrix);
   }
 
   /**
    * Stores sets of gates that can be applied together. Only one- and
-   * two-qubit gates will get fused. Gates fused with this method are not
-   * multiplied together until ApplyFusedGate is called on the output.
+   * two-qubit gates will get fused.
    * @param param Options for gate fusion.
    * @param num_qubits The number of qubits acted on by 'gates'.
    * @param gates The gates (or pointers to the gates) to be fused.
@@ -85,43 +85,44 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
    * @param times_to_split_at Ordered list of time steps (boundaries) at which
    *   to separate fused gates. Each element of the output will contain gates
    *   from a single 'window' in this list.
+   * @param fuse_matrix If true, multiply gate matrices together.
    * @return A vector of fused gate objects. Each element is a set of gates
    *   acting on a specific pair of qubits which can be applied as a group.
    */
   static std::vector<GateFused> FuseGates(
       const Parameter& param,
       unsigned num_qubits, const std::vector<Gate>& gates,
-      const std::vector<unsigned>& times_to_split_at) {
+      const std::vector<unsigned>& times_to_split_at,
+      bool fuse_matrix = true) {
     return FuseGates(param, num_qubits, gates.cbegin(), gates.cend(),
-                     times_to_split_at);
+                     times_to_split_at, fuse_matrix);
   }
 
   /**
    * Stores sets of gates that can be applied together. Only one- and
-   * two-qubit gates will get fused. Gates fused with this method are not
-   * multiplied together until ApplyFusedGate is called on the output.
-   * To respect specific time boundaries while fusing gates, use the other
-   * version of this method below.
+   * two-qubit gates will get fused. To respect specific time boundaries while
+   * fusing gates, use the other version of this method below.
    * @param param Options for gate fusion.
    * @param num_qubits The number of qubits acted on by gates.
    * @param gfirst, glast The iterator range [gfirst, glast) to fuse gates
    *   (or pointers to gates) in. Gate times of the gates that act on the same
    *   qubits should be ordered. Gates that are out of time order should not
    *   cross the time boundaries set by measurement gates.
+   * @param fuse_matrix If true, multiply gate matrices together.
    * @return A vector of fused gate objects. Each element is a set of gates
    *   acting on a specific pair of qubits which can be applied as a group.
    */
   static std::vector<GateFused> FuseGates(
       const Parameter& param, unsigned num_qubits,
       typename std::vector<Gate>::const_iterator gfirst,
-      typename std::vector<Gate>::const_iterator glast) {
-    return FuseGates(param, num_qubits, gfirst, glast, {});
+      typename std::vector<Gate>::const_iterator glast,
+      bool fuse_matrix = true) {
+    return FuseGates(param, num_qubits, gfirst, glast, {}, fuse_matrix);
   }
 
   /**
    * Stores sets of gates that can be applied together. Only one- and
-   * two-qubit gates will get fused. Gates fused with this method are not
-   * multiplied together until ApplyFusedGate is called on the output.
+   * two-qubit gates will get fused.
    * @param param Options for gate fusion.
    * @param num_qubits The number of qubits acted on by gates.
    * @param gfirst, glast The iterator range [gfirst, glast) to fuse gates
@@ -132,14 +133,16 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
    * @param times_to_split_at Ordered list of time steps (boundaries) at which
    *   to separate fused gates. Each element of the output will contain gates
    *   from a single 'window' in this list.
+   * @param fuse_matrix If true, multiply gate matrices together.
    * @return A vector of fused gate objects. Each element is a set of gates
    *   acting on a specific pair of qubits which can be applied as a group.
    */
   static std::vector<GateFused> FuseGates(
       const Parameter& param, unsigned num_qubits,
       typename std::vector<Gate>::const_iterator gfirst,
       typename std::vector<Gate>::const_iterator glast,
-      const std::vector<unsigned>& times_to_split_at) {
+      const std::vector<unsigned>& times_to_split_at,
+      bool fuse_matrix = true) {
     std::vector<GateFused> gates_fused;
 
     if (gfirst >= glast) return gates_fused;
@@ -243,11 +246,11 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
           }
 
           gates_fused.push_back({pgate->kind, pgate->time, pgate->qubits,
-                                 pgate, {pgate}});
+                                 pgate, {pgate}, {}});
         } else if (pgate->qubits.size() == 1) {
           unsigned q0 = pgate->qubits[0];
 
-          GateFused gate_f = {pgate->kind, pgate->time, {q0}, pgate, {}};
+          GateFused gate_f = {pgate->kind, pgate->time, {q0}, pgate, {}, {}};
 
           last[q0] = Advance(last[q0], gates_lat[q0], gate_f.gates);
           gate_f.gates.push_back(gates_lat[q0][last[q0]]);
@@ -260,7 +263,8 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
 
           if (Done(last[q0], pgate->time, gates_lat[q0])) continue;
 
-          GateFused gate_f = {pgate->kind, pgate->time, {q0, q1}, pgate, {}};
+          GateFused gate_f =
+              {pgate->kind, pgate->time, {q0, q1}, pgate, {}, {}};
 
           do {
             last[q0] = Advance(last[q0], gates_lat[q0], gate_f.gates);
@@ -290,7 +294,7 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
 
         const auto& mea_gates_at_time = measurement_gates[pgate->time];
 
-        GateFused gate_f = {pgate->kind, pgate->time, {}, pgate, {}};
+        GateFused gate_f = {pgate->kind, pgate->time, {}, pgate, {}, {}};
         gate_f.gates.reserve(mea_gates_at_time.size());
 
         // Fuse measurement gates with equal times.
@@ -307,6 +311,14 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
       if (gate_it == glast) break;
     }
 
+    if (fuse_matrix) {
+      for (auto& gate_f : gates_fused) {
+        if (gate_f.kind != gate::kMeasurement && gate_f.kind != gate::kDecomp) {
+          CalculateFusedMatrix(gate_f);
+        }
+      }
+    }
+
     return gates_fused;
   }
 
@@ -338,7 +350,7 @@ class BasicGateFuser final : public Fuser<IO, Gate> {
                                    std::vector<GateFused>& gates_fused) {
     auto pgate = gates_lat[q][k];
 
-    GateFused gate_f = {pgate->kind, pgate->time, {q}, pgate, {}};
+    GateFused gate_f = {pgate->kind, pgate->time, {q}, pgate, {}, {}};
     gate_f.gates.push_back(pgate);
 
     k = Advance(k + 1, gates_lat[q], gate_f.gates);