diff --git a/qiskit/circuit/add_control.py b/qiskit/circuit/add_control.py
index 3849136fbffe..43e2b8e7c816 100644
--- a/qiskit/circuit/add_control.py
+++ b/qiskit/circuit/add_control.py
@@ -99,57 +99,75 @@ def control(operation: Union[Gate, ControlledGate],
     from math import pi
     # pylint: disable=cyclic-import
     import qiskit.circuit.controlledgate as controlledgate
-    # pylint: disable=unused-import
-    import qiskit.circuit.library.standard_gates.multi_control_rotation_gates
 
     q_control = QuantumRegister(num_ctrl_qubits, name='control')
     q_target = QuantumRegister(operation.num_qubits, name='target')
     q_ancillae = None  # TODO: add
-    qc = QuantumCircuit(q_control, q_target)
-
+    controlled_circ = QuantumCircuit(q_control, q_target,
+                                     name='c_{}'.format(operation.name))
+    global_phase = 0
     if operation.name == 'x' or (
             isinstance(operation, controlledgate.ControlledGate) and
             operation.base_gate.name == 'x'):
-        qc.mct(q_control[:] + q_target[:-1], q_target[-1], q_ancillae)
-    elif operation.name == 'rx':
-        qc.mcrx(operation.definition.data[0][0].params[0], q_control, q_target[0],
-                use_basis_gates=True)
-    elif operation.name == 'ry':
-        qc.mcry(operation.definition.data[0][0].params[0], q_control, q_target[0],
-                q_ancillae, mode='noancilla', use_basis_gates=True)
-    elif operation.name == 'rz':
-        qc.mcrz(operation.definition.data[0][0].params[0], q_control, q_target[0],
-                use_basis_gates=True)
+        controlled_circ.mct(q_control[:] + q_target[:-1], q_target[-1], q_ancillae)
+        if operation.definition is not None and operation.definition.global_phase:
+            global_phase += operation.definition.global_phase
     else:
-        bgate = _unroll_gate(operation, ['u1', 'u3', 'cx'])
-        # now we have a bunch of single qubit rotation gates and cx
-        for rule in bgate.definition.data:
-            if rule[0].name == 'u3':
-                theta, phi, lamb = rule[0].params
+        basis = ['u1', 'u3', 'x', 'rx', 'ry', 'rz', 'cx']
+        unrolled_gate = _unroll_gate(operation, basis_gates=basis)
+        for gate, qreg, _ in unrolled_gate.definition.data:
+            if gate.name == 'x':
+                controlled_circ.mct(q_control, q_target[qreg[0].index],
+                                    q_ancillae)
+            elif gate.name == 'rx':
+                controlled_circ.mcrx(gate.definition.data[0][0].params[0],
+                                     q_control, q_target[qreg[0].index],
+                                     use_basis_gates=True)
+            elif gate.name == 'ry':
+                controlled_circ.mcry(gate.definition.data[0][0].params[0],
+                                     q_control, q_target[qreg[0].index],
+                                     q_ancillae, mode='noancilla',
+                                     use_basis_gates=True)
+            elif gate.name == 'rz':
+                controlled_circ.mcrz(gate.definition.data[0][0].params[0],
+                                     q_control, q_target[qreg[0].index],
+                                     use_basis_gates=True)
+            elif gate.name == 'u1':
+                controlled_circ.mcu1(gate.params[0], q_control, q_target[qreg[0].index])
+            elif gate.name == 'cx':
+                controlled_circ.mct(q_control[:] + [q_target[qreg[0].index]],
+                                    q_target[qreg[1].index],
+                                    q_ancillae)
+            elif gate.name == 'u3':
+                theta, phi, lamb = gate.params
                 if phi == -pi / 2 and lamb == pi / 2:
-                    qc.mcrx(theta, q_control, q_target[rule[1][0].index],
-                            use_basis_gates=True)
+                    controlled_circ.mcrx(theta, q_control, q_target[qreg[0].index],
+                                         use_basis_gates=True)
                 elif phi == 0 and lamb == 0:
-                    qc.mcry(theta, q_control, q_target[rule[1][0].index],
-                            q_ancillae, use_basis_gates=True)
+                    controlled_circ.mcry(theta, q_control, q_target[qreg[0].index],
+                                         q_ancillae, use_basis_gates=True)
                 elif theta == 0 and phi == 0:
-                    qc.mcrz(lamb, q_control, q_target[rule[1][0].index],
-                            use_basis_gates=True)
+                    controlled_circ.mcrz(lamb, q_control, q_target[qreg[0].index],
+                                         use_basis_gates=True)
                 else:
-                    qc.mcrz(lamb, q_control, q_target[rule[1][0].index],
-                            use_basis_gates=True)
-                    qc.mcry(theta, q_control, q_target[rule[1][0].index],
-                            q_ancillae, use_basis_gates=True)
-                    qc.mcrz(phi, q_control, q_target[rule[1][0].index],
-                            use_basis_gates=True)
-            elif rule[0].name == 'u1':
-                qc.mcu1(rule[0].params[0], q_control, q_target[rule[1][0].index])
-            elif rule[0].name == 'cx':
-                qc.mct(q_control[:] + [q_target[rule[1][0].index]], q_target[rule[1][1].index],
-                       q_ancillae)
+                    controlled_circ.mcrz(lamb, q_control, q_target[qreg[0].index],
+                                         use_basis_gates=True)
+                    controlled_circ.mcry(theta, q_control, q_target[qreg[0].index],
+                                         q_ancillae, use_basis_gates=True)
+                    controlled_circ.mcrz(phi, q_control, q_target[qreg[0].index],
+                                         use_basis_gates=True)
             else:
-                raise CircuitError('gate contains non-controllable instructions')
-
+                raise CircuitError('gate contains non-controllable instructions: {}'.format(
+                    gate.name))
+            if gate.definition is not None and gate.definition.global_phase:
+                global_phase += gate.definition.global_phase
+    # apply controlled global phase
+    if ((operation.definition is not None and operation.definition.global_phase) or global_phase):
+        if len(q_control) < 2:
+            controlled_circ.u1(operation.definition.global_phase + global_phase, q_control)
+        else:
+            controlled_circ.mcu1(operation.definition.global_phase + global_phase,
+                                 q_control[:-1], q_control[-1])
     if isinstance(operation, controlledgate.ControlledGate):
         new_num_ctrl_qubits = num_ctrl_qubits + operation.num_ctrl_qubits
         new_ctrl_state = operation.ctrl_state << num_ctrl_qubits | ctrl_state
@@ -170,11 +188,11 @@ def control(operation: Union[Gate, ControlledGate],
         ctrl_substr = ('{0}' * new_num_ctrl_qubits).format('c')
     new_name = '{0}{1}'.format(ctrl_substr, base_name)
     cgate = controlledgate.ControlledGate(new_name,
-                                          qc.num_qubits,
+                                          controlled_circ.num_qubits,
                                           operation.params,
                                           label=label,
                                           num_ctrl_qubits=new_num_ctrl_qubits,
-                                          definition=qc,
+                                          definition=controlled_circ,
                                           ctrl_state=new_ctrl_state)
     cgate.base_gate = base_gate
     return cgate
@@ -196,7 +214,7 @@ def _gate_to_circuit(operation):
 
 def _gate_to_dag(operation):
     from qiskit.converters.circuit_to_dag import circuit_to_dag
-    if hasattr(operation, 'definition') and operation.definition:
+    if hasattr(operation, 'definition') and operation.definition is not None:
         return circuit_to_dag(operation.definition)
     else:
         qr = QuantumRegister(operation.num_qubits)
@@ -210,5 +228,5 @@ def _unroll_gate(operation, basis_gates):
     from qiskit.transpiler.passes import Unroller
     unroller = Unroller(basis_gates)
     dag = _gate_to_dag(operation)
-    qc = dag_to_circuit(unroller.run(dag))
-    return qc.to_gate()
+    opqc = dag_to_circuit(unroller.run(dag))
+    return opqc.to_gate()
diff --git a/qiskit/circuit/controlledgate.py b/qiskit/circuit/controlledgate.py
index 5353108616d6..24a466e3d54c 100644
--- a/qiskit/circuit/controlledgate.py
+++ b/qiskit/circuit/controlledgate.py
@@ -207,4 +207,5 @@ def __eq__(self, other) -> bool:
 
     def inverse(self) -> 'ControlledGate':
         """Invert this gate by calling inverse on the base gate."""
-        return self.base_gate.inverse().control(self.num_ctrl_qubits)
+        return self.base_gate.inverse().control(self.num_ctrl_qubits,
+                                                ctrl_state=self.ctrl_state)
diff --git a/qiskit/circuit/equivalence.py b/qiskit/circuit/equivalence.py
index b1701714e5b4..d4e905e085bd 100644
--- a/qiskit/circuit/equivalence.py
+++ b/qiskit/circuit/equivalence.py
@@ -258,7 +258,6 @@ def _raise_if_shape_mismatch(gate, circuit):
 
 def _rebind_equiv(equiv, query_params):
     equiv_params, equiv_circuit = equiv
-
     param_map = dict(zip(equiv_params, query_params))
     equiv = equiv_circuit.assign_parameters(param_map, inplace=False)
 
diff --git a/qiskit/circuit/library/standard_gates/equivalence_library.py b/qiskit/circuit/library/standard_gates/equivalence_library.py
index c019fa3f3e2d..fabc105ea258 100644
--- a/qiskit/circuit/library/standard_gates/equivalence_library.py
+++ b/qiskit/circuit/library/standard_gates/equivalence_library.py
@@ -36,6 +36,7 @@
     RZGate,
     CRZGate,
     RZZGate,
+    RZXGate,
     SGate,
     SdgGate,
     SwapGate,
@@ -158,7 +159,7 @@
         (HGate(), [q[0]], []),
         (HGate(), [q[1]], []),
         (CXGate(), [q[0], q[1]], []),
-        (U1Gate(theta), [q[1]], []),
+        (RZGate(theta), [q[1]], []),
         (CXGate(), [q[0], q[1]], []),
         (HGate(), [q[1]], []),
         (HGate(), [q[0]], []),
@@ -166,6 +167,22 @@
     def_rxx.append(inst, qargs, cargs)
 _sel.add_equivalence(RXXGate(theta), def_rxx)
 
+# RZXGate
+
+q = QuantumRegister(2, 'q')
+theta = Parameter('theta')
+def_rzx = QuantumCircuit(q)
+for inst, qargs, cargs in [
+        (HGate(), [q[1]], []),
+        (CXGate(), [q[0], q[1]], []),
+        (RZGate(theta), [q[1]], []),
+        (CXGate(), [q[0], q[1]], []),
+        (HGate(), [q[1]], []),
+]:
+    def_rzx.append(inst, qargs, cargs)
+_sel.add_equivalence(RZXGate(theta), def_rzx)
+
+
 # RYGate
 
 q = QuantumRegister(1, 'q')
@@ -209,7 +226,7 @@
 
 q = QuantumRegister(1, 'q')
 theta = Parameter('theta')
-def_rz = QuantumCircuit(q)
+def_rz = QuantumCircuit(q, global_phase=-theta / 2)
 def_rz.append(U1Gate(theta), [q[0]], [])
 _sel.add_equivalence(RZGate(theta), def_rz)
 
@@ -242,7 +259,7 @@
 def_rzz = QuantumCircuit(q)
 for inst, qargs, cargs in [
         (CXGate(), [q[0], q[1]], []),
-        (U1Gate(theta), [q[1]], []),
+        (RZGate(theta), [q[1]], []),
         (CXGate(), [q[0], q[1]], [])
 ]:
     def_rzz.append(inst, qargs, cargs)
@@ -338,7 +355,7 @@
 def_tdg.append(U1Gate(-pi / 4), [q[0]], [])
 _sel.add_equivalence(TdgGate(), def_tdg)
 
-# U1Gate
+# U2Gate
 
 q = QuantumRegister(1, 'q')
 phi = Parameter('phi')
@@ -376,7 +393,7 @@
 theta = Parameter('theta')
 phi = Parameter('phi')
 lam = Parameter('lam')
-u3_qasm_def = QuantumCircuit(q)
+u3_qasm_def = QuantumCircuit(q, global_phase=(lam + phi) / 2)
 u3_qasm_def.rz(lam, 0)
 u3_qasm_def.rx(pi/2, 0)
 u3_qasm_def.rz(theta+pi, 0)
@@ -427,7 +444,7 @@
 _sel.add_equivalence(CXGate(), cx_to_cz)
 
 q = QuantumRegister(2, 'q')
-cx_to_iswap = QuantumCircuit(q)
+cx_to_iswap = QuantumCircuit(q, global_phase=3*pi/4)
 for inst, qargs, cargs in [
         (HGate(), [q[0]], []),
         (XGate(), [q[1]], []),
diff --git a/qiskit/circuit/library/standard_gates/ms.py b/qiskit/circuit/library/standard_gates/ms.py
index efae92180708..47ca3c99ee48 100644
--- a/qiskit/circuit/library/standard_gates/ms.py
+++ b/qiskit/circuit/library/standard_gates/ms.py
@@ -39,11 +39,12 @@ def _define(self):
         # pylint: disable=cyclic-import
         from qiskit.circuit.quantumcircuit import QuantumCircuit
         from .rxx import RXXGate
+        theta = self.params[0]
         q = QuantumRegister(self.num_qubits, 'q')
         qc = QuantumCircuit(q, name=self.name)
         rules = []
         for i in range(self.num_qubits):
             for j in range(i + 1, self.num_qubits):
-                rules += [(RXXGate(self.params[0]), [q[i], q[j]], [])]
+                rules += [(RXXGate(theta), [q[i], q[j]], [])]
         qc._data = rules
         self.definition = qc
diff --git a/qiskit/circuit/library/standard_gates/rx.py b/qiskit/circuit/library/standard_gates/rx.py
index 2c996d7cd7d9..76d178ce809d 100644
--- a/qiskit/circuit/library/standard_gates/rx.py
+++ b/qiskit/circuit/library/standard_gates/rx.py
@@ -200,18 +200,23 @@ def inverse(self):
         """Return inverse RX gate (i.e. with the negative rotation angle)."""
         return CRXGate(-self.params[0])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the CRX gate."""
-    #    half_theta = self.params[0] / 2
-    #    cos = numpy.cos(half_theta)
-    #    isin = 1j * numpy.sin(half_theta)
-    #    return numpy.array([[1,     0, 0,     0],
-    #                        [0,   cos, 0, -isin],
-    #                        [0,     0, 1,     0],
-    #                        [0, -isin, 0,   cos]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the CRX gate."""
+        half_theta = self.params[0] / 2
+        cos = numpy.cos(half_theta)
+        isin = 1j * numpy.sin(half_theta)
+        if self.ctrl_state:
+            return numpy.array([[1, 0, 0, 0],
+                                [0, cos, 0, -isin],
+                                [0, 0, 1, 0],
+                                [0, -isin, 0, cos]],
+                               dtype=complex)
+        else:
+            return numpy.array([[cos, 0, -isin, 0],
+                                [0, 1, 0, 0],
+                                [-isin, 0, cos, 0],
+                                [0, 0, 0, 1]],
+                               dtype=complex)
 
 
 class CrxGate(CRXGate, metaclass=CRXMeta):
diff --git a/qiskit/circuit/library/standard_gates/rxx.py b/qiskit/circuit/library/standard_gates/rxx.py
index f08fa34deaac..5c4539c79f4b 100644
--- a/qiskit/circuit/library/standard_gates/rxx.py
+++ b/qiskit/circuit/library/standard_gates/rxx.py
@@ -77,16 +77,16 @@ def _define(self):
         # pylint: disable=cyclic-import
         from qiskit.circuit.quantumcircuit import QuantumCircuit
         from .x import CXGate
-        from .u1 import U1Gate
         from .h import HGate
+        from .rz import RZGate
+        theta = self.params[0]
         q = QuantumRegister(2, 'q')
         qc = QuantumCircuit(q, name=self.name)
-        theta = self.params[0]
         rules = [
             (HGate(), [q[0]], []),
             (HGate(), [q[1]], []),
             (CXGate(), [q[0], q[1]], []),
-            (U1Gate(theta), [q[1]], []),
+            (RZGate(theta), [q[1]], []),
             (CXGate(), [q[0], q[1]], []),
             (HGate(), [q[1]], []),
             (HGate(), [q[0]], []),
@@ -98,14 +98,14 @@ def inverse(self):
         """Return inverse RXX gate (i.e. with the negative rotation angle)."""
         return RXXGate(-self.params[0])
 
-    # NOTE: we should use the following as the canonical matrix
-    # definition but we don't include it yet since it differs from
-    # the circuit decomposition matrix by a global phase
-    # def to_matrix(self):
-    #   """Return a Numpy.array for the RXX gate."""
-    #    theta = float(self.params[0])
-    #    return np.array([
-    #        [np.cos(theta / 2), 0, 0, -1j * np.sin(theta / 2)],
-    #        [0, np.cos(theta / 2), -1j * np.sin(theta / 2), 0],
-    #        [0, -1j * np.sin(theta / 2), np.cos(theta / 2), 0],
-    #        [-1j * np.sin(theta / 2), 0, 0, np.cos(theta / 2)]], dtype=complex)
+    def to_matrix(self):
+        """Return a Numpy.array for the RXX gate."""
+        import numpy
+        theta2 = float(self.params[0]) / 2
+        cos = numpy.cos(theta2)
+        isin = 1j * numpy.sin(theta2)
+        return numpy.array([
+            [cos, 0, 0, -isin],
+            [0, cos, -isin, 0],
+            [0, -isin, cos, 0],
+            [-isin, 0, 0, cos]], dtype=complex)
diff --git a/qiskit/circuit/library/standard_gates/ry.py b/qiskit/circuit/library/standard_gates/ry.py
index 5a85ce740074..d0b1b372f410 100644
--- a/qiskit/circuit/library/standard_gates/ry.py
+++ b/qiskit/circuit/library/standard_gates/ry.py
@@ -195,18 +195,23 @@ def inverse(self):
         """Return inverse RY gate (i.e. with the negative rotation angle)."""
         return CRYGate(-self.params[0])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the CRY gate."""
-    #    half_theta = self.params[0] / 2
-    #    cos = numpy.cos(half_theta)
-    #    sin = numpy.sin(half_theta)
-    #    return numpy.array([[1,     0, 0,    0],
-    #                        [0,   cos, 0, -sin],
-    #                        [0,     0, 1,    0],
-    #                        [0,   sin, 0,  cos]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the CRY gate."""
+        half_theta = self.params[0] / 2
+        cos = numpy.cos(half_theta)
+        sin = numpy.sin(half_theta)
+        if self.ctrl_state:
+            return numpy.array([[1, 0, 0, 0],
+                                [0, cos, 0, -sin],
+                                [0, 0, 1, 0],
+                                [0, sin, 0, cos]],
+                               dtype=complex)
+        else:
+            return numpy.array([[cos, 0, -sin, 0],
+                                [0, 1, 0, 0],
+                                [sin, 0, cos, 0],
+                                [0, 0, 0, 1]],
+                               dtype=complex)
 
 
 class CryGate(CRYGate, metaclass=CRYMeta):
diff --git a/qiskit/circuit/library/standard_gates/ryy.py b/qiskit/circuit/library/standard_gates/ryy.py
index 04869f7117eb..842995570c26 100644
--- a/qiskit/circuit/library/standard_gates/ryy.py
+++ b/qiskit/circuit/library/standard_gates/ryy.py
@@ -82,8 +82,8 @@ def _define(self):
         from .rz import RZGate
 
         q = QuantumRegister(2, 'q')
-        qc = QuantumCircuit(q, name=self.name)
         theta = self.params[0]
+        qc = QuantumCircuit(q, name=self.name)
         rules = [
             (RXGate(np.pi / 2), [q[0]], []),
             (RXGate(np.pi / 2), [q[1]], []),
@@ -100,14 +100,14 @@ def inverse(self):
         """Return inverse RYY gate (i.e. with the negative rotation angle)."""
         return RYYGate(-self.params[0])
 
-    # TODO: this is the correct matrix and is equal to the definition above,
-    # however the control mechanism cannot distinguish U1 and RZ yet.
-    # def to_matrix(self):
-    #     """Return a numpy.array for the RYY gate."""
-    #     theta = self.params[0]
-    #     return np.exp(0.5j * theta) * np.array([
-    #         [np.cos(theta / 2), 0, 0, 1j * np.sin(theta / 2)],
-    #         [0, np.cos(theta / 2), -1j * np.sin(theta / 2), 0],
-    #         [0, -1j * np.sin(theta / 2), np.cos(theta / 2), 0],
-    #         [1j * np.sin(theta / 2), 0, 0, np.cos(theta / 2)]
-    #     ], dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the RYY gate."""
+        theta = self.params[0]
+        cos = np.cos(theta / 2)
+        isin = 1j * np.sin(theta / 2)
+        return np.array([
+            [cos, 0, 0, isin],
+            [0, cos, -isin, 0],
+            [0, -isin, cos, 0],
+            [isin, 0, 0, cos]
+        ], dtype=complex)
diff --git a/qiskit/circuit/library/standard_gates/rz.py b/qiskit/circuit/library/standard_gates/rz.py
index 03942c4758bd..1e9837652a67 100644
--- a/qiskit/circuit/library/standard_gates/rz.py
+++ b/qiskit/circuit/library/standard_gates/rz.py
@@ -68,9 +68,10 @@ def _define(self):
         from qiskit.circuit.quantumcircuit import QuantumCircuit
         from .u1 import U1Gate
         q = QuantumRegister(1, 'q')
-        qc = QuantumCircuit(q, name=self.name)
+        theta = self.params[0]
+        qc = QuantumCircuit(q, name=self.name, global_phase=-theta / 2)
         rules = [
-            (U1Gate(self.params[0]), [q[0]], [])
+            (U1Gate(theta), [q[0]], [])
         ]
         qc._data = rules
         self.definition = qc
@@ -100,14 +101,12 @@ def inverse(self):
         """
         return RZGate(-self.params[0])
 
-    # TODO: this is the correct matrix however the control mechanism
-    # cannot distinguish U1 and RZ yet.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the RZ gate."""
-    #    import numpy
-    #    lam = float(self.params[0])
-    #    return numpy.array([[numpy.exp(-1j * lam / 2), 0],
-    #                        [0, numpy.exp(1j * lam / 2)]], dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the RZ gate."""
+        import numpy as np
+        ilam2 = 0.5j * float(self.params[0])
+        return np.array([[np.exp(-ilam2), 0],
+                         [0, np.exp(ilam2)]], dtype=complex)
 
 
 class CRZMeta(type):
@@ -213,16 +212,22 @@ def inverse(self):
         """Return inverse RZ gate (i.e. with the negative rotation angle)."""
         return CRZGate(-self.params[0])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the CRZ gate."""
-    #    arg = 1j * self.params[0] / 2
-    #    return numpy.array([[1,               0, 0,              0],
-    #                        [0, numpy.exp(-arg), 0,              0],
-    #                        [0,               0, 1,              0],
-    #                        [0,               0, 0, numpy.exp(arg)]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the CRZ gate."""
+        import numpy
+        arg = 1j * self.params[0] / 2
+        if self.ctrl_state:
+            return numpy.array([[1, 0, 0, 0],
+                                [0, numpy.exp(-arg), 0, 0],
+                                [0, 0, 1, 0],
+                                [0, 0, 0, numpy.exp(arg)]],
+                               dtype=complex)
+        else:
+            return numpy.array([[numpy.exp(-arg), 0, 0, 0],
+                                [0, 1, 0, 0],
+                                [0, 0, numpy.exp(arg), 0],
+                                [0, 0, 0, 1]],
+                               dtype=complex)
 
 
 class CrzGate(CRZGate, metaclass=CRZMeta):
diff --git a/qiskit/circuit/library/standard_gates/rzx.py b/qiskit/circuit/library/standard_gates/rzx.py
index 037b437df047..fbb32dd22068 100644
--- a/qiskit/circuit/library/standard_gates/rzx.py
+++ b/qiskit/circuit/library/standard_gates/rzx.py
@@ -16,9 +16,6 @@
 
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.quantumregister import QuantumRegister
-from .rz import RZGate
-from .h import HGate
-from .x import CXGate
 
 
 class RZXGate(Gate):
@@ -127,12 +124,16 @@ def _define(self):
         """
         # pylint: disable=cyclic-import
         from qiskit.circuit.quantumcircuit import QuantumCircuit
+        from .h import HGate
+        from .x import CXGate
+        from .rz import RZGate
+        theta = self.params[0]
         q = QuantumRegister(2, 'q')
         qc = QuantumCircuit(q, name=self.name)
         rules = [
             (HGate(), [q[1]], []),
             (CXGate(), [q[0], q[1]], []),
-            (RZGate(self.params[0]), [q[1]], []),
+            (RZGate(theta), [q[1]], []),
             (CXGate(), [q[0], q[1]], []),
             (HGate(), [q[1]], [])
         ]
@@ -143,15 +144,14 @@ def inverse(self):
         """Return inverse RZX gate (i.e. with the negative rotation angle)."""
         return RZXGate(-self.params[0])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the RZX gate."""
-    #    half_theta = self.params[0] / 2
-    #    cos = numpy.cos(half_theta)
-    #    isin = 1j * numpy.sin(half_theta)
-    #    return numpy.array([[    cos,      0,   -isin,      0],
-    #                        [      0,    cos,       0,   isin],
-    #                        [-1j*sin,      0,     cos,      0],
-    #                        [      0,   isin,       0,    cos]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the RZX gate."""
+        import numpy
+        half_theta = self.params[0] / 2
+        cos = numpy.cos(half_theta)
+        isin = 1j * numpy.sin(half_theta)
+        return numpy.array([[cos, 0, -isin, 0],
+                            [0, cos, 0, isin],
+                            [-isin, 0, cos, 0],
+                            [0, isin, 0, cos]],
+                           dtype=complex)
diff --git a/qiskit/circuit/library/standard_gates/rzz.py b/qiskit/circuit/library/standard_gates/rzz.py
index 57611d1bd212..5f98aea69d17 100644
--- a/qiskit/circuit/library/standard_gates/rzz.py
+++ b/qiskit/circuit/library/standard_gates/rzz.py
@@ -91,13 +91,14 @@ def _define(self):
         """
         # pylint: disable=cyclic-import
         from qiskit.circuit.quantumcircuit import QuantumCircuit
-        from .u1 import U1Gate
         from .x import CXGate
+        from .rz import RZGate
         q = QuantumRegister(2, 'q')
+        theta = self.params[0]
         qc = QuantumCircuit(q, name=self.name)
         rules = [
             (CXGate(), [q[0], q[1]], []),
-            (U1Gate(self.params[0]), [q[1]], []),
+            (RZGate(theta), [q[1]], []),
             (CXGate(), [q[0], q[1]], [])
         ]
         qc._data = rules
@@ -107,13 +108,11 @@ def inverse(self):
         """Return inverse RZZ gate (i.e. with the negative rotation angle)."""
         return RZZGate(-self.params[0])
 
-    # TODO: this is the correct matrix and is equal to the definition above,
-    # however the control mechanism cannot distinguish U1 and RZ yet.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the RZZ gate."""
-    #    import numpy
-    #    theta = float(self.params[0])
-    #    return numpy.array([[numpy.exp(-1j*theta/2), 0, 0, 0],
-    #                     [0, numpy.exp(1j*theta/2), 0, 0],
-    #                     [0, 0, numpy.exp(1j*theta/2), 0],
-    #                     [0, 0, 0, numpy.exp(-1j*theta/2)]], dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the RZZ gate."""
+        import numpy
+        itheta2 = 1j * float(self.params[0]) / 2
+        return numpy.array([[numpy.exp(-itheta2), 0, 0, 0],
+                            [0, numpy.exp(itheta2), 0, 0],
+                            [0, 0, numpy.exp(itheta2), 0],
+                            [0, 0, 0, numpy.exp(-itheta2)]], dtype=complex)
diff --git a/qiskit/circuit/library/standard_gates/u1.py b/qiskit/circuit/library/standard_gates/u1.py
index 4d77e21a2a43..c45554150295 100644
--- a/qiskit/circuit/library/standard_gates/u1.py
+++ b/qiskit/circuit/library/standard_gates/u1.py
@@ -221,16 +221,22 @@ def inverse(self):
         r"""Return inverted CU1 gate (:math:`CU1(\lambda){\dagger} = CU1(-\lambda)`)"""
         return CU1Gate(-self.params[0])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the CU1 gate."""
-    #    eith = numpy.exp(1j * self.params[0])
-    #    return numpy.array([[1, 0, 0,    0],
-    #                        [0, 1, 0,    0],
-    #                        [0, 0, 1,    0],
-    #                        [0, 0, 0, eith]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the CU1 gate."""
+
+        eith = numpy.exp(1j * float(self.params[0]))
+        if self.ctrl_state:
+            return numpy.array([[1, 0, 0, 0],
+                                [0, 1, 0, 0],
+                                [0, 0, 1, 0],
+                                [0, 0, 0, eith]],
+                               dtype=complex)
+        else:
+            return numpy.array([[1, 0, 0, 0],
+                                [0, 1, 0, 0],
+                                [0, 0, eith, 0],
+                                [0, 0, 0, 1]],
+                               dtype=complex)
 
 
 class Cu1Gate(CU1Gate, metaclass=CU1Meta):
diff --git a/qiskit/circuit/library/standard_gates/u3.py b/qiskit/circuit/library/standard_gates/u3.py
index 913cb008a4a2..5d11fa528d90 100644
--- a/qiskit/circuit/library/standard_gates/u3.py
+++ b/qiskit/circuit/library/standard_gates/u3.py
@@ -93,15 +93,11 @@ def to_matrix(self):
         """Return a Numpy.array for the U3 gate."""
         theta, phi, lam = self.params
         theta, phi, lam = float(theta), float(phi), float(lam)
+        cos = numpy.cos(theta / 2)
+        sin = numpy.sin(theta / 2)
         return numpy.array([
-            [
-                numpy.cos(theta / 2),
-                -numpy.exp(1j * lam) * numpy.sin(theta / 2)
-            ],
-            [
-                numpy.exp(1j * phi) * numpy.sin(theta / 2),
-                numpy.exp(1j * (phi + lam)) * numpy.cos(theta / 2)
-            ]
+            [cos, -numpy.exp(1j * lam) * sin],
+            [numpy.exp(1j * phi) * sin, numpy.exp(1j * (phi + lam)) * cos]
         ], dtype=complex)
 
 
@@ -142,7 +138,7 @@ class CU3Gate(ControlledGate, metaclass=CU3Meta):
             U3(\theta,\phi,\lambda) \otimes |1\rangle\langle 1| =
             \begin{pmatrix}
                 1 & 0                   & 0 & 0 \\
-                0 & \cos(\th)           & 0 & e^{-i\lambda}\sin(\th) \\
+                0 & \cos(\th)           & 0 & -e^{i\lambda}\sin(\th) \\
                 0 & 0                   & 1 & 0 \\
                 0 & e^{i\phi}\sin(\th)  & 0 & e^{i(\phi+\lambda)\cos(\th)}
             \end{pmatrix}
@@ -169,7 +165,7 @@ class CU3Gate(ControlledGate, metaclass=CU3Meta):
                 \begin{pmatrix}
                     1 & 0   & 0                  & 0 \\
                     0 & 1   & 0                  & 0 \\
-                    0 & 0   & \cos(\th)          & e^{-i\lambda}\sin(\th) \\
+                    0 & 0   & \cos(\th)          & -e^{i\lambda}\sin(\th) \\
                     0 & 0   & e^{i\phi}\sin(\th) & e^{i(\phi+\lambda)\cos(\th)}
                 \end{pmatrix}
     """
@@ -215,18 +211,26 @@ def inverse(self):
         """
         return CU3Gate(-self.params[0], -self.params[2], -self.params[1])
 
-    # TODO: this is the correct definition but has a global phase with respect
-    # to the decomposition above. Restore after allowing phase on circuits.
-    # def to_matrix(self):
-    #    """Return a numpy.array for the CRY gate."""
-    #    theta, phi, lam = self.params
-    #    cos = numpy.cos(theta / 2)
-    #    sin = numpy.sin(theta / 2)
-    #    return numpy.array([[1,0, 0, 0],
-    #                        [0, cos, 0, numpy.exp(-1j * lam) * sin],
-    #                        [0, 0, 1, 0],
-    #                        [0, numpy.exp(1j * phi) * sin, 0, numpy.exp(1j * (phi+lam)) * cos]],
-    #                       dtype=complex)
+    def to_matrix(self):
+        """Return a numpy.array for the CRY gate."""
+        theta, phi, lam = self.params
+        theta, phi, lam = float(theta), float(phi), float(lam)
+        cos = numpy.cos(theta / 2)
+        sin = numpy.sin(theta / 2)
+        if self.ctrl_state:
+            return numpy.array(
+                [[1, 0, 0, 0],
+                 [0, cos, 0, -numpy.exp(1j * lam) * sin],
+                 [0, 0, 1, 0],
+                 [0, numpy.exp(1j * phi) * sin, 0, numpy.exp(1j * (phi+lam)) * cos]],
+                dtype=complex)
+        else:
+            return numpy.array(
+                [[cos, 0, -numpy.exp(1j * lam) * sin, 0],
+                 [0, 1, 0, 0],
+                 [numpy.exp(1j * phi) * sin, 0, numpy.exp(1j * (phi+lam)) * cos, 0],
+                 [0, 0, 0, 1]],
+                dtype=complex)
 
 
 class Cu3Gate(CU3Gate, metaclass=CU3Meta):
diff --git a/qiskit/circuit/library/standard_gates/x.py b/qiskit/circuit/library/standard_gates/x.py
index 545d9df5e678..ebf17cc94cb3 100644
--- a/qiskit/circuit/library/standard_gates/x.py
+++ b/qiskit/circuit/library/standard_gates/x.py
@@ -749,6 +749,12 @@ def inverse(self):
         """Invert this gate. The C4X is its own inverse."""
         return C4XGate()
 
+    def to_matrix(self):
+        """Return a numpy.array for the C4X gate."""
+        return _compute_control_matrix(self.base_gate.to_matrix(),
+                                       self.num_ctrl_qubits,
+                                       ctrl_state=self.ctrl_state)
+
 
 class MCXGate(ControlledGate):
     """The general, multi-controlled X gate."""
diff --git a/qiskit/circuit/quantumcircuit.py b/qiskit/circuit/quantumcircuit.py
index ea4c8ff6d1cc..e63d96afca19 100644
--- a/qiskit/circuit/quantumcircuit.py
+++ b/qiskit/circuit/quantumcircuit.py
@@ -79,6 +79,7 @@ class QuantumCircuit:
 
         name (str): the name of the quantum circuit. If not set, an
             automatically generated string will be assigned.
+        global_phase (float): The global phase of the circuit in radians.
 
     Raises:
         CircuitError: if the circuit name, if given, is not valid.
@@ -136,7 +137,7 @@ class QuantumCircuit:
     header = "OPENQASM 2.0;"
     extension_lib = "include \"qelib1.inc\";"
 
-    def __init__(self, *regs, name=None):
+    def __init__(self, *regs, name=None, global_phase=0):
         if any([not isinstance(reg, (QuantumRegister, ClassicalRegister)) for reg in regs]):
             try:
                 regs = tuple(int(reg) for reg in regs)
@@ -172,6 +173,8 @@ def __init__(self, *regs, name=None):
         self._parameter_table = ParameterTable()
 
         self._layout = None
+        self._global_phase = 0
+        self.global_phase = global_phase
 
     @property
     def data(self):
@@ -362,7 +365,7 @@ def inverse(self):
                  └───────────┘
         """
         inverse_circ = QuantumCircuit(*self.qregs, *self.cregs,
-                                      name=self.name + '_dg')
+                                      name=self.name + '_dg', global_phase=-self.global_phase)
 
         for inst, qargs, cargs in reversed(self._data):
             inverse_circ._append(inst.inverse(), qargs, cargs)
@@ -378,7 +381,8 @@ def repeat(self, reps):
             QuantumCircuit: A circuit containing ``reps`` repetitions of this circuit.
         """
         repeated_circ = QuantumCircuit(*self.qregs, *self.cregs,
-                                       name=self.name + '**{}'.format(reps))
+                                       name=self.name + '**{}'.format(reps),
+                                       global_phase=reps * self.global_phase)
 
         # benefit of appending instructions: decomposing shows the subparts, i.e. the power
         # is actually `reps` times this circuit, and it is currently much faster than `compose`.
@@ -496,6 +500,7 @@ def combine(self, rhs):
         circuit = QuantumCircuit(*combined_qregs, *combined_cregs)
         for instruction_context in itertools.chain(self.data, rhs.data):
             circuit._append(*instruction_context)
+        circuit.global_phase = self.global_phase + rhs.global_phase
         return circuit
 
     def extend(self, rhs):
@@ -535,6 +540,7 @@ def extend(self, rhs):
         # Add new gates
         for instruction_context in data:
             self._append(*instruction_context)
+        self.global_phase += rhs.global_phase
         return self
 
     def compose(self, other, qubits=None, clbits=None, front=False, inplace=False):
@@ -645,6 +651,8 @@ def compose(self, other, qubits=None, clbits=None, front=False, inplace=False):
         for instr, _, _ in mapped_instrs:
             dest._update_parameter_table(instr)
 
+        dest.global_phase += other.global_phase
+
         if inplace:
             return None
 
@@ -1697,6 +1705,30 @@ def from_qasm_str(qasm_str):
         qasm = Qasm(data=qasm_str)
         return _circuit_from_qasm(qasm)
 
+    @property
+    def global_phase(self):
+        """Return the global phase of the circuit in radians."""
+        return self._global_phase
+
+    @global_phase.setter
+    def global_phase(self, angle):
+        """Set the phase of the circuit.
+
+        Args:
+            angle (float, ParameterExpression): radians
+        """
+        if isinstance(angle, ParameterExpression):
+            self._global_phase = angle
+        else:
+            # Set the phase to the [-2 * pi, 2 * pi] interval
+            angle = float(angle)
+            if not angle:
+                self._global_phase = 0
+            elif angle < 0:
+                self._global_phase = angle % (-2 * np.pi)
+            else:
+                self._global_phase = angle % (2 * np.pi)
+
     @property
     def parameters(self):
         """Convenience function to get the parameters defined in the parameter table."""
@@ -1836,6 +1868,10 @@ def _bind_parameter(self, parameter, value):
             # instructions), search the definition for instances of the
             # parameter which also need to be bound.
             self._rebind_definition(instr, parameter, value)
+        # bind circuit's phase
+        if (isinstance(self.global_phase, ParameterExpression) and
+                parameter in self.global_phase.parameters):
+            self.global_phase = self.global_phase.bind({parameter: value})
 
     def _substitute_parameter(self, old_parameter, new_parameter_expr):
         """Substitute an existing parameter in all circuit instructions and the parameter table."""
@@ -1847,6 +1883,9 @@ def _substitute_parameter(self, old_parameter, new_parameter_expr):
         entry = self._parameter_table.pop(old_parameter)
         for new_parameter in new_parameter_expr.parameters:
             self._parameter_table[new_parameter] = entry
+        if (isinstance(self.global_phase, ParameterExpression)
+                and old_parameter in self.global_phase.parameters):
+            self.global_phase = self.global_phase.subs({old_parameter: new_parameter_expr})
 
     def _rebind_definition(self, instruction, parameter, value):
         if instruction._definition:
diff --git a/qiskit/converters/circuit_to_dag.py b/qiskit/converters/circuit_to_dag.py
index 6134d02b1dc8..0d337add13fe 100644
--- a/qiskit/converters/circuit_to_dag.py
+++ b/qiskit/converters/circuit_to_dag.py
@@ -47,6 +47,8 @@ def circuit_to_dag(circuit):
     """
     dagcircuit = DAGCircuit()
     dagcircuit.name = circuit.name
+    dagcircuit.global_phase = circuit.global_phase
+
     for register in circuit.qregs:
         dagcircuit.add_qreg(register)
     for register in circuit.cregs:
diff --git a/qiskit/converters/circuit_to_gate.py b/qiskit/converters/circuit_to_gate.py
index bf4075edcb83..2fd7be869c61 100644
--- a/qiskit/converters/circuit_to_gate.py
+++ b/qiskit/converters/circuit_to_gate.py
@@ -14,7 +14,6 @@
 
 
 """Helper function for converting a circuit to a gate"""
-
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.quantumregister import QuantumRegister, Qubit
 from qiskit.exceptions import QiskitError
@@ -103,7 +102,7 @@ def find_bit_position(bit):
                    list(map(lambda y: q[find_bit_position(y)], x[1])),
                    []),
         rules))
-    qc = QuantumCircuit(q, name=gate.name)
+    qc = QuantumCircuit(q, name=gate.name, global_phase=target.global_phase)
     qc._data = rules
     gate.definition = qc
     return gate
diff --git a/qiskit/converters/circuit_to_instruction.py b/qiskit/converters/circuit_to_instruction.py
index f878fef1797e..10fed116b218 100644
--- a/qiskit/converters/circuit_to_instruction.py
+++ b/qiskit/converters/circuit_to_instruction.py
@@ -124,6 +124,9 @@ def find_bit_position(bit):
 
     qc = QuantumCircuit(*regs, name=instruction.name)
     qc._data = definition
+    if circuit.global_phase:
+        qc.global_phase = circuit.global_phase
+
     instruction.definition = qc
 
     return instruction
diff --git a/qiskit/converters/dag_to_circuit.py b/qiskit/converters/dag_to_circuit.py
index 32e429559cf5..ed9046fdb3f8 100644
--- a/qiskit/converters/dag_to_circuit.py
+++ b/qiskit/converters/dag_to_circuit.py
@@ -48,7 +48,8 @@ def dag_to_circuit(dag):
     """
 
     name = dag.name or None
-    circuit = QuantumCircuit(*dag.qregs.values(), *dag.cregs.values(), name=name)
+    circuit = QuantumCircuit(*dag.qregs.values(), *dag.cregs.values(), name=name,
+                             global_phase=dag.global_phase)
 
     for node in dag.topological_op_nodes():
         # Get arguments for classical control (if any)
diff --git a/qiskit/dagcircuit/dagcircuit.py b/qiskit/dagcircuit/dagcircuit.py
index d7ecbe25a01d..a17af95f3476 100644
--- a/qiskit/dagcircuit/dagcircuit.py
+++ b/qiskit/dagcircuit/dagcircuit.py
@@ -26,6 +26,7 @@
 import copy
 import itertools
 import warnings
+import math
 
 import retworkx as rx
 import networkx as nx
@@ -90,6 +91,8 @@ def __call__(self):
         self._qubits = DummyCallableList()  # TODO: make these a regular empty list [] after the
         self._clbits = DummyCallableList()  # DeprecationWarning period, and remove name underscore.
 
+        self._global_phase = 0
+
     def to_networkx(self):
         """Returns a copy of the DAGCircuit in networkx format."""
         G = nx.MultiDiGraph()
@@ -153,6 +156,31 @@ def node_counter(self):
         """
         return len(self._multi_graph)
 
+    @property
+    def global_phase(self):
+        """Return the global phase of the circuit."""
+        return self._global_phase
+
+    @global_phase.setter
+    def global_phase(self, angle):
+        """Set the global phase of the circuit.
+
+        Args:
+            angle (float, ParameterExpression)
+        """
+        from qiskit.circuit.parameterexpression import ParameterExpression  # needed?
+        if isinstance(angle, ParameterExpression):
+            self._global_phase = angle
+        else:
+            # Set the phase to the [-2 * pi, 2 * pi] interval
+            angle = float(angle)
+            if not angle:
+                self._global_phase = 0
+            elif angle < 0:
+                self._global_phase = angle % (-2 * math.pi)
+            else:
+                self._global_phase = angle % (2 * math.pi)
+
     def remove_all_ops_named(self, opname):
         """Remove all operation nodes with the given name."""
         for n in self.named_nodes(opname):
@@ -571,6 +599,7 @@ def compose(self, other, edge_map=None, qubits=None, clbits=None, front=False, i
             dag = self
         else:
             dag = copy.deepcopy(self)
+        dag.global_phase += other.global_phase
 
         for nd in other.topological_nodes():
             if nd.type == "in":
diff --git a/qiskit/extensions/quantum_initializer/initializer.py b/qiskit/extensions/quantum_initializer/initializer.py
index 80058ba041ae..baaad25688b1 100644
--- a/qiskit/extensions/quantum_initializer/initializer.py
+++ b/qiskit/extensions/quantum_initializer/initializer.py
@@ -25,7 +25,7 @@
 from qiskit.circuit import Instruction
 from qiskit.circuit.library.standard_gates.x import CXGate
 from qiskit.circuit.library.standard_gates.ry import RYGate
-from qiskit.circuit.library.standard_gates.rz import RZGate
+from qiskit.circuit.library.standard_gates.u1 import U1Gate
 from qiskit.circuit.reset import Reset
 
 _EPS = 1e-10  # global variable used to chop very small numbers to zero
@@ -113,7 +113,7 @@ def gates_to_uncompute(self):
                 add_last_cnot = False
 
             if np.linalg.norm(phis) != 0:
-                rz_mult = self._multiplex(RZGate, phis, last_cnot=add_last_cnot)
+                rz_mult = self._multiplex(U1Gate, phis, last_cnot=add_last_cnot)
                 circuit.append(rz_mult.to_instruction(), q[i:self.num_qubits])
 
             if np.linalg.norm(thetas) != 0:
diff --git a/qiskit/extensions/unitary.py b/qiskit/extensions/unitary.py
index 14add1df1e12..c4b606ba4619 100644
--- a/qiskit/extensions/unitary.py
+++ b/qiskit/extensions/unitary.py
@@ -132,13 +132,26 @@ def control(self, num_ctrl_qubits=1, label=None, ctrl_state=None):
             UnitaryGate: controlled version of gate.
 
         Raises:
-            QiskitError: invalid ctrl_state
+            QiskitError: Invalid ctrl_state.
+            ExtensionError: Non-unitary controlled unitary.
         """
-        cmat = _compute_control_matrix(self.to_matrix(), num_ctrl_qubits)
+        cmat = _compute_control_matrix(self.to_matrix(), num_ctrl_qubits, ctrl_state=ctrl_state)
         iso = isometry.Isometry(cmat, 0, 0)
         cunitary = ControlledGate('c-unitary', num_qubits=self.num_qubits+num_ctrl_qubits,
                                   params=[cmat], label=label, num_ctrl_qubits=num_ctrl_qubits,
                                   definition=iso.definition, ctrl_state=ctrl_state)
+
+        from qiskit.quantum_info import Operator
+        # hack to correct global phase; should fix to prevent need for correction here
+        pmat = (Operator(iso.inverse()).data @ cmat)
+        diag = numpy.diag(pmat)
+        if not numpy.allclose(diag, diag[0]):
+            raise ExtensionError('controlled unitary generation failed')
+        phase = numpy.angle(diag[0])
+        if phase:
+            qreg = cunitary.definition.qregs[0]
+            cunitary.definition.u3(numpy.pi, phase, phase - numpy.pi, qreg[0])
+            cunitary.definition.u3(numpy.pi, 0, numpy.pi, qreg[0])
         cunitary.base_gate = self.copy()
         cunitary.base_gate.label = self.label
         return cunitary
diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py
index 281229395298..495b9561200e 100644
--- a/qiskit/quantum_info/operators/operator.py
+++ b/qiskit/quantum_info/operators/operator.py
@@ -485,11 +485,12 @@ def _einsum_matmul(cls, tensor, mat, indices, shift=0, right_mul=False):
     @classmethod
     def _init_instruction(cls, instruction):
         """Convert a QuantumCircuit or Instruction to an Operator."""
+        # Initialize an identity operator of the correct size of the circuit
+        dimension = 2 ** instruction.num_qubits
+        op = Operator(np.eye(dimension))
         # Convert circuit to an instruction
         if isinstance(instruction, QuantumCircuit):
             instruction = instruction.to_instruction()
-        # Initialize an identity operator of the correct size of the circuit
-        op = Operator(np.eye(2 ** instruction.num_qubits))
         op._append_instruction(instruction)
         return op
 
@@ -511,6 +512,7 @@ def _instruction_to_matrix(cls, obj):
     def _append_instruction(self, obj, qargs=None):
         """Update the current Operator by apply an instruction."""
         from qiskit.circuit.barrier import Barrier
+        from .scalar_op import ScalarOp
 
         mat = self._instruction_to_matrix(obj)
         if mat is not None:
@@ -530,6 +532,12 @@ def _append_instruction(self, obj, qargs=None):
                 raise QiskitError('Instruction "{}" '
                                   'definition is {} but expected QuantumCircuit.'.format(
                                       obj.name, type(obj.definition)))
+            if obj.definition.global_phase:
+                dimension = 2 ** self.num_qubits
+                op = self.compose(
+                    ScalarOp(dimension, np.exp(1j * float(obj.definition.global_phase))),
+                    qargs=qargs)
+                self._data = op.data
             flat_instr = obj.definition.to_instruction()
             for instr, qregs, cregs in flat_instr.definition.data:
                 if cregs:
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index 4526b68da124..4e1b6bacc0fb 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -677,6 +677,8 @@ def _evolve_instruction(statevec, obj, qargs=None):
         if not isinstance(obj.definition, QuantumCircuit):
             raise QiskitError('{0} instruction definition is {1}; expected QuantumCircuit'.format(
                 obj.name, type(obj.definition)))
+        if obj.definition.global_phase:
+            statevec._data *= np.exp(1j * obj.definition.global_phase)
         for instr, qregs, cregs in obj.definition:
             if cregs:
                 raise QiskitError(
diff --git a/qiskit/quantum_info/synthesis/ion_decompose.py b/qiskit/quantum_info/synthesis/ion_decompose.py
index 94ee549a9f85..c42ecb15b106 100644
--- a/qiskit/quantum_info/synthesis/ion_decompose.py
+++ b/qiskit/quantum_info/synthesis/ion_decompose.py
@@ -49,8 +49,7 @@ def cnot_rxx_decompose(plus_ry=True, plus_rxx=True):
         sgn_rxx = 1
     else:
         sgn_rxx = -1
-
-    circuit = QuantumCircuit(2)
+    circuit = QuantumCircuit(2, global_phase=-sgn_ry * sgn_rxx * np.pi / 4)
     circuit.append(RYGate(sgn_ry * np.pi / 2), [0])
     circuit.append(RXXGate(sgn_rxx * np.pi / 2), [0, 1])
     circuit.append(RXGate(-sgn_rxx * np.pi / 2), [0])
diff --git a/qiskit/transpiler/passes/basis/unroller.py b/qiskit/transpiler/passes/basis/unroller.py
index 7a8158e96d50..51bae5b6b9e9 100644
--- a/qiskit/transpiler/passes/basis/unroller.py
+++ b/qiskit/transpiler/passes/basis/unroller.py
@@ -77,6 +77,8 @@ def run(self, dag):
             # different that the width of the node.
             while rule and len(rule) == 1 and len(node.qargs) == len(rule[0][1]):
                 if rule[0][0].name in self.basis:
+                    if node.op.definition and node.op.definition.global_phase:
+                        dag.global_phase += node.op.definition.global_phase
                     dag.substitute_node(node, rule[0][0], inplace=True)
                     break
                 try:
@@ -95,5 +97,10 @@ def run(self, dag):
                                       (str(self.basis), node.op.name))
                 decomposition = circuit_to_dag(node.op.definition)
                 unrolled_dag = self.run(decomposition)  # recursively unroll ops
+                if node.op.definition and node.op.definition.global_phase:
+                    dag.global_phase += node.op.definition.global_phase
+                if unrolled_dag.global_phase:
+                    dag.global_phase += unrolled_dag.global_phase
+                    unrolled_dag.global_phase = 0
                 dag.substitute_node_with_dag(node, unrolled_dag)
         return dag
diff --git a/releasenotes/notes/add_quantumcircuit_phase-5006d1e930348d2e.yaml b/releasenotes/notes/add_quantumcircuit_phase-5006d1e930348d2e.yaml
new file mode 100644
index 000000000000..255913a55734
--- /dev/null
+++ b/releasenotes/notes/add_quantumcircuit_phase-5006d1e930348d2e.yaml
@@ -0,0 +1,17 @@
+features:
+  - |
+    Adds a `global_phase` attribute to `QuantumCircuit` for tracking global
+    phase. This allows, for instance, to have the `to_matrix` method
+    of a gate to exactly correspond to its decompositions instead of
+    just up to a global phase. The same attribute has also been
+    extended to DAGCircuit so the phase can be tracked when converting
+    between QuantumCircuit and DAGCircuit. For example::
+
+      from qiskit import QuantumCircuit
+      circ = QuantumCircuit(1, phase=math.pi)
+      circ.u1(0)
+
+    The global phase may also be set or queried with the `circ.global_phase`
+    property. In either case the setting is in radians. If the circuit
+    is converted to an instruction or gate the global phase is
+    represented by two single qubit rotations on the first qubit.
diff --git a/test/python/circuit/test_circuit_operations.py b/test/python/circuit/test_circuit_operations.py
index 440159dc52d1..28dec4a066e0 100644
--- a/test/python/circuit/test_circuit_operations.py
+++ b/test/python/circuit/test_circuit_operations.py
@@ -574,6 +574,21 @@ def test_cnot_alias(self):
         expected.cx(0, 1)
         self.assertEqual(qc, expected)
 
+    def test_inverse(self):
+        """Test inverse circuit."""
+        qr = QuantumRegister(2)
+        qc = QuantumCircuit(qr, global_phase=0.5)
+        qc.h(0)
+        qc.barrier(qr)
+        qc.t(1)
+
+        expected = QuantumCircuit(qr)
+        expected.tdg(1)
+        expected.barrier(qr)
+        expected.h(0)
+        expected.global_phase = -0.5
+        self.assertEqual(qc.inverse(), expected)
+
 
 class TestCircuitBuilding(QiskitTestCase):
     """QuantumCircuit tests."""
diff --git a/test/python/circuit/test_compose.py b/test/python/circuit/test_compose.py
index c2bc3bd1a83e..16d9a7644493 100644
--- a/test/python/circuit/test_compose.py
+++ b/test/python/circuit/test_compose.py
@@ -451,6 +451,16 @@ def test_compose_one_liner(self):
 
         self.assertEqual(circ, expected)
 
+    def test_compose_global_phase(self):
+        """Composing with global phase."""
+        circ1 = QuantumCircuit(1, global_phase=1)
+        circ1.rz(0.5, 0)
+        circ2 = QuantumCircuit(1, global_phase=2)
+        circ3 = QuantumCircuit(1, global_phase=3)
+        circ4 = circ1.compose(circ2).compose(circ3)
+        self.assertEqual(circ4.global_phase,
+                         circ1.global_phase + circ2.global_phase + circ3.global_phase)
+
     def test_compose_front_circuit(self):
         """Test composing a circuit at the front of a circuit.
         """
diff --git a/test/python/circuit/test_controlled_gate.py b/test/python/circuit/test_controlled_gate.py
index aa88e8adbda1..f17ce5bc0e35 100644
--- a/test/python/circuit/test_controlled_gate.py
+++ b/test/python/circuit/test_controlled_gate.py
@@ -578,7 +578,7 @@ def test_inverse_x(self, num_ctrl_qubits):
         np.testing.assert_array_almost_equal(result.data, np.identity(result.dim[0]))
 
     @data(1, 2, 3)
-    def test_inverse_circuit(self, num_ctrl_qubits):
+    def test_inverse_gate(self, num_ctrl_qubits):
         """Test inverting a controlled gate based on a circuit definition."""
         qc = QuantumCircuit(3)
         qc.h(0)
@@ -591,6 +591,19 @@ def test_inverse_circuit(self, num_ctrl_qubits):
         result = Operator(cgate).compose(Operator(inv_cgate))
         np.testing.assert_array_almost_equal(result.data, np.identity(result.dim[0]))
 
+    @data(1, 2, 3)
+    def test_inverse_circuit(self, num_ctrl_qubits):
+        """Test inverting a controlled gate based on a circuit definition."""
+        qc = QuantumCircuit(3)
+        qc.h(0)
+        qc.cx(0, 1)
+        qc.cx(1, 2)
+        qc.rx(np.pi / 4, [0, 1, 2])
+        cqc = qc.control(num_ctrl_qubits)
+        cqc_inv = cqc.inverse()
+        result = Operator(cqc).compose(Operator(cqc_inv))
+        np.testing.assert_array_almost_equal(result.data, np.identity(result.dim[0]))
+
     @data(1, 2, 3, 4, 5)
     def test_controlled_unitary(self, num_ctrl_qubits):
         """Test the matrix data of an Operator, which is based on a controlled gate."""
@@ -621,13 +634,14 @@ def test_controlled_random_unitary(self, num_ctrl_qubits):
         cop_mat = _compute_control_matrix(base_mat, num_ctrl_qubits)
         self.assertTrue(matrix_equal(cop_mat, test_op.data, ignore_phase=True))
 
-    @combine(num_ctrl_qubits=[1, 2, 3], ctrl_state=[None, 0])
+    @combine(num_ctrl_qubits=[1, 2, 3], ctrl_state=[None])
     def test_open_controlled_unitary_z(self, num_ctrl_qubits, ctrl_state):
         """Test that UnitaryGate with control returns params."""
         umat = np.array([[1, 0], [0, -1]])
-        ugate = UnitaryGate(umat).control(num_ctrl_qubits, ctrl_state=ctrl_state)
+        ugate = UnitaryGate(umat)
+        cugate = ugate.control(num_ctrl_qubits, ctrl_state=ctrl_state)
         ref_mat = _compute_control_matrix(umat, num_ctrl_qubits, ctrl_state=ctrl_state)
-        self.assertTrue(matrix_equal(Operator(ugate).data, ref_mat))
+        self.assertEqual(Operator(cugate), Operator(ref_mat))
 
     @data(1, 2, 3)
     def test_open_controlled_unitary_matrix(self, num_ctrl_qubits):
@@ -765,7 +779,6 @@ def test_all_inverses(self, gate):
             try:
                 numargs = len(_get_free_params(gate))
                 args = [2] * numargs
-
                 gate = gate(*args)
                 self.assertEqual(gate.inverse().control(2), gate.control(2).inverse())
             except AttributeError:
@@ -887,6 +900,58 @@ def test_open_controlled_equality(self):
             XGate().control(1, ctrl_state='0'),
             XGate().control(1, ctrl_state='1'))
 
+    def test_cx_global_phase(self):
+        """
+        Test controlling CX with global phase
+        """
+        theta = pi/2
+        circ = QuantumCircuit(2, global_phase=theta)
+        circ.cx(0, 1)
+        cx = circ.to_gate()
+        self.assertNotEqual(Operator(CXGate()), Operator(cx))
+
+        ccx = cx.control(1)
+        base_mat = Operator(cx).data
+        target = _compute_control_matrix(base_mat, 1)
+        self.assertEqual(Operator(ccx), Operator(target))
+
+        expected = QuantumCircuit(*ccx.definition.qregs)
+        expected.ccx(0, 1, 2)
+        expected.u1(theta, 0)
+        self.assertEqual(ccx.definition, expected)
+
+    @data(1, 2)
+    def test_controlled_global_phase(self, num_ctrl_qubits):
+        """
+        Test controlled global phase on base gate.
+        """
+        theta = pi/4
+        circ = QuantumCircuit(2, global_phase=theta)
+        base_gate = circ.to_gate()
+        base_mat = Operator(base_gate).data
+        target = _compute_control_matrix(base_mat, num_ctrl_qubits)
+        cgate = base_gate.control(num_ctrl_qubits)
+        ccirc = circ.control(num_ctrl_qubits)
+        self.assertEqual(Operator(cgate), Operator(target))
+        self.assertEqual(Operator(ccirc), Operator(target))
+
+    @data(1, 2)
+    def test_rz_composite_global_phase(self, num_ctrl_qubits):
+        """
+        Test controlling CX with global phase
+        """
+        theta = pi/4
+        circ = QuantumCircuit(2, global_phase=theta)
+        circ.rz(0.1, 0)
+        circ.rz(0.2, 1)
+        ccirc = circ.control(num_ctrl_qubits)
+        base_gate = circ.to_gate()
+        cgate = base_gate.control(num_ctrl_qubits)
+        base_mat = Operator(base_gate).data
+        target = _compute_control_matrix(base_mat, num_ctrl_qubits)
+        self.assertEqual(Operator(cgate), Operator(target))
+        self.assertEqual(Operator(ccirc), Operator(target))
+
 
 @ddt
 class TestOpenControlledToMatrix(QiskitTestCase):
@@ -1009,6 +1074,7 @@ def test_controlled_standard_gates(self, num_ctrl_qubits, gate_class):
             args = [5]
 
         gate = gate_class(*args)
+
         for ctrl_state in {ctrl_state_ones, ctrl_state_zeros, ctrl_state_mixed}:
             with self.subTest(i='{0}, ctrl_state={1}'.format(gate_class.__name__,
                                                              ctrl_state)):
@@ -1029,7 +1095,7 @@ def test_controlled_standard_gates(self, num_ctrl_qubits, gate_class):
                     base_mat = Operator(gate).data
                 target_mat = _compute_control_matrix(base_mat, num_ctrl_qubits,
                                                      ctrl_state=ctrl_state)
-                self.assertTrue(matrix_equal(Operator(cgate).data, target_mat, ignore_phase=True))
+                self.assertEqual(Operator(cgate), Operator(target_mat))
 
 
 @ddt
diff --git a/test/python/circuit/test_extensions_standard.py b/test/python/circuit/test_extensions_standard.py
index 0ec21175a9a0..50d09b6df71e 100644
--- a/test/python/circuit/test_extensions_standard.py
+++ b/test/python/circuit/test_extensions_standard.py
@@ -1388,6 +1388,45 @@ def test_to_matrix(self):
             self.assertTrue(matrix_equal(definition_unitary, gate_matrix))
             self.assertTrue(is_unitary_matrix(gate_matrix))
 
+    def test_to_matrix_op(self):
+        """test gates implementing to_matrix generate matrix which matches
+        definition using Operator."""
+        from qiskit.quantum_info import Operator
+        from qiskit.circuit.library.standard_gates.ms import MSGate
+
+        params = [0.1 * i for i in range(10)]
+        gate_class_list = Gate.__subclasses__() + ControlledGate.__subclasses__()
+        for gate_class in gate_class_list:
+            sig = signature(gate_class)
+            if gate_class == MSGate:
+                # due to the signature (num_qubits, theta, *, n_qubits=Noe) the signature detects
+                # 3 arguments but really its only 2. This if can be removed once the deprecated
+                # n_qubits argument is no longer supported.
+                free_params = 2
+            else:
+                free_params = len(set(sig.parameters) - {'label'})
+            try:
+                gate = gate_class(*params[0:free_params])
+            except (CircuitError, QiskitError, AttributeError):
+                self.log.info(
+                    'Cannot init gate with params only. Skipping %s',
+                    gate_class)
+                continue
+            if gate.name in ['U', 'CX']:
+                continue
+            try:
+                gate_matrix = gate.to_matrix()
+            except CircuitError:
+                # gate doesn't implement to_matrix method: skip
+                self.log.info('to_matrix method FAILED for "%s" gate',
+                              gate.name)
+                continue
+            if not hasattr(gate, 'definition') or not gate.definition:
+                continue
+            definition_unitary = Operator(gate.definition).data
+            self.assertTrue(matrix_equal(definition_unitary, gate_matrix))
+            self.assertTrue(is_unitary_matrix(gate_matrix))
+
 
 @ddt
 class TestQubitKeywordArgRenaming(QiskitTestCase):
diff --git a/test/python/circuit/test_gate_definitions.py b/test/python/circuit/test_gate_definitions.py
index 7fd98d9ff5b3..8d568a5da5a7 100644
--- a/test/python/circuit/test_gate_definitions.py
+++ b/test/python/circuit/test_gate_definitions.py
@@ -15,12 +15,13 @@
 
 """Test hardcoded decomposition rules and matrix definitions for standard gates."""
 
+import numpy as np
 from ddt import ddt, data
 
 from qiskit import QuantumCircuit
 from qiskit.quantum_info import Operator
 from qiskit.test import QiskitTestCase
-from qiskit.circuit import ParameterVector
+from qiskit.circuit import ParameterVector, Gate, ControlledGate
 
 
 from qiskit.circuit.library import (
@@ -39,7 +40,8 @@
 
 class TestGateDefinitions(QiskitTestCase):
     """Test the decomposition of a gate in terms of other gates
-    yields the same matrix as the hardcoded matrix definition."""
+    yields the equivalent matrix as the hardcoded matrix definition
+    up to a global phase."""
 
     def test_ch_definition(self):  # TODO: expand this to all gates
         """Test ch gate matrix and definition.
@@ -114,6 +116,45 @@ def test_cx_definition(self):
         self.assertTrue(Operator(circ).equiv(Operator(decomposed_circ)))
 
 
+class TestGateEquivalenceEqual(QiskitTestCase):
+    """Test the decomposition of a gate in terms of other gates
+    yields the same matrix as the hardcoded matrix definition."""
+
+    @classmethod
+    def setUpClass(cls):
+        class_list = Gate.__subclasses__() + ControlledGate.__subclasses__()
+        exclude = {'ControlledGate', 'DiagonalGate', 'UCGate', 'MCGupDiag',
+                   'MCU1Gate', 'UnitaryGate', 'HamiltonianGate',
+                   'UCPauliRotGate', 'SingleQubitUnitary', 'MCXGate',
+                   'VariadicZeroParamGate'}
+        cls._gate_classes = []
+        for aclass in class_list:
+            if aclass.__name__ not in exclude:
+                cls._gate_classes.append(aclass)
+
+    def test_equivalence_phase(self):
+        """Test that the equivalent circuits from the equivalency_library
+        have equal matrix representations"""
+        for gate_class in self._gate_classes:
+            with self.subTest(i=gate_class):
+                n_params = len(_get_free_params(gate_class))
+                params = [0.1 * i for i in range(1, n_params+1)]
+                if gate_class.__name__ == 'RXXGate':
+                    params = [np.pi/2]
+                params = [-np.pi/2 * i for i in range(1, n_params+1)]
+                if gate_class.__name__ in ['MSGate']:
+                    params[0] = 2
+                elif gate_class in ['MCU1Gate']:
+                    params[1] = 2
+                gate = gate_class(*params)
+                equiv_lib_list = std_eqlib.get_entry(gate)
+                for ieq, equivalency in enumerate(equiv_lib_list):
+                    with self.subTest(msg=gate.name + '_' + str(ieq)):
+                        op1 = Operator(gate)
+                        op2 = Operator(equivalency)
+                        self.assertEqual(op1, op2)
+
+
 @ddt
 class TestStandardEquivalenceLibrary(QiskitTestCase):
     """Standard Extension Test."""
diff --git a/test/python/circuit/test_gate_power.py b/test/python/circuit/test_gate_power.py
index 483da2177af1..bd3b74a94d7b 100644
--- a/test/python/circuit/test_gate_power.py
+++ b/test/python/circuit/test_gate_power.py
@@ -114,19 +114,31 @@ def test_composite_sqrt(self):
         """Test composite Gate.power(1/2) method.
         """
         circ = QuantumCircuit(1, name='my_gate')
-        circ.rz(0.1, 0)
-        circ.rx(0.2, 0)
+        import numpy as np
+        thetaz = 0.1
+        thetax = 0.2
+        circ.rz(thetaz, 0)
+        circ.rx(thetax, 0)
         gate = circ.to_gate()
 
-        expected = array([[0.99874948+6.25390559e-05j, 0.00374609-4.98542083e-02j],
-                          [-0.00124974-4.99791301e-02j, 0.99750443+4.98542083e-02j]])
-
         result = gate.power(1 / 2)
 
+        iden = Operator.from_label('I')
+        xgen = Operator.from_label('X')
+        zgen = Operator.from_label('Z')
+
+        def rzgate(theta):
+            return np.cos(0.5 * theta) * iden - 1j * np.sin(0.5 * theta) * zgen
+
+        def rxgate(theta):
+            return np.cos(0.5 * theta) * iden - 1j * np.sin(0.5 * theta) * xgen
+
+        rxrz = rxgate(thetax) * rzgate(thetaz)
+
         self.assertEqual(result.label, 'my_gate^0.5')
         self.assertEqual(len(result.definition), 1)
         self.assertIsInstance(result, Gate)
-        self.assertEqual(Operator(result), Operator(expected))
+        self.assertEqual(Operator(result) @ Operator(result), rxrz)
 
 
 @ddt
diff --git a/test/python/compiler/test_transpiler.py b/test/python/compiler/test_transpiler.py
index c7cc6ef851a4..f093041c0481 100644
--- a/test/python/compiler/test_transpiler.py
+++ b/test/python/compiler/test_transpiler.py
@@ -454,7 +454,6 @@ def test_parameter_expression_circuit_for_simulator(self):
 
         expected_qc = QuantumCircuit(qr)
         expected_qc.u1(square, qr[0])
-
         self.assertEqual(expected_qc, transpiled_qc)
 
     def test_parameter_expression_circuit_for_device(self):
@@ -473,7 +472,6 @@ def test_parameter_expression_circuit_for_device(self):
         qr = QuantumRegister(14, 'q')
         expected_qc = QuantumCircuit(qr)
         expected_qc.u1(square, qr[0])
-
         self.assertEqual(expected_qc, transpiled_qc)
 
     def test_final_measurement_barrier_for_devices(self):
diff --git a/test/python/quantum_info/states/test_statevector.py b/test/python/quantum_info/states/test_statevector.py
index eeed6f49365d..fa3b36f9f607 100644
--- a/test/python/quantum_info/states/test_statevector.py
+++ b/test/python/quantum_info/states/test_statevector.py
@@ -266,6 +266,17 @@ def test_evolve_subsystem(self):
             target = Statevector(np.dot(op_full.data, vec))
             self.assertEqual(state.evolve(op, qargs=[2, 1, 0]), target)
 
+    def test_evolve_global_phase(self):
+        """Test evolve circuit with global phase."""
+        state_i = Statevector([1, 0])
+        qr = QuantumRegister(2)
+        phase = np.pi / 4
+        circ = QuantumCircuit(qr, global_phase=phase)
+        circ.x(0)
+        state_f = state_i.evolve(circ, qargs=[0])
+        target = Statevector([0, 1]) * np.exp(1j * phase)
+        self.assertEqual(state_f, target)
+
     def test_conjugate(self):
         """Test conjugate method."""
         for _ in range(10):
diff --git a/test/python/transpiler/test_unroller.py b/test/python/transpiler/test_unroller.py
index 35ae020afd5c..d7d0e8439e89 100644
--- a/test/python/transpiler/test_unroller.py
+++ b/test/python/transpiler/test_unroller.py
@@ -92,6 +92,7 @@ def test_unroll_1q_chain_conditional(self):
         ref_circuit.u3(pi, pi/2, pi/2, qr[0]).c_if(cr, 1)
         ref_circuit.u1(pi, qr[0]).c_if(cr, 1)
         ref_dag = circuit_to_dag(ref_circuit)
+
         self.assertEqual(unrolled_dag, ref_dag)
 
     def test_unroll_no_basis(self):
@@ -107,143 +108,6 @@ def test_unroll_no_basis(self):
         with self.assertRaises(QiskitError):
             pass_.run(dag)
 
-    def test_unroll_all_instructions(self):
-        """Test unrolling a circuit containing all standard instructions.
-        """
-        qr = QuantumRegister(3, 'qr')
-        cr = ClassicalRegister(3, 'cr')
-        circuit = QuantumCircuit(qr, cr)
-        circuit.crx(0.5, qr[1], qr[2])
-        circuit.cry(0.5, qr[1], qr[2])
-        circuit.ccx(qr[0], qr[1], qr[2])
-        circuit.ch(qr[0], qr[2])
-        circuit.crz(0.5, qr[1], qr[2])
-        circuit.cswap(qr[1], qr[0], qr[2])
-        circuit.cu1(0.1, qr[0], qr[2])
-        circuit.cu3(0.2, 0.1, 0.0, qr[1], qr[2])
-        circuit.cx(qr[1], qr[0])
-        circuit.cy(qr[1], qr[2])
-        circuit.cz(qr[2], qr[0])
-        circuit.h(qr[1])
-        circuit.i(qr[0])
-        circuit.rx(0.1, qr[0])
-        circuit.ry(0.2, qr[1])
-        circuit.rz(0.3, qr[2])
-        circuit.rzz(0.6, qr[1], qr[0])
-        circuit.s(qr[0])
-        circuit.sdg(qr[1])
-        circuit.swap(qr[1], qr[2])
-        circuit.t(qr[2])
-        circuit.tdg(qr[0])
-        circuit.u1(0.1, qr[1])
-        circuit.u2(0.2, -0.1, qr[0])
-        circuit.u3(0.3, 0.0, -0.1, qr[2])
-        circuit.x(qr[2])
-        circuit.y(qr[1])
-        circuit.z(qr[0])
-        circuit.snapshot('0')
-        circuit.measure(qr, cr)
-        dag = circuit_to_dag(circuit)
-        pass_ = Unroller(basis=['u3', 'cx', 'id'])
-        unrolled_dag = pass_.run(dag)
-
-        ref_circuit = QuantumCircuit(qr, cr)
-        ref_circuit.u3(0, 0, pi/2, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(-0.25, 0, 0, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0.25, -pi/2, 0, qr[2])
-        ref_circuit.u3(0.25, 0, 0, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(-0.25, 0, 0, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0, 0, -pi/4, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0, 0, pi/4, qr[1])
-        ref_circuit.u3(0, 0, -pi/4, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.cx(qr[0], qr[1])
-        ref_circuit.u3(0, 0, pi/4, qr[0])
-        ref_circuit.u3(0, 0, -pi/4, qr[1])
-        ref_circuit.cx(qr[0], qr[1])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.u3(0, 0, pi/2, qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, -pi/4, qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.u3(0, 0, -pi/2, qr[2])
-        ref_circuit.u3(0, 0, 0.25, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0, 0, -0.25, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.cx(qr[2], qr[0])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, -pi/4, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, pi/4, qr[0])
-        ref_circuit.u3(0, 0, -pi/4, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.cx(qr[1], qr[0])
-        ref_circuit.u3(0, 0, -pi/4, qr[0])
-        ref_circuit.u3(0, 0, pi/4, qr[1])
-        ref_circuit.cx(qr[1], qr[0])
-        ref_circuit.u3(0, 0, 0.05, qr[1])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[2])
-        ref_circuit.cx(qr[2], qr[0])
-        ref_circuit.u3(0, 0, 0.05, qr[0])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, -0.05, qr[2])
-        ref_circuit.cx(qr[0], qr[2])
-        ref_circuit.u3(0, 0, 0.05, qr[2])
-        ref_circuit.u3(0, 0, -0.05, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(-0.1, 0, -0.05, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.cx(qr[1], qr[0])
-        ref_circuit.u3(pi/2, 0, pi, qr[0])
-        ref_circuit.u3(0.1, 0.1, 0, qr[2])
-        ref_circuit.u3(0, 0, -pi/2, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(pi/2, 0, pi, qr[1])
-        ref_circuit.u3(0.2, 0, 0, qr[1])
-        ref_circuit.u3(0, 0, pi/2, qr[2])
-        ref_circuit.cx(qr[2], qr[0])
-        ref_circuit.u3(pi/2, 0, pi, qr[0])
-        ref_circuit.i(qr[0])
-        ref_circuit.u3(0.1, -pi/2, pi/2, qr[0])
-        ref_circuit.cx(qr[1], qr[0])
-        ref_circuit.u3(0, 0, 0.6, qr[0])
-        ref_circuit.cx(qr[1], qr[0])
-        ref_circuit.u3(0, 0, pi/2, qr[0])
-        ref_circuit.u3(0, 0, -pi/4, qr[0])
-        ref_circuit.u3(pi/2, 0.2, -0.1, qr[0])
-        ref_circuit.u3(0, 0, pi, qr[0])
-        ref_circuit.u3(0, 0, -pi/2, qr[1])
-        ref_circuit.u3(0, 0, 0.3, qr[2])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.cx(qr[2], qr[1])
-        ref_circuit.cx(qr[1], qr[2])
-        ref_circuit.u3(0, 0, 0.1, qr[1])
-        ref_circuit.u3(pi, pi/2, pi/2, qr[1])
-        ref_circuit.u3(0, 0, pi/4, qr[2])
-        ref_circuit.u3(0.3, 0.0, -0.1, qr[2])
-        ref_circuit.u3(pi, 0, pi, qr[2])
-        ref_circuit.snapshot('0')
-        ref_circuit.measure(qr, cr)
-        ref_dag = circuit_to_dag(ref_circuit)
-        self.assertEqual(unrolled_dag, ref_dag)
-
     def test_simple_unroll_parameterized_without_expressions(self):
         """Verify unrolling parameterized gates without expressions."""
         qr = QuantumRegister(1)
@@ -254,7 +118,7 @@ def test_simple_unroll_parameterized_without_expressions(self):
         qc.rz(theta, qr[0])
         dag = circuit_to_dag(qc)
 
-        unrolled_dag = Unroller(['u1', 'cx']).run(dag)
+        unrolled_dag = Unroller(['u1', 'u3', 'cx']).run(dag)
 
         expected = QuantumCircuit(qr)
         expected.u1(theta, qr[0])
@@ -273,7 +137,7 @@ def test_simple_unroll_parameterized_with_expressions(self):
         qc.rz(sum_, qr[0])
         dag = circuit_to_dag(qc)
 
-        unrolled_dag = Unroller(['u1', 'cx']).run(dag)
+        unrolled_dag = Unroller(['u1', 'u3', 'cx']).run(dag)
 
         expected = QuantumCircuit(qr)
         expected.u1(sum_, qr[0])
@@ -314,16 +178,16 @@ def test_unrolling_parameterized_composite_gates(self):
         qc.append(subqc.to_instruction(), [qr2[2], qr2[3]])
 
         dag = circuit_to_dag(qc)
-        out_dag = Unroller(['u1', 'cx']).run(dag)
+        out_dag = Unroller(['u1', 'u3', 'cx']).run(dag)
 
         expected = QuantumCircuit(qr2)
         expected.u1(theta, qr2[0])
-        expected.u1(theta, qr2[2])
         expected.cx(qr2[0], qr2[1])
-        expected.cx(qr2[2], qr2[3])
         expected.u1(theta, qr2[1])
-        expected.u1(theta, qr2[3])
 
+        expected.u1(theta, qr2[2])
+        expected.cx(qr2[2], qr2[3])
+        expected.u1(theta, qr2[3])
         self.assertEqual(circuit_to_dag(expected), out_dag)
 
         # Expanding across register with shared parameter
@@ -336,14 +200,275 @@ def test_unrolling_parameterized_composite_gates(self):
         qc.append(subqc.to_instruction({theta: gamma}), [qr2[2], qr2[3]])
 
         dag = circuit_to_dag(qc)
-        out_dag = Unroller(['u1', 'cx']).run(dag)
+        out_dag = Unroller(['u1', 'u3', 'cx']).run(dag)
 
         expected = QuantumCircuit(qr2)
         expected.u1(phi, qr2[0])
-        expected.u1(gamma, qr2[2])
         expected.cx(qr2[0], qr2[1])
-        expected.cx(qr2[2], qr2[3])
         expected.u1(phi, qr2[1])
+
+        expected.u1(gamma, qr2[2])
+        expected.cx(qr2[2], qr2[3])
         expected.u1(gamma, qr2[3])
 
         self.assertEqual(circuit_to_dag(expected), out_dag)
+
+
+class TestUnrollAllInstructions(QiskitTestCase):
+    """Test unrolling a circuit containing all standard instructions."""
+
+    def setUp(self):
+        qr = self.qr = QuantumRegister(3, 'qr')
+        cr = self.cr = ClassicalRegister(3, 'cr')
+        self.circuit = QuantumCircuit(qr, cr)
+        self.ref_circuit = QuantumCircuit(qr, cr)
+        self.pass_ = Unroller(basis=['u3', 'cx', 'id'])
+
+    def compare_dags(self):
+        """compare dags in class tests"""
+        dag = circuit_to_dag(self.circuit)
+        unrolled_dag = self.pass_.run(dag)
+        ref_dag = circuit_to_dag(self.ref_circuit)
+        self.assertEqual(unrolled_dag, ref_dag)
+
+    def test_unroll_crx(self):
+        """test unroll crx"""
+        self.circuit.crx(0.5, 1, 2)
+        self.ref_circuit.u3(0, 0, pi/2, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(-0.25, 0, 0, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0.25, -pi/2, 0, 2)
+        self.compare_dags()
+
+    def test_unroll_cry(self):
+        """test unroll cry"""
+        self.circuit.cry(0.5, 1, 2)
+        self.ref_circuit.u3(0.25, 0, 0, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(-0.25, 0, 0, 2)
+        self.ref_circuit.cx(1, 2)
+        self.compare_dags()
+
+    def test_unroll_ccx(self):
+        """test unroll ccx"""
+        self.circuit.ccx(0, 1, 2)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0, 0, -pi/4, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0, 0, pi/4, 1)
+        self.ref_circuit.u3(0, 0, -pi/4, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.cx(0, 1)
+        self.ref_circuit.u3(0, 0, pi/4, 0)
+        self.ref_circuit.u3(0, 0, -pi/4, 1)
+        self.ref_circuit.cx(0, 1)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.compare_dags()
+
+    def test_unroll_ch(self):
+        """test unroll ch"""
+        self.circuit.ch(0, 2)
+        self.ref_circuit.u3(0, 0, pi/2, 2)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, -pi/4, 2)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.ref_circuit.u3(0, 0, -pi/2, 2)
+        self.compare_dags()
+
+    def test_unroll_crz(self):
+        """test unroll crz"""
+        self.circuit.crz(0.5, 1, 2)
+        self.ref_circuit.u3(0, 0, 0.25, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0, 0, -0.25, 2)
+        self.ref_circuit.cx(1, 2)
+
+    def test_unroll_cswap(self):
+        """test unroll cswap"""
+        self.circuit.cswap(1, 0, 2)
+        self.ref_circuit.cx(2, 0)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, -pi/4, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, pi/4, 0)
+        self.ref_circuit.u3(0, 0, -pi/4, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.cx(1, 0)
+        self.ref_circuit.u3(0, 0, -pi/4, 0)
+        self.ref_circuit.u3(0, 0, pi/4, 1)
+        self.ref_circuit.cx(1, 0)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.ref_circuit.u3(pi/2, 0, pi, 2)
+        self.ref_circuit.cx(2, 0)
+        self.compare_dags()
+
+    def test_unroll_cu1(self):
+        """test unroll cu1"""
+        self.circuit.cu1(0.1, 0, 2)
+        self.ref_circuit.u3(0, 0, 0.05, 0)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, -0.05, 2)
+        self.ref_circuit.cx(0, 2)
+        self.ref_circuit.u3(0, 0, 0.05, 2)
+        self.compare_dags()
+
+    def test_unroll_cu3(self):
+        """test unroll cu3"""
+        self.circuit.cu3(0.2, 0.1, 0.0, 1, 2)
+        self.ref_circuit.u3(0, 0, 0.05, 1)
+        self.ref_circuit.u3(0, 0, -0.05, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(-0.1, 0, -0.05, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0.1, 0.1, 0, 2)
+        self.compare_dags()
+
+    def test_unroll_cx(self):
+        """test unroll cx"""
+        self.circuit.cx(1, 0)
+        self.ref_circuit.cx(1, 0)
+        self.compare_dags()
+
+    def test_unroll_cy(self):
+        """test unroll cy"""
+        self.circuit.cy(1, 2)
+        self.ref_circuit.u3(0, 0, -pi/2, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.u3(0, 0, pi/2, 2)
+        self.compare_dags()
+
+    def test_unroll_cz(self):
+        """test unroll cz"""
+        self.circuit.cz(2, 0)
+        self.ref_circuit.u3(pi/2, 0, pi, 0)
+        self.ref_circuit.cx(2, 0)
+        self.ref_circuit.u3(pi/2, 0, pi, 0)
+        self.compare_dags()
+
+    def test_unroll_h(self):
+        """test unroll h"""
+        self.circuit.h(1)
+        self.ref_circuit.u3(pi/2, 0, pi, 1)
+        self.compare_dags()
+
+    def test_unroll_i(self):
+        """test unroll i"""
+        self.circuit.i(0)
+        self.ref_circuit.i(0)
+        self.compare_dags()
+
+    def test_unroll_rx(self):
+        """test unroll rx"""
+        self.circuit.rx(0.1, 0)
+        self.ref_circuit.u3(0.1, -pi/2, pi/2, 0)
+        self.compare_dags()
+
+    def test_unroll_ry(self):
+        """test unroll ry"""
+        self.circuit.ry(0.2, 1)
+        self.ref_circuit.u3(0.2, 0, 0, 1)
+        self.compare_dags()
+
+    def test_unroll_rz(self):
+        """test unroll rz"""
+        self.circuit.rz(0.3, 2)
+        self.ref_circuit.u3(0, 0, 0.3, 2)
+        self.compare_dags()
+
+    def test_unroll_rzz(self):
+        """test unroll rzz"""
+        self.circuit.rzz(0.6, 1, 0)
+        self.ref_circuit.cx(1, 0)
+        self.ref_circuit.u3(0, 0, 0.6, 0)
+        self.ref_circuit.cx(1, 0)
+        self.compare_dags()
+
+    def test_unroll_s(self):
+        """test unroll s"""
+        self.circuit.s(0)
+        self.ref_circuit.u3(0, 0, pi/2, 0)
+        self.compare_dags()
+
+    def test_unroll_sdg(self):
+        """test unroll sdg"""
+        self.circuit.sdg(1)
+        self.ref_circuit.u3(0, 0, -pi/2, 1)
+        self.compare_dags()
+
+    def test_unroll_swap(self):
+        """test unroll swap"""
+        self.circuit.swap(1, 2)
+        self.ref_circuit.cx(1, 2)
+        self.ref_circuit.cx(2, 1)
+        self.ref_circuit.cx(1, 2)
+        self.compare_dags()
+
+    def test_unroll_t(self):
+        """test unroll t"""
+        self.circuit.t(2)
+        self.ref_circuit.u3(0, 0, pi/4, 2)
+        self.compare_dags()
+
+    def test_unroll_tdg(self):
+        """test unroll tdg"""
+        self.circuit.tdg(0)
+        self.ref_circuit.u3(0, 0, -pi/4, 0)
+        self.compare_dags()
+
+    def test_unroll_u1(self):
+        """test unroll u1"""
+        self.circuit.u1(0.1, 1)
+        self.ref_circuit.u3(0, 0, 0.1, 1)
+        self.compare_dags()
+
+    def test_unroll_u2(self):
+        """test unroll u2"""
+        self.circuit.u2(0.2, -0.1, 0)
+        self.ref_circuit.u3(pi/2, 0.2, -0.1, 0)
+        self.compare_dags()
+
+    def test_unroll_u3(self):
+        """test unroll u3"""
+        self.circuit.u3(0.3, 0.0, -0.1, 2)
+        self.ref_circuit.u3(0.3, 0.0, -0.1, 2)
+        self.compare_dags()
+
+    def test_unroll_x(self):
+        """test unroll x"""
+        self.circuit.x(2)
+        self.ref_circuit.u3(pi, 0, pi, 2)
+        self.compare_dags()
+
+    def test_unroll_y(self):
+        """test unroll y"""
+        self.circuit.y(1)
+        self.ref_circuit.u3(pi, pi/2, pi/2, 1)
+        self.compare_dags()
+
+    def test_unroll_z(self):
+        """test unroll z"""
+        self.circuit.z(0)
+        self.ref_circuit.u3(0, 0, pi, 0)
+        self.compare_dags()
+
+    def test_unroll_snapshot(self):
+        """test unroll snapshot"""
+        self.circuit.snapshot('0')
+        self.ref_circuit.snapshot('0')
+        self.compare_dags()
+
+    def test_unroll_measure(self):
+        """test unroll measure"""
+        self.circuit.measure(self.qr, self.cr)
+        self.ref_circuit.measure(self.qr, self.cr)
+        self.compare_dags()