diff --git a/pyquil/quilbase.py b/pyquil/quilbase.py
index 896fa59d5..e54f6b16c 100644
--- a/pyquil/quilbase.py
+++ b/pyquil/quilbase.py
@@ -221,6 +221,16 @@ def __init__(
         self.qubits = qubits_list
         self.modifiers: List[str] = []
 
+    @property
+    def modified_name(self) -> str:
+        """ If there's a modifier on this gate then the 'official' name
+            of the gate needs to change, otherwise gate lookups for
+            matrices don't work.
+        """
+        if "CONTROLLED" in self.modifiers:
+            return "C" + self.name
+        return self.name
+
     def get_qubits(self, indices: bool = True) -> Set[QubitDesignator]:
         return {_extract_qubit_index(q, indices) for q in self.qubits}
 
diff --git a/pyquil/simulation/_numpy.py b/pyquil/simulation/_numpy.py
index b27525634..5914bf810 100644
--- a/pyquil/simulation/_numpy.py
+++ b/pyquil/simulation/_numpy.py
@@ -150,9 +150,9 @@ def _get_gate_tensor_and_qubits(gate: Gate) -> Tuple[np.ndarray, List[int]]:
     :return: tensor, qubit_inds.
     """
     if len(gate.params) > 0:
-        matrix = QUANTUM_GATES[gate.name](*gate.params)
+        matrix = QUANTUM_GATES[gate.modified_name](*gate.params)
     else:
-        matrix = QUANTUM_GATES[gate.name]
+        matrix = QUANTUM_GATES[gate.modified_name]
 
     qubit_inds = [q.index for q in gate.qubits]
 
diff --git a/pyquil/simulation/matrices.py b/pyquil/simulation/matrices.py
index 1a55fc7a8..12b095e66 100644
--- a/pyquil/simulation/matrices.py
+++ b/pyquil/simulation/matrices.py
@@ -107,6 +107,8 @@
 
 import numpy as np
 
+INV_ROOT2 = 1.0 / np.sqrt(2.0)
+
 I = np.array([[1.0, 0.0], [0.0, 1.0]])
 
 X = np.array([[0.0, 1.0], [1.0, 0.0]])
@@ -115,7 +117,16 @@
 
 Z = np.array([[1.0, 0.0], [0.0, -1.0]])
 
-H = (1.0 / np.sqrt(2.0)) * np.array([[1.0, 1.0], [1.0, -1.0]])
+H = INV_ROOT2 * np.array([[1.0, 1.0], [1.0, -1.0]])
+
+CH = np.array(
+    [
+        [1.0, 0.0, 0.0, 0.0],
+        [0.0, 1.0, 0.0, 0.0],
+        [0.0, 0.0, INV_ROOT2, INV_ROOT2],
+        [0.0, 0.0, INV_ROOT2, -INV_ROOT2],
+    ]
+)
 
 S = np.array([[1.0, 0.0], [0.0, 1.0j]])
 
@@ -236,9 +247,11 @@ def BARENCO(alpha: float, phi: float, theta: float) -> np.ndarray:
 QUANTUM_GATES = {
     "I": I,
     "X": X,
+    "CX": CNOT,
     "Y": Y,
     "Z": Z,
     "H": H,
+    "CH": CH,
     "S": S,
     "T": T,
     "PHASE": PHASE,
diff --git a/pyquil/tests/test_pyqvm.py b/pyquil/tests/test_pyqvm.py
new file mode 100644
index 000000000..e08c1d592
--- /dev/null
+++ b/pyquil/tests/test_pyqvm.py
@@ -0,0 +1,30 @@
+import pytest
+
+from pyquil import Program
+from pyquil.gates import H, X, CNOT, Z, SWAP
+from pyquil.pyqvm import PyQVM
+
+
+@pytest.mark.parametrize(
+    "first_gate,second_gate,expected",
+    [
+        (X(0), H(1), "X 0\nH 1\n"),
+        (X(0), H(1).controlled(0), "X 0\nCONTROLLED H 0 1\n"),
+        (H(0), X(1), "H 0\nX 1\n"),
+        (H(0), X(1).controlled(0), "H 0\nCONTROLLED X 0 1\n"),
+        (H(0), Z(1), "H 0\nZ 1\n"),
+        (H(0), Z(1).controlled(0), "H 0\nCONTROLLED Z 0 1\n"),
+        (X(0), X(1), "X 0\nX 1\n"),
+        (X(0), X(1).controlled(0), "X 0\nCONTROLLED X 0 1\n"),
+        (H(0), SWAP(0, 1), "H 0\nSWAP 0 1\n"),
+        (X(0), CNOT(0, 1), "X 0\nCNOT 0 1\n"),
+    ],
+)
+def test_pyqvm_controlled_gates(first_gate, second_gate, expected):
+    """ Unit-test based on the bug report in
+        https://github.com/rigetti/pyquil/issues/1259
+    """
+    p = Program(first_gate, second_gate)
+    qvm = PyQVM(n_qubits=2)
+    result = qvm.execute(p)
+    assert result.program.out() == expected