Repository: PennyLaneAI/pennylane
Pull request #5069 (2024-01-15 21:54:44)

Make parametric controlled operators inherit from `ControlledOp` 
**Context:**

Many controlled operations have been added to PennyLane over the years, such as CNOT, CRY, CZ, etc. For code tidiness, we would like to make them all inherit from the [ControlledOp](https://docs.pennylane.ai/en/stable/code/api/pennylane.ops.op_math.ControlledOp.html) base.

**Description of the Change:**

`CRX`, `CRY`, `CRZ`, `CRot`, `ControlledPhaseShift` are now:

1. Moved from `pennylane.ops.qubit` to `pennylane.ops.op_math`
2. Inherit from `ControlledOp`

Calling `qml.ctrl` on `RX`, `RY`, `RZ`, `Rot`, and `PhaseShift` with a single control wire will return gates of types `CRX`, `CRY`, etc. as opposed to a general `Controlled` operator.

The decomposition of an operator created with calling `qml.ctrl` on `RX`, `RY`, `RZ`, `Rot`, `PhaseShift` with a single control wire will now be the full decomposition instead of a single controlled gate. For example:
  ```
  >>> qml.ctrl(qml.RX(0.123, wires=1), control=0).decomposition()
  [CRX(0.123, wires=[0, 1])]
  ```
  Will become:
  ```
  >>> qml.ctrl(qml.RX(0.123, wires=1), control=0).decomposition()
  [
    RZ(1.5707963267948966, wires=[1]), 
    RY(0.0615, wires=[1]), 
    CNOT(wires=[0, 1]), 
    RY(-0.0615, wires=[1]), 
    CNOT(wires=[0, 1]), 
    RZ(-1.5707963267948966, wires=[1])
  ]
  ```

**Benefits:**

Code cleanliness.

**Related GitHub Issues:**

https://github.com/PennyLaneAI/pennylane/issues/1447


github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
trbromley: Given the size of the diff, could it make sense to split this up into smaller PRs?
astralcai: > Given the size of the diff, could it make sense to split this up into smaller PRs?

The size of the diff is actually smaller than it looks. We're moving these controlled operators from `ops.qubit` to `ops.op_math`, which is probably why it looks like a lot.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: `9 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`c01cb4b`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/c01cb4b96f758a9f1a0c4e422115829b80e113cb?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.67% compared to head [(`446090e`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) 99.62%.
> Report is 20 commits behind head on master.

| [Files](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [pennylane/ops/op\_math/controlled\_parametric\_ops.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL29wcy9vcF9tYXRoL2NvbnRyb2xsZWRfcGFyYW1ldHJpY19vcHMucHk=) | 94.96% | [7 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |
| [pennylane/workflow/jacobian\_products.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3dvcmtmbG93L2phY29iaWFuX3Byb2R1Y3RzLnB5) | 91.30% | [2 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #5069      +/-   ##
==========================================
- Coverage   99.67%   99.62%   -0.06%     
==========================================
  Files         394      394              
  Lines       35670    35432     -238     
==========================================
- Hits        35554    35298     -256     
- Misses        116      134      +18     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/5069?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

astralcai: codecov seems to be failing again for all branches of `if interface == "tensorflow"` in `controlled_ops.py`. 
astralcai: > So my general question for all these operations is:
> 
> Which methods and properties can we simply inherit from `Controlled` and still get the same behavior?
> 
> * `label`
> * `basis`
> * `grad_method`
> * `compute_matrix`? This might be based on performance considerations.
> * `decomposition`/ `compute_decomposition`
> * `compute_eigvals`
> 
> `parameter_frequencies` is probably worth hard-coding because otherwise it can get a bit expensive to compute.

`label`, `basis`, `grad_method` are quite straight forward to remove from the sub-classes, but the rest seems a bit complicated. The main source of issue is that `Controlled` does not define static methods `compute_matrix` and `compute_decomposition`, instead it implements `matrix` and `decomposition` directly. So here are some observations:

1. Currently all controlled operators define their own `compute_matrix` but does not override `matrix`, so if `qml.CY.compute_matrix()` is called, we do get the hard-coded matrix, but if we do `qml.CY(0).matrix()`, the method resolves to `Controlled.matrix`, which computes the matrix from scratch and never uses `CY.compute_matrix`
2. For rotational gates, their `compute_matrix` and `compute_decomposition` are parametrized functions, so if we want to keep supporting things like `qml.CRX.compute_decomposition(1.2, wires=(0,1))`, we would have to keep it in `CRX`. If we remove the `decomposition`/ `compute_decomposition` combo, `qml.CRX.compute_decomposition` is going to resolve to `Operator.compute_decomposition`, which is undefined. I don't see a clear path forward to solve this, since the method signature for `CRX.compute_decomposition`: `def compute_decomposition(phi, wires)` does not make sense for the general `Controlled` class.

I think cleaning this up would probably have to be a separate PR, as @mudit2812 mentioned.
Alex-Preciado: Linking to corresponding Shortcut story: [sc-53529]

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index dba94eba7..3231f7d40 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -85,6 +85,9 @@
 * Raise a more informative error when calling `adjoint_jacobian` with trainable state-prep operations.
   [(#5026)](https://github.com/PennyLaneAI/pennylane/pull/5026)
 
+* `CRX`, `CRY`, `CRZ`, `CROT`, and `ControlledPhaseShift` (i.e. `CPhaseShift`) now inherit from `ControlledOp`, giving them additional properties such as `control_wire` and `control_values`. Calling `qml.ctrl` on `RX`, `RY`, `RZ`, `Rot`, and `PhaseShift` with a single control wire will return gates of types `CRX`, `CRY`, etc. as opposed to a general `Controlled` operator.
+  [(#5069)](https://github.com/PennyLaneAI/pennylane/pull/5069)
+
 <h4>Community contributions 🥳</h4>
 
 * The transform `split_non_commuting` now accepts measurements of type `probs`, `sample` and `counts` which accept both wires and observables. 
@@ -113,6 +116,20 @@
   (with potentially negative eigenvalues) has been implemented.
   [(#5048)](https://github.com/PennyLaneAI/pennylane/pull/5048)
 
+* The decomposition of an operator created with calling `qml.ctrl` on a parametric operator (specifically `RX`, `RY`, `RZ`, `Rot`, `PhaseShift`) with a single control wire will now be the full decomposition instead of a single controlled gate. For example:
+  ```
+  >>> qml.ctrl(qml.RX(0.123, wires=1), control=0).decomposition()
+  [
+    RZ(1.5707963267948966, wires=[1]), 
+    RY(0.0615, wires=[1]), 
+    CNOT(wires=[0, 1]), 
+    RY(-0.0615, wires=[1]), 
+    CNOT(wires=[0, 1]), 
+    RZ(-1.5707963267948966, wires=[1])
+  ]
+  ```
+  [(#5069)](https://github.com/PennyLaneAI/pennylane/pull/5069)
+
 * `QuantumScript.is_sampled` and `QuantumScript.all_sampled` have been removed. Users should now
   validate these properties manually.
   [(#5072)](https://github.com/PennyLaneAI/pennylane/pull/5072)
diff --git a/pennylane/data/attributes/operator/operator.py b/pennylane/data/attributes/operator/operator.py
index 8dd84c105..e80b4b917 100644
--- a/pennylane/data/attributes/operator/operator.py
+++ b/pennylane/data/attributes/operator/operator.py
@@ -74,8 +74,6 @@ class DatasetOperator(Generic[Op], DatasetAttribute[HDF5Group, Op, Op]):
                 qml.S,
                 qml.SX,
                 qml.CNOT,
-                qml.CZ,
-                qml.CY,
                 qml.CH,
                 qml.SWAP,
                 qml.ECR,
@@ -87,15 +85,6 @@ class DatasetOperator(Generic[Op], DatasetAttribute[HDF5Group, Op, Op]):
                 # pennylane/ops/qubit/observables.py
                 qml.Hermitian,
                 qml.Projector,
-                # pennylane/ops/qubit/parametric_ops_controlled.py
-                qml.ControlledPhaseShift,
-                qml.CPhaseShift00,
-                qml.CPhaseShift01,
-                qml.CPhaseShift10,
-                qml.CRX,
-                qml.CRY,
-                qml.CRZ,
-                qml.CRot,
                 # pennylane/ops/qubit/parametric_ops_multi_qubit.py
                 qml.MultiRZ,
                 qml.IsingXX,
@@ -103,6 +92,9 @@ class DatasetOperator(Generic[Op], DatasetAttribute[HDF5Group, Op, Op]):
                 qml.IsingZZ,
                 qml.IsingXY,
                 qml.PSWAP,
+                qml.CPhaseShift00,
+                qml.CPhaseShift01,
+                qml.CPhaseShift10,
                 # pennylane/ops/qubit/parametric_ops_single_qubit.py
                 qml.RX,
                 qml.RY,
@@ -177,6 +169,15 @@ class DatasetOperator(Generic[Op], DatasetAttribute[HDF5Group, Op, Op]):
                 qml.FockStateProjector,
                 # pennylane/ops/identity.py
                 qml.Identity,
+                # pennylane/ops/op_math/controlled_ops.py
+                qml.ControlledQubitUnitary,
+                qml.ControlledPhaseShift,
+                qml.CRX,
+                qml.CRY,
+                qml.CRZ,
+                qml.CRot,
+                qml.CZ,
+                qml.CY,
             )
         )
 
diff --git a/pennylane/operation.py b/pennylane/operation.py
index c4404fea0..03c35ef57 100644
--- a/pennylane/operation.py
+++ b/pennylane/operation.py
@@ -397,8 +397,6 @@ def _process_data(op):
     # valued QNode (one that returns qml.state) requires complex typed inputs.
     if op.name in ("RX", "RY", "RZ", "PhaseShift", "Rot"):
         mod_val = 2 * np.pi
-    elif op.name in ("CRX", "CRY", "CRZ", "CRot"):
-        mod_val = 4 * np.pi
     else:
         mod_val = None
 
diff --git a/pennylane/ops/functions/bind_new_parameters.py b/pennylane/ops/functions/bind_new_parameters.py
index e976e43a6..96b93b29c 100644
--- a/pennylane/ops/functions/bind_new_parameters.py
+++ b/pennylane/ops/functions/bind_new_parameters.py
@@ -16,7 +16,7 @@ This module contains the qml.bind_new_parameters function.
 """
 # pylint: disable=missing-docstring
 
-from typing import Sequence
+from typing import Sequence, Union
 import copy
 from functools import singledispatch
 
@@ -111,18 +111,24 @@ def bind_new_parameters_composite_op(op: CompositeOp, params: Sequence[TensorLik
     return op.__class__(*new_operands)
 
 
-@bind_new_parameters.register
-def bind_new_parameters_cy(
-    op: qml.CY, params: Sequence[TensorLike]
+@bind_new_parameters.register(qml.CY)
+@bind_new_parameters.register(qml.CZ)
+def bind_new_parameters_copy(
+    op: Union[qml.CY, qml.CZ], params: Sequence[TensorLike]
 ):  # pylint:disable=unused-argument
     return copy.copy(op)
 
 
-@bind_new_parameters.register
-def bind_new_parameters_cz(
-    op: qml.CZ, params: Sequence[TensorLike]
-):  # pylint:disable=unused-argument
-    return copy.copy(op)
+@bind_new_parameters.register(qml.CRX)
+@bind_new_parameters.register(qml.CRY)
+@bind_new_parameters.register(qml.CRZ)
+@bind_new_parameters.register(qml.CRot)
+@bind_new_parameters.register(qml.ControlledPhaseShift)
+def bind_new_parameters_parametric_controlled_ops(
+    op: Union[qml.CRX, qml.CRY, qml.CRZ, qml.CRot, qml.ControlledPhaseShift],
+    params: Sequence[TensorLike],
+):
+    return op.__class__(*params, wires=op.wires)
 
 
 @bind_new_parameters.register
diff --git a/pennylane/ops/op_math/__init__.py b/pennylane/ops/op_math/__init__.py
index 6f57455e6..e15ab1cf6 100644
--- a/pennylane/ops/op_math/__init__.py
+++ b/pennylane/ops/op_math/__init__.py
@@ -64,6 +64,11 @@ Controlled Operator Classes
     ~ControlledQubitUnitary
     ~CY
     ~CZ
+    ~CRX
+    ~CRY
+    ~CRZ
+    ~CRot
+    ~ControlledPhaseShift
 
 Decompositions
 ~~~~~~~~~~~~~~
@@ -94,7 +99,17 @@ from .adjoint import Adjoint, adjoint
 from .composite import CompositeOp
 from .condition import cond, Conditional
 from .controlled import Controlled, ControlledOp, ctrl
-from .controlled_ops import ControlledQubitUnitary, CY, CZ
+from .controlled_ops import (
+    ControlledQubitUnitary,
+    ControlledPhaseShift,
+    CPhase,
+    CRot,
+    CRX,
+    CRY,
+    CRZ,
+    CY,
+    CZ,
+)
 from .decompositions import one_qubit_decomposition, two_qubit_decomposition, sk_decomposition
 from .evolution import Evolution
 from .exp import Exp, exp
@@ -105,4 +120,14 @@ from .sum import Sum, sum
 from .symbolicop import ScalarSymbolicOp, SymbolicOp
 from .controlled_decompositions import ctrl_decomp_zyz, ctrl_decomp_bisect
 
-controlled_qubit_ops = {"ControlledQubitUnitary", "CY", "CZ"}
+controlled_qubit_ops = {
+    "ControlledQubitUnitary",
+    "CY",
+    "CZ",
+    "CRX",
+    "CRY",
+    "CRZ",
+    "CRot",
+    "ControlledPhaseShift",
+    "CPhase",
+}
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index 766a683a7..ccef2b9dc 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -127,19 +127,24 @@ def ctrl(op, control, control_values=None, work_wires=None):
         ops_loader = available_eps[active_jit]["ops"].load()
         return ops_loader.ctrl(op, control, control_values=control_values, work_wires=work_wires)
 
-    custom_controlled_ops = {
+    custom_ops = {
         (qml.PauliZ, 1): qml.CZ,
         (qml.PauliY, 1): qml.CY,
         (qml.PauliX, 1): qml.CNOT,
         (qml.PauliX, 2): qml.Toffoli,
+        (qml.RX, 1): qml.CRX,
+        (qml.RY, 1): qml.CRY,
+        (qml.RZ, 1): qml.CRZ,
+        (qml.Rot, 1): qml.CRot,
+        (qml.PhaseShift, 1): qml.ControlledPhaseShift,
     }
     control_values = [control_values] if isinstance(control_values, (int, bool)) else control_values
     control = qml.wires.Wires(control)
     custom_key = (type(op), len(control))
 
-    if custom_key in custom_controlled_ops and (control_values is None or all(control_values)):
+    if custom_key in custom_ops and (control_values is None or all(control_values)):
         qml.QueuingManager.remove(op)
-        return custom_controlled_ops[custom_key](control + op.wires)
+        return custom_ops[custom_key](*op.data, control + op.wires)
     if isinstance(op, qml.PauliX):
         qml.QueuingManager.remove(op)
         control_string = (
@@ -338,7 +343,12 @@ class Controlled(SymbolicOp):
         # these gates do not consider global phases in their hash
         if self.base.name in ("RX", "RY", "RZ", "Rot"):
             base_params = str(
-                [qml.math.round(qml.math.real(d) % (4 * np.pi), 10) for d in self.base.data]
+                [
+                    id(d)
+                    if qml.math.is_abstract(d)
+                    else qml.math.round(qml.math.real(d) % (4 * np.pi), 10)
+                    for d in self.base.data
+                ]
             )
             base_hash = hash(
                 (
@@ -589,7 +599,7 @@ class Controlled(SymbolicOp):
             for op in base_pow
         ]
 
-    def simplify(self) -> "Controlled":
+    def simplify(self) -> "Operator":
         if isinstance(self.base, Controlled):
             base = self.base.base.simplify()
             return ctrl(
@@ -599,8 +609,12 @@ class Controlled(SymbolicOp):
                 work_wires=self.work_wires + self.base.work_wires,
             )
 
+        simplified_base = self.base.simplify()
+        if isinstance(simplified_base, qml.Identity):
+            return simplified_base
+
         return ctrl(
-            op=self.base.simplify(),
+            op=simplified_base,
             control=self.control_wires,
             control_values=self.control_values,
             work_wires=self.work_wires,
diff --git a/pennylane/ops/op_math/controlled_ops.py b/pennylane/ops/op_math/controlled_ops.py
index b9b8c601f..22844dc64 100644
--- a/pennylane/ops/op_math/controlled_ops.py
+++ b/pennylane/ops/op_math/controlled_ops.py
@@ -1,4 +1,4 @@
-# Copyright 2018-2022 Xanadu Quantum Technologies Inc.
+# Copyright 2018-2024 Xanadu Quantum Technologies Inc.
 
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -12,7 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 """
-This submodule contains controlled operators based on the Controlled and ControlledOp class.
+This submodule contains controlled operators based on the ControlledOp class.
 """
 
 import warnings
@@ -23,8 +23,8 @@ import numpy as np
 
 import pennylane as qml
 from pennylane.operation import AnyWires
-from pennylane.ops.qubit.non_parametric_ops import PauliY, S, PauliZ, CCZ
 from pennylane.ops.qubit.matrix_ops import QubitUnitary
+from pennylane.ops.qubit.parametric_ops_single_qubit import stack_last
 from .controlled import ControlledOp
 
 
@@ -39,7 +39,8 @@ class ControlledQubitUnitary(ControlledOp):
     * ``control_wires``: wires that act as control for the operation
     * ``control_values``: the state on which to apply the controlled operation (see below)
     * ``target_wires``: the wires the unitary matrix will be applied to
-    * ``active_wires``: Wires modified by the operator. This is the control wires followed by the target wires.
+    * ``active_wires``: Wires modified by the operator. This is the control wires followed
+        by the target wires.
 
     **Details:**
 
@@ -49,13 +50,17 @@ class ControlledQubitUnitary(ControlledOp):
     * Gradient recipe: None
 
     Args:
-        base (Union[array[complex], QubitUnitary]): square unitary matrix or a QubitUnitary operation. If passing a matrix,
-          this will be used to construct a QubitUnitary operator that will be used as the base operator. If providing
-          a ``qml.QubitUnitary``, this will be used as the base directly.
+        base (Union[array[complex], QubitUnitary]): square unitary matrix or a QubitUnitary
+            operation. If passing a matrix, this will be used to construct a QubitUnitary
+            operator that will be used as the base operator. If providing a ``qml.QubitUnitary``,
+            this will be used as the base directly.
         control_wires (Union[Wires, Sequence[int], or int]): the control wire(s)
-        wires (Union[Wires, Sequence[int], or int]): the wire(s) the unitary acts on (optional if U is provided as a QubitUnitary)
-        control_values (List[int, bool]): a list providing the state of the control qubits to control on (default is the all 1s state)
-        unitary_check (bool): whether to check whether an array U is unitary when creating the operator (default False)
+        wires (Union[Wires, Sequence[int], or int]): the wire(s) the unitary acts on
+            (optional if U is provided as a QubitUnitary)
+        control_values (List[int, bool]): a list providing the state of the control qubits to
+            control on (default is the all 1s state)
+        unitary_check (bool): whether to check whether an array U is unitary when creating the
+            operator (default False)
 
     **Example**
 
@@ -190,8 +195,7 @@ class CY(ControlledOp):
     ndim_params = ()
     """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
 
-    grad_method = None
-    """Gradient computation method."""
+    name = "CY"
 
     def _flatten(self):
         return tuple(), (self.wires,)
@@ -202,10 +206,7 @@ class CY(ControlledOp):
 
     def __init__(self, wires, id=None):
         control_wire, wire = wires
-        base = PauliY(wire)
-
-        super().__init__(base, control_wire, id=id)
-        self._name = "CY"
+        super().__init__(qml.PauliY(wire), control_wire, id=id)
 
     @staticmethod
     @lru_cache()
@@ -260,7 +261,7 @@ class CY(ControlledOp):
         [CRY(3.141592653589793, wires=[0, 1]), S(wires=[0])]
 
         """
-        return [qml.CRY(np.pi, wires=wires), S(wires=wires[0])]
+        return [qml.CRY(np.pi, wires=wires), qml.S(wires=wires[0])]
 
     def decomposition(self):
         return self.compute_decomposition(self.wires)
@@ -296,6 +297,8 @@ class CZ(ControlledOp):
     ndim_params = ()
     """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
 
+    name = "CZ"
+
     def _flatten(self):
         return tuple(), (self.wires,)
 
@@ -305,9 +308,7 @@ class CZ(ControlledOp):
 
     def __init__(self, wires):
         control_wire, wire = wires
-        base = PauliZ(wires=wire)
-        super().__init__(base, control_wire)
-        self._name = "CZ"
+        super().__init__(qml.PauliZ(wires=wire), control_wire)
 
     @staticmethod
     @lru_cache()
@@ -333,4 +334,871 @@ class CZ(ControlledOp):
         return np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])
 
     def _controlled(self, wire):
-        return CCZ(wires=wire + self.wires)
+        return qml.CCZ(wires=wire + self.wires)
+
+
+class CRX(ControlledOp):
+    r"""The controlled-RX operator
+
+    .. math::
+
+        \begin{align}
+            CR_x(\phi) &=
+            \begin{bmatrix}
+            & 1 & 0 & 0 & 0 \\
+            & 0 & 1 & 0 & 0\\
+            & 0 & 0 & \cos(\phi/2) & -i\sin(\phi/2)\\
+            & 0 & 0 & -i\sin(\phi/2) & \cos(\phi/2)
+            \end{bmatrix}.
+        \end{align}
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe: The controlled-RX operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
+
+      .. math::
+
+          \frac{d}{d\phi}f(CR_x(\phi)) = c_+ \left[f(CR_x(\phi+a)) - f(CR_x(\phi-a))\right] - c_- \left[f(CR_x(\phi+b)) - f(CR_x(\phi-b))\right]
+
+      where :math:`f` is an expectation value depending on :math:`CR_x(\phi)`, and
+
+      - :math:`a = \pi/2`
+      - :math:`b = 3\pi/2`
+      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    """int: Number of wires that the operation acts on."""
+
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "CRX"
+    parameter_frequencies = [(0.5, 1.0)]
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(qml.RX(phi, wires=wires[1:]), control_wires=wires[0], id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        base = data[0]
+        control_wires = metadata[0]
+        return cls(*base.data, wires=control_wires + base.wires)
+
+    @staticmethod
+    def compute_matrix(theta):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CRX.matrix`
+
+        Args:
+            theta (tensor_like or float): rotation angle
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CRX.compute_matrix(torch.tensor(0.5))
+        tensor([[1.0+0.0j, 0.0+0.0j,    0.0+0.0j,    0.0+0.0j],
+                [0.0+0.0j, 1.0+0.0j,    0.0+0.0j,    0.0+0.0j],
+                [0.0+0.0j, 0.0+0.0j, 0.9689+0.0j, 0.0-0.2474j],
+                [0.0+0.0j, 0.0+0.0j, 0.0-0.2474j, 0.9689+0.0j]])
+        """
+        interface = qml.math.get_interface(theta)
+
+        c = qml.math.cos(theta / 2)
+        s = qml.math.sin(theta / 2)
+
+        if interface == "tensorflow":
+            c = qml.math.cast_like(c, 1j)
+            s = qml.math.cast_like(s, 1j)
+
+        # The following avoids casting an imaginary quantity to reals when back propagating
+        c = (1 + 0j) * c
+        js = -1j * s
+        ones = qml.math.ones_like(js)
+        zeros = qml.math.zeros_like(js)
+        matrix = [
+            [ones, zeros, zeros, zeros],
+            [zeros, ones, zeros, zeros],
+            [zeros, zeros, c, js],
+            [zeros, zeros, js, c],
+        ]
+
+        return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+
+    @staticmethod
+    def compute_decomposition(phi, wires):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+        .. seealso:: :meth:`~.CRot.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): the wires the operation acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CRX.compute_decomposition(1.2, wires=(0,1))
+        [RZ(1.5707963267948966, wires=[1]),
+        RY(0.6, wires=[1]),
+        CNOT(wires=[0, 1]),
+        RY(-0.6, wires=[1]),
+        CNOT(wires=[0, 1]),
+        RZ(-1.5707963267948966, wires=[1])]
+
+        """
+        pi_half = qml.math.ones_like(phi) * (np.pi / 2)
+        return [
+            qml.RZ(pi_half, wires=wires[1]),
+            qml.RY(phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.RY(-phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.RZ(-pi_half, wires=wires[1]),
+        ]
+
+    def decomposition(self):
+        r"""Representation of the operator as a product of other operators.
+
+        .. math:: O = O_1 O_2 \dots O_n
+
+        A ``DecompositionUndefinedError`` is raised if no representation by decomposition is defined.
+
+        .. seealso:: :meth:`~.Operator.compute_decomposition`.
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
+
+
+class CRY(ControlledOp):
+    r"""The controlled-RY operator
+
+    .. math::
+
+        \begin{align}
+            CR_y(\phi) &=
+            \begin{bmatrix}
+                1 & 0 & 0 & 0 \\
+                0 & 1 & 0 & 0\\
+                0 & 0 & \cos(\phi/2) & -\sin(\phi/2)\\
+                0 & 0 & \sin(\phi/2) & \cos(\phi/2)
+            \end{bmatrix}.
+        \end{align}
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe: The controlled-RY operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
+
+      .. math::
+
+          \frac{d}{d\phi}f(CR_y(\phi)) = c_+ \left[f(CR_y(\phi+a)) - f(CR_y(\phi-a))\right] - c_- \left[f(CR_y(\phi+b)) - f(CR_y(\phi-b))\right]
+
+      where :math:`f` is an expectation value depending on :math:`CR_y(\phi)`, and
+
+      - :math:`a = \pi/2`
+      - :math:`b = 3\pi/2`
+      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    """int: Number of wires that the operation acts on."""
+
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "CRY"
+    parameter_frequencies = [(0.5, 1.0)]
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(qml.RY(phi, wires=wires[1:]), control_wires=wires[0], id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        base = data[0]
+        control_wires = metadata[0]
+        return cls(*base.data, wires=control_wires + base.wires)
+
+    @staticmethod
+    def compute_matrix(theta):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CRY.matrix`
+
+
+        Args:
+            theta (tensor_like or float): rotation angle
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CRY.compute_matrix(torch.tensor(0.5))
+        tensor([[ 1.0000,  0.0000,  0.0000,  0.0000],
+                [ 0.0000,  1.0000,  0.0000,  0.0000],
+                [ 0.0000,  0.0000,  0.9689, -0.2474],
+                [ 0.0000,  0.0000,  0.2474,  0.9689]], dtype=torch.float64)
+        """
+        interface = qml.math.get_interface(theta)
+
+        c = qml.math.cos(theta / 2)
+        s = qml.math.sin(theta / 2)
+
+        if interface == "tensorflow":
+            c = qml.math.cast_like(c, 1j)
+            s = qml.math.cast_like(s, 1j)
+
+        # The following avoids casting an imaginary quantity to reals when back propagating
+        c = (1 + 0j) * c
+        s = (1 + 0j) * s
+        ones = qml.math.ones_like(s)
+        zeros = qml.math.zeros_like(s)
+        matrix = [
+            [ones, zeros, zeros, zeros],
+            [zeros, ones, zeros, zeros],
+            [zeros, zeros, c, -s],
+            [zeros, zeros, s, c],
+        ]
+
+        return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+
+    @staticmethod
+    def compute_decomposition(phi, wires):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+        .. seealso:: :meth:`~.CRY.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CRY.compute_decomposition(1.2, wires=(0,1))
+        [RY(0.6, wires=[1]),
+        CNOT(wires=[0, 1]),
+        RY(-0.6, wires=[1]),
+        CNOT(wires=[0, 1])]
+
+        """
+        return [
+            qml.RY(phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.RY(-phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+        ]
+
+    def decomposition(self):
+        r"""Representation of the operator as a product of other operators.
+
+        .. math:: O = O_1 O_2 \dots O_n
+
+        A ``DecompositionUndefinedError`` is raised if no representation by decomposition is defined.
+
+        .. seealso:: :meth:`~.Operator.compute_decomposition`.
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
+
+
+class CRZ(ControlledOp):
+    r"""The controlled-RZ operator
+
+    .. math::
+
+        \begin{align}
+             CR_z(\phi) &=
+             \begin{bmatrix}
+                1 & 0 & 0 & 0 \\
+                0 & 1 & 0 & 0\\
+                0 & 0 & e^{-i\phi/2} & 0\\
+                0 & 0 & 0 & e^{i\phi/2}
+            \end{bmatrix}.
+        \end{align}
+
+
+    .. note:: The subscripts of the operations in the formula refer to the wires they act on, e.g. 1 corresponds
+        to the first element in ``wires`` that is the **control qubit**.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe: The controlled-RZ operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
+
+      .. math::
+
+          \frac{d}{d\phi}f(CR_z(\phi)) = c_+ \left[f(CR_z(\phi+a)) - f(CR_z(\phi-a))\right] - c_- \left[f(CR_z(\phi+b)) - f(CR_z(\phi-b))\right]
+
+      where :math:`f` is an expectation value depending on :math:`CR_z(\phi)`, and
+
+      - :math:`a = \pi/2`
+      - :math:`b = 3\pi/2`
+      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+
+    """
+    num_wires = 2
+    """int: Number of wires that the operation acts on."""
+
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "CRZ"
+    parameter_frequencies = [(0.5, 1.0)]
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(qml.RZ(phi, wires=wires[1:]), control_wires=wires[0], id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        base = data[0]
+        control_wires = metadata[0]
+        return cls(*base.data, wires=control_wires + base.wires)
+
+    @staticmethod
+    def compute_matrix(theta):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CRZ.matrix`
+
+        Args:
+            theta (tensor_like or float): rotation angle
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CRZ.compute_matrix(torch.tensor(0.5))
+        tensor([[1.0+0.0j, 0.0+0.0j,       0.0+0.0j,       0.0+0.0j],
+                [0.0+0.0j, 1.0+0.0j,       0.0+0.0j,       0.0+0.0j],
+                [0.0+0.0j, 0.0+0.0j, 0.9689-0.2474j,       0.0+0.0j],
+                [0.0+0.0j, 0.0+0.0j,       0.0+0.0j, 0.9689+0.2474j]])
+        """
+        if qml.math.get_interface(theta) == "tensorflow":
+            p = qml.math.exp(-0.5j * qml.math.cast_like(theta, 1j))
+            if qml.math.ndim(p) == 0:
+                return qml.math.diag([1, 1, p, qml.math.conj(p)])
+
+            ones = qml.math.ones_like(p)
+            diags = stack_last([ones, ones, p, qml.math.conj(p)])
+            return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+
+        signs = qml.math.array([0, 0, 1, -1], like=theta)
+        arg = -0.5j * theta
+
+        if qml.math.ndim(arg) == 0:
+            return qml.math.diag(qml.math.exp(arg * signs))
+
+        diags = qml.math.exp(qml.math.outer(arg, signs))
+        return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+
+    @staticmethod
+    def compute_eigvals(theta, **_):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.CRZ.eigvals`
+
+
+        Args:
+            theta (tensor_like or float): rotation angle
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.CRZ.compute_eigvals(torch.tensor(0.5))
+        tensor([1.0000+0.0000j, 1.0000+0.0000j, 0.9689-0.2474j, 0.9689+0.2474j])
+        """
+        if qml.math.get_interface(theta) == "tensorflow":
+            phase = qml.math.exp(-0.5j * qml.math.cast_like(theta, 1j))
+            ones = qml.math.ones_like(phase)
+            return stack_last([ones, ones, phase, qml.math.conj(phase)])
+
+        prefactors = qml.math.array([0, 0, -0.5j, 0.5j], like=theta)
+        if qml.math.ndim(theta) == 0:
+            product = theta * prefactors
+        else:
+            product = qml.math.outer(theta, prefactors)
+        return qml.math.exp(product)
+
+    def eigvals(self):
+        return self.compute_eigvals(*self.parameters)
+
+    @staticmethod
+    def compute_decomposition(phi, wires):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+        .. seealso:: :meth:`~.CRZ.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CRZ.compute_decomposition(1.2, wires=(0,1))
+        [PhaseShift(0.6, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(-0.6, wires=[1]),
+        CNOT(wires=[0, 1])]
+
+        """
+        return [
+            qml.PhaseShift(phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.PhaseShift(-phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+        ]
+
+    def decomposition(self):
+        r"""Representation of the operator as a product of other operators.
+
+        .. math:: O = O_1 O_2 \dots O_n
+
+        A ``DecompositionUndefinedError`` is raised if no representation by decomposition is defined.
+
+        .. seealso:: :meth:`~.Operator.compute_decomposition`.
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
+
+
+class CRot(ControlledOp):
+    r"""The controlled-Rot operator
+
+    .. math:: CR(\phi, \theta, \omega) = \begin{bmatrix}
+            1 & 0 & 0 & 0 \\
+            0 & 1 & 0 & 0\\
+            0 & 0 & e^{-i(\phi+\omega)/2}\cos(\theta/2) & -e^{i(\phi-\omega)/2}\sin(\theta/2)\\
+            0 & 0 & e^{-i(\phi-\omega)/2}\sin(\theta/2) & e^{i(\phi+\omega)/2}\cos(\theta/2)
+        \end{bmatrix}.
+
+    .. note:: The first wire provided corresponds to the **control qubit**.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 3
+    * Number of dimensions per parameter: (0, 0, 0)
+    * Gradient recipe: The controlled-Rot operator satisfies a four-term parameter-shift rule
+      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
+
+      .. math::
+
+          \frac{d}{d\mathbf{x}_i}f(CR(\mathbf{x}_i)) = c_+ \left[f(CR(\mathbf{x}_i+a)) - f(CR(\mathbf{x}_i-a))\right] - c_- \left[f(CR(\mathbf{x}_i+b)) - f(CR(\mathbf{x}_i-b))\right]
+
+      where :math:`f` is an expectation value depending on :math:`CR(\mathbf{x}_i)`, and
+
+      - :math:`\mathbf{x} = (\phi, \theta, \omega)` and `i` is an index to :math:`\mathbf{x}`
+      - :math:`a = \pi/2`
+      - :math:`b = 3\pi/2`
+      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        theta (float): rotation angle :math:`\theta`
+        omega (float): rotation angle :math:`\omega`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+
+    """
+    num_wires = 2
+    """int: Number of wires this operator acts on."""
+
+    num_params = 3
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0, 0, 0)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "CRot"
+    parameter_frequencies = [(0.5, 1.0), (0.5, 1.0), (0.5, 1.0)]
+
+    def __init__(self, phi, theta, omega, wires, id=None):  # pylint: disable=too-many-arguments
+        super().__init__(qml.Rot(phi, theta, omega, wires=wires[1]), control_wires=wires[0], id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        base = data[0]
+        control_wires = metadata[0]
+        return cls(*base.data, wires=control_wires + base.wires)
+
+    @staticmethod
+    def compute_matrix(phi, theta, omega):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CRot.matrix`
+
+
+        Args:
+            phi(tensor_like or float): first rotation angle
+            theta (tensor_like or float): second rotation angle
+            omega (tensor_like or float): third rotation angle
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+         >>> qml.CRot.compute_matrix(torch.tensor(0.1), torch.tensor(0.2), torch.tensor(0.3))
+         tensor([[ 1.0+0.0j,  0.0+0.0j,        0.0+0.0j,        0.0+0.0j],
+                [ 0.0+0.0j,  1.0+0.0j,        0.0+0.0j,        0.0+0.0j],
+                [ 0.0+0.0j,  0.0+0.0j,  0.9752-0.1977j, -0.0993+0.0100j],
+                [ 0.0+0.0j,  0.0+0.0j,  0.0993+0.0100j,  0.9752+0.1977j]])
+        """
+        # It might be that they are in different interfaces, e.g.,
+        # CRot(0.2, 0.3, tf.Variable(0.5), wires=[0, 1])
+        # So we need to make sure the matrix comes out having the right type
+        interface = qml.math.get_interface(phi, theta, omega)
+
+        c = qml.math.cos(theta / 2)
+        s = qml.math.sin(theta / 2)
+
+        # If anything is not tensorflow, it has to be casted
+        if interface == "tensorflow":
+            phi = qml.math.cast_like(qml.math.asarray(phi, like=interface), 1j)
+            omega = qml.math.cast_like(qml.math.asarray(omega, like=interface), 1j)
+            c = qml.math.cast_like(qml.math.asarray(c, like=interface), 1j)
+            s = qml.math.cast_like(qml.math.asarray(s, like=interface), 1j)
+
+        # The following variable is used to assert the all terms to be stacked have same shape
+        one = qml.math.ones_like(phi) * qml.math.ones_like(omega)
+        c = c * one
+        s = s * one
+
+        o = qml.math.ones_like(c)
+        z = qml.math.zeros_like(c)
+        mat = [
+            [o, z, z, z],
+            [z, o, z, z],
+            [
+                z,
+                z,
+                qml.math.exp(-0.5j * (phi + omega)) * c,
+                -qml.math.exp(0.5j * (phi - omega)) * s,
+            ],
+            [
+                z,
+                z,
+                qml.math.exp(-0.5j * (phi - omega)) * s,
+                qml.math.exp(0.5j * (phi + omega)) * c,
+            ],
+        ]
+
+        return qml.math.stack([stack_last(row) for row in mat], axis=-2)
+
+    @staticmethod
+    def compute_decomposition(phi, theta, omega, wires):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+        .. seealso:: :meth:`~.CRot.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            theta (float): rotation angle :math:`\theta`
+            omega (float): rotation angle :math:`\omega`
+            wires (Iterable, Wires): the wires the operation acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.PhaseShift.compute_decomposition(1.234, wires=0)
+        [RZ(-1.1, wires=[1]),
+        CNOT(wires=[0, 1]),
+        RZ(-2.3, wires=[1]),
+        RY(-1.15, wires=[1]),
+        CNOT(wires=[0, 1]),
+        RY(1.15, wires=[1]),
+        RZ(3.4, wires=[1])]
+
+        """
+        return [
+            qml.RZ((phi - omega) / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.RZ(-(phi + omega) / 2, wires=wires[1]),
+            qml.RY(-theta / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.RY(theta / 2, wires=wires[1]),
+            qml.RZ(omega, wires=wires[1]),
+        ]
+
+    def decomposition(self):
+        r"""Representation of the operator as a product of other operators.
+
+        .. math:: O = O_1 O_2 \dots O_n
+
+        A ``DecompositionUndefinedError`` is raised if no representation by decomposition is defined.
+
+        .. seealso:: :meth:`~.Operator.compute_decomposition`.
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
+
+
+class ControlledPhaseShift(ControlledOp):
+    r"""A qubit controlled phase shift.
+
+    .. math:: CR_\phi(\phi) = \begin{bmatrix}
+                1 & 0 & 0 & 0 \\
+                0 & 1 & 0 & 0 \\
+                0 & 0 & 1 & 0 \\
+                0 & 0 & 0 & e^{i\phi}
+            \end{bmatrix}.
+
+    .. note:: The first wire provided corresponds to the **control qubit**.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe: :math:`\frac{d}{d\phi}f(CR_\phi(\phi)) = \frac{1}{2}\left[f(CR_\phi(\phi+\pi/2)) - f(CR_\phi(\phi-\pi/2))\right]`
+        where :math:`f` is an expectation value depending on :math:`CR_{\phi}(\phi)`.
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+
+    """
+    num_wires = 2
+    """int: Number of wires the operator acts on."""
+
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    name = "ControlledPhaseShift"
+    parameter_frequencies = [(1,)]
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(qml.PhaseShift(phi, wires=wires[1:]), control_wires=wires[0], id=id)
+
+    @classmethod
+    def _unflatten(cls, data, metadata):
+        base = data[0]
+        control_wires = metadata[0]
+        return cls(*base.data, wires=control_wires + base.wires)
+
+    @staticmethod
+    def compute_matrix(phi):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.ControlledPhaseShift.matrix`
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.ControlledPhaseShift.compute_matrix(torch.tensor(0.5))
+            tensor([[1.0+0.0j, 0.0+0.0j, 0.0+0.0j, 0.0000+0.0000j],
+                    [0.0+0.0j, 1.0+0.0j, 0.0+0.0j, 0.0000+0.0000j],
+                    [0.0+0.0j, 0.0+0.0j, 1.0+0.0j, 0.0000+0.0000j],
+                    [0.0+0.0j, 0.0+0.0j, 0.0+0.0j, 0.8776+0.4794j]])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            p = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
+            if qml.math.ndim(p) == 0:
+                return qml.math.diag([1, 1, 1, p])
+
+            ones = qml.math.ones_like(p)
+            diags = stack_last([ones, ones, ones, p])
+            return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+
+        signs = qml.math.array([0, 0, 0, 1], like=phi)
+        arg = 1j * phi
+
+        if qml.math.ndim(arg) == 0:
+            return qml.math.diag(qml.math.exp(arg * signs))
+
+        diags = qml.math.exp(qml.math.outer(arg, signs))
+        return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
+
+    @staticmethod
+    def compute_eigvals(phi, **_):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.ControlledPhaseShift.eigvals`
+
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.ControlledPhaseShift.compute_eigvals(torch.tensor(0.5))
+        tensor([1.0000+0.0000j, 1.0000+0.0000j, 1.0000+0.0000j, 0.8776+0.4794j])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phase = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
+            ones = qml.math.ones_like(phase)
+            return stack_last([ones, ones, ones, phase])
+
+        prefactors = qml.math.array([0, 0, 0, 1j], like=phi)
+        if qml.math.ndim(phi) == 0:
+            product = phi * prefactors
+        else:
+            product = qml.math.outer(phi, prefactors)
+        return qml.math.exp(product)
+
+    def eigvals(self):
+        return self.compute_eigvals(*self.parameters)
+
+    @staticmethod
+    def compute_decomposition(phi, wires):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+        .. seealso:: :meth:`~.ControlledPhaseShift.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.ControlledPhaseShift.compute_decomposition(1.234, wires=(0,1))
+        [PhaseShift(0.617, wires=[0]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(-0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(0.617, wires=[1])]
+
+        """
+        return [
+            qml.PhaseShift(phi / 2, wires=wires[0]),
+            qml.CNOT(wires=wires),
+            qml.PhaseShift(-phi / 2, wires=wires[1]),
+            qml.CNOT(wires=wires),
+            qml.PhaseShift(phi / 2, wires=wires[1]),
+        ]
+
+    def decomposition(self):
+        r"""Representation of the operator as a product of other operators.
+
+        .. math:: O = O_1 O_2 \dots O_n
+
+        A ``DecompositionUndefinedError`` is raised if no representation by decomposition is defined.
+
+        .. seealso:: :meth:`~.Operator.compute_decomposition`.
+
+        Returns:
+            list[Operator]: decomposition of the operator
+        """
+        return self.compute_decomposition(*self.parameters, wires=self.wires)
+
+
+CPhase = ControlledPhaseShift
diff --git a/pennylane/ops/qubit/__init__.py b/pennylane/ops/qubit/__init__.py
index c539ba378..d9484a05c 100644
--- a/pennylane/ops/qubit/__init__.py
+++ b/pennylane/ops/qubit/__init__.py
@@ -24,7 +24,6 @@ The operations are divided into the following files:
   located in ``non_parameteric_ops.py`` instead.
 * ``parametric_ops_single_qubit.py``: Core single qubit parametric operations.
 * ``parametric_ops_multi_qubit.py``: Core multi-qubit parametric operations.
-* ``parametric_ops_controlled.py``: Controlled parametric operations.
 * ``qchem_ops.py``: Operations for quantum chemistry applications.
 * ``state_preparation.py``: Operations that initialize the state.
 * ``special_unitary.py``: The ``SpecialUnitary`` operation.
@@ -36,7 +35,6 @@ from .non_parametric_ops import *
 from .observables import *
 from .parametric_ops_single_qubit import *
 from .parametric_ops_multi_qubit import *
-from .parametric_ops_controlled import *
 from .qchem_ops import *
 from .state_preparation import *
 from .special_unitary import SpecialUnitary
@@ -72,16 +70,10 @@ __ops__ = {
     "RZ",
     "PhaseShift",
     "PCPhase",
-    "ControlledPhaseShift",
     "CPhaseShift00",
     "CPhaseShift01",
     "CPhaseShift10",
-    "CPhase",
     "Rot",
-    "CRX",
-    "CRY",
-    "CRZ",
-    "CRot",
     "U1",
     "U2",
     "U3",
diff --git a/pennylane/ops/qubit/parametric_ops_controlled.py b/pennylane/ops/qubit/parametric_ops_controlled.py
deleted file mode 100644
index a037c83cc..000000000
--- a/pennylane/ops/qubit/parametric_ops_controlled.py
+++ /dev/null
@@ -1,1463 +0,0 @@
-# Copyright 2018-2023 Xanadu Quantum Technologies Inc.
-
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-
-#     http://www.apache.org/licenses/LICENSE-2.0
-
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-# pylint: disable=too-many-arguments
-"""
-This submodule contains the discrete-variable quantum operations that are the
-core parameterized gates.
-"""
-# pylint:disable=abstract-method,arguments-differ,protected-access,invalid-overridden-method
-import numpy as np
-
-import pennylane as qml
-from pennylane.operation import Operation
-from pennylane.wires import Wires
-
-from .non_parametric_ops import Hadamard, PauliX, PauliY, PauliZ
-from .parametric_ops_single_qubit import _can_replace, stack_last, RY, RZ, PhaseShift
-
-
-class ControlledPhaseShift(Operation):
-    r"""
-    A qubit controlled phase shift.
-
-    .. math:: CR_\phi(\phi) = \begin{bmatrix}
-                1 & 0 & 0 & 0 \\
-                0 & 1 & 0 & 0 \\
-                0 & 0 & 1 & 0 \\
-                0 & 0 & 0 & e^{i\phi}
-            \end{bmatrix}.
-
-    .. note:: The first wire provided corresponds to the **control qubit**.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe: :math:`\frac{d}{d\phi}f(CR_\phi(\phi)) = \frac{1}{2}\left[f(CR_\phi(\phi+\pi/2)) - f(CR_\phi(\phi-\pi/2))\right]`
-      where :math:`f` is an expectation value depending on :math:`CR_{\phi}(\phi)`.
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    basis = "Z"
-    grad_method = "A"
-    parameter_frequencies = [(1,)]
-
-    def generator(self):
-        return qml.Projector(np.array([1, 1]), wires=self.wires)
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label=base_label or "Rϕ", cache=cache)
-
-    @staticmethod
-    def compute_matrix(phi):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.ControlledPhaseShift.matrix`
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.ControlledPhaseShift.compute_matrix(torch.tensor(0.5))
-            tensor([[1.0+0.0j, 0.0+0.0j, 0.0+0.0j, 0.0000+0.0000j],
-                    [0.0+0.0j, 1.0+0.0j, 0.0+0.0j, 0.0000+0.0000j],
-                    [0.0+0.0j, 0.0+0.0j, 1.0+0.0j, 0.0000+0.0000j],
-                    [0.0+0.0j, 0.0+0.0j, 0.0+0.0j, 0.8776+0.4794j]])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            p = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
-            if qml.math.ndim(p) == 0:
-                return qml.math.diag([1, 1, 1, p])
-
-            ones = qml.math.ones_like(p)
-            diags = stack_last([ones, ones, ones, p])
-            return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
-
-        signs = qml.math.array([0, 0, 0, 1], like=phi)
-        arg = 1j * phi
-
-        if qml.math.ndim(arg) == 0:
-            return qml.math.diag(qml.math.exp(arg * signs))
-
-        diags = qml.math.exp(qml.math.outer(arg, signs))
-        return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
-
-    @staticmethod
-    def compute_eigvals(phi):  # pylint: disable=arguments-differ
-        r"""Eigenvalues of the operator in the computational basis (static method).
-
-        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
-        the operator can be reconstructed as
-
-        .. math:: O = U \Sigma U^{\dagger},
-
-        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
-
-        Otherwise, no particular order for the eigenvalues is guaranteed.
-
-        .. seealso:: :meth:`~.ControlledPhaseShift.eigvals`
-
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: eigenvalues
-
-        **Example**
-
-        >>> qml.ControlledPhaseShift.compute_eigvals(torch.tensor(0.5))
-        tensor([1.0000+0.0000j, 1.0000+0.0000j, 1.0000+0.0000j, 0.8776+0.4794j])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phase = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
-            ones = qml.math.ones_like(phase)
-            return stack_last([ones, ones, ones, phase])
-
-        prefactors = qml.math.array([0, 0, 0, 1j], like=phi)
-        if qml.math.ndim(phi) == 0:
-            product = phi * prefactors
-        else:
-            product = qml.math.outer(phi, prefactors)
-        return qml.math.exp(product)
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-
-        .. seealso:: :meth:`~.ControlledPhaseShift.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.ControlledPhaseShift.compute_decomposition(1.234, wires=(0,1))
-        [PhaseShift(0.617, wires=[0]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(-0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(0.617, wires=[1])]
-
-        """
-        decomp_ops = [
-            PhaseShift(phi / 2, wires=wires[0]),
-            qml.CNOT(wires=wires),
-            PhaseShift(-phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            PhaseShift(phi / 2, wires=wires[1]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return ControlledPhaseShift(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [ControlledPhaseShift(self.data[0] * z, wires=self.wires)]
-
-    def simplify(self):
-        phi = self.data[0] % (2 * np.pi)
-
-        if _can_replace(phi, 0):
-            return qml.Identity(wires=self.wires[0])
-
-        return ControlledPhaseShift(phi, wires=self.wires)
-
-    @property
-    def control_wires(self):
-        return Wires(self.wires[0])
-
-
-CPhase = ControlledPhaseShift
-
-
-class CPhaseShift00(Operation):
-    r"""
-    A qubit controlled phase shift.
-
-    .. math:: CR_{00}(\phi) = \begin{bmatrix}
-                e^{i\phi} & 0 & 0 & 0 \\
-                0 & 1 & 0 & 0 \\
-                0 & 0 & 1 & 0 \\
-                0 & 0 & 0 & 1
-            \end{bmatrix}.
-
-    .. note:: The first wire provided corresponds to the **control qubit** and controls
-        on the zero state :math:`|0\rangle`.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe:
-
-    .. math::
-        \frac{d}{d \phi} CR_{00}(\phi)
-        = \frac{1}{2} \left[ CR_{00}(\phi + \pi / 2)
-            - CR_{00}(\phi - \pi / 2) \right]
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    grad_method = "A"
-    parameter_frequencies = [(1,)]
-
-    def generator(self):
-        return qml.Projector(np.array([0, 0]), wires=self.wires)
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label="Rϕ(00)", cache=cache)
-
-    @staticmethod
-    def compute_matrix(phi):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CPhaseShift00.matrix`
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CPhaseShift00.compute_matrix(torch.tensor(0.5))
-            tensor([[0.8776+0.4794j, 0.0+0.0j, 0.0+0.0j, 0.0+0.0j],
-                    [0.0000+0.0000j, 1.0+0.0j, 0.0+0.0j, 0.0+0.0j],
-                    [0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j, 0.0+0.0j],
-                    [0.0000+0.0000j, 0.0+0.0j, 0.0+0.0j, 1.0+0.0j]])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
-            matrix = [
-                [exp_part, zeros, zeros, zeros],
-                [zeros, ones, zeros, zeros],
-                [zeros, zeros, ones, zeros],
-                [zeros, zeros, zeros, ones],
-            ]
-
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
-
-        return qml.math.diag([exp_part, 1, 1, 1])
-
-    @staticmethod
-    def compute_eigvals(phi):  # pylint: disable=arguments-differ
-        r"""Eigenvalues of the operator in the computational basis (static method).
-
-        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
-        the operator can be reconstructed as
-
-        .. math:: O = U \Sigma U^{\dagger},
-
-        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
-
-        Otherwise, no particular order for the eigenvalues is guaranteed.
-
-        .. seealso:: :meth:`~.CPhaseShift00.eigvals`
-
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: eigenvalues
-
-        **Example**
-
-        >>> qml.CPhaseShift00.compute_eigvals(torch.tensor(0.5))
-        tensor([0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j, 1.0000+0.0000j])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
-        return stack_last([exp_part, ones, ones, ones])
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-
-        .. seealso:: :meth:`~.CPhaseShift00.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CPhaseShift00.compute_decomposition(1.234, wires=(0,1))
-        [PauliX(wires=[0]),
-        PauliX(wires=[1]),
-        PhaseShift(0.617, wires=[0]),
-        PhaseShift(0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(-0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PauliX(wires=[1]),
-        PauliX(wires=[0])]
-
-        """
-        decomp_ops = [
-            PauliX(wires[0]),
-            PauliX(wires[1]),
-            PhaseShift(phi / 2, wires=[wires[0]]),
-            PhaseShift(phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PhaseShift(-phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PauliX(wires[1]),
-            PauliX(wires[0]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CPhaseShift00(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CPhaseShift00(self.data[0] * z, wires=self.wires)]
-
-    @property
-    def control_values(self):
-        """str: The control values of the operation"""
-        return "0"
-
-    @property
-    def control_wires(self):
-        return self.wires[0:1]
-
-
-class CPhaseShift01(Operation):
-    r"""
-    A qubit controlled phase shift.
-
-    .. math:: CR_{01\phi}(\phi) = \begin{bmatrix}
-                1 & 0 & 0 & 0 \\
-                0 & e^{i\phi} & 0 & 0 \\
-                0 & 0 & 1 & 0 \\
-                0 & 0 & 0 & 1
-            \end{bmatrix}.
-
-    .. note:: The first wire provided corresponds to the **control qubit** and controls
-        on the zero state :math:`|0\rangle`.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe:
-
-    .. math::
-        \frac{d}{d \phi} CR_{01}(\phi)
-        = \frac{1}{2} \left[ CR_{01}(\phi + \pi / 2)
-            - CR_{01}(\phi - \pi / 2) \right]
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    grad_method = "A"
-    parameter_frequencies = [(1,)]
-
-    def generator(self):
-        return qml.Projector(np.array([0, 1]), wires=self.wires)
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label="Rϕ(01)", cache=cache)
-
-    @staticmethod
-    def compute_matrix(phi):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CPhaseShift01.matrix`
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CPhaseShift01.compute_matrix(torch.tensor(0.5))
-            tensor([[1.0+0.0j, 0.0000+0.0000j, 0.0+0.0j, 0.0+0.0j],
-                    [0.0+0.0j, 0.8776+0.4794j, 0.0+0.0j, 0.0+0.0j],
-                    [0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j, 0.0+0.0j],
-                    [0.0+0.0j, 0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j]])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
-            matrix = [
-                [ones, zeros, zeros, zeros],
-                [zeros, exp_part, zeros, zeros],
-                [zeros, zeros, ones, zeros],
-                [zeros, zeros, zeros, ones],
-            ]
-
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
-
-        return qml.math.diag([1, exp_part, 1, 1])
-
-    @staticmethod
-    def compute_eigvals(phi):  # pylint: disable=arguments-differ
-        r"""Eigenvalues of the operator in the computational basis (static method).
-
-        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
-        the operator can be reconstructed as
-
-        .. math:: O = U \Sigma U^{\dagger},
-
-        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
-
-        Otherwise, no particular order for the eigenvalues is guaranteed.
-
-        .. seealso:: :meth:`~.CPhaseShift01.eigvals`
-
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: eigenvalues
-
-        **Example**
-
-        >>> qml.CPhaseShift01.compute_eigvals(torch.tensor(0.5))
-        tensor([1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
-        return stack_last([ones, exp_part, ones, ones])
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-        .. seealso:: :meth:`~.CPhaseShift01.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CPhaseShift01.compute_decomposition(1.234, wires=(0,1))
-        [PauliX(wires=[0]),
-        PhaseShift(0.617, wires=[0]),
-        PhaseShift(0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(-0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PauliX(wires=[0])]
-
-        """
-        decomp_ops = [
-            PauliX(wires[0]),
-            PhaseShift(phi / 2, wires=[wires[0]]),
-            PhaseShift(phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PhaseShift(-phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PauliX(wires[0]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CPhaseShift01(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CPhaseShift01(self.data[0] * z, wires=self.wires)]
-
-    @property
-    def control_values(self):
-        """str: The control values of the operation"""
-        return "0"
-
-    @property
-    def control_wires(self):
-        return self.wires[0:1]
-
-
-class CPhaseShift10(Operation):
-    r"""
-    A qubit controlled phase shift.
-
-    .. math:: CR_{10\phi}(\phi) = \begin{bmatrix}
-                1 & 0 & 0 & 0 \\
-                0 & 1 & 0 & 0 \\
-                0 & 0 & e^{i\phi} & 0 \\
-                0 & 0 & 0 & 1
-            \end{bmatrix}.
-
-    .. note:: The first wire provided corresponds to the **control qubit**.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe:
-
-    .. math::
-        \frac{d}{d \phi} CR_{10}(\phi)
-        = \frac{1}{2} \left[ CR_{10}(\phi + \pi / 2)
-            - CR_{10}(\phi - \pi / 2) \right]
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Any, Wires): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    grad_method = "A"
-    parameter_frequencies = [(1,)]
-
-    def generator(self):
-        return qml.Projector(np.array([1, 0]), wires=self.wires)
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label="Rϕ(10)", cache=cache)
-
-    @staticmethod
-    def compute_matrix(phi):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CPhaseShift10.matrix`
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CPhaseShift10.compute_matrix(torch.tensor(0.5))
-            tensor([[1.0+0.0j, 0.0+0.0j, 0.0000+0.0000j, 0.0+0.0j],
-                    [0.0+0.0j, 1.0+0.0j, 0.0000+0.0000j, 0.0+0.0j],
-                    [0.0+0.0j, 0.0+0.0j, 0.8776+0.4794j, 0.0+0.0j],
-                    [0.0+0.0j, 0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j]])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-
-        if qml.math.ndim(phi) > 0:
-            ones = qml.math.ones_like(exp_part)
-            zeros = qml.math.zeros_like(exp_part)
-            matrix = [
-                [ones, zeros, zeros, zeros],
-                [zeros, ones, zeros, zeros],
-                [zeros, zeros, exp_part, zeros],
-                [zeros, zeros, zeros, ones],
-            ]
-
-            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
-
-        return qml.math.diag([1, 1, exp_part, 1])
-
-    @staticmethod
-    def compute_eigvals(phi):  # pylint: disable=arguments-differ
-        r"""Eigenvalues of the operator in the computational basis (static method).
-
-        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
-        the operator can be reconstructed as
-
-        .. math:: O = U \Sigma U^{\dagger},
-
-        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
-
-        Otherwise, no particular order for the eigenvalues is guaranteed.
-
-        .. seealso:: :meth:`~.CPhaseShift10.eigvals`
-
-
-        Args:
-            phi (tensor_like or float): phase shift
-
-        Returns:
-            tensor_like: eigenvalues
-
-        **Example**
-
-        >>> qml.CPhaseShift10.compute_eigvals(torch.tensor(0.5))
-        tensor([1.0000+0.0000j, 1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j])
-        """
-        if qml.math.get_interface(phi) == "tensorflow":
-            phi = qml.math.cast_like(phi, 1j)
-
-        exp_part = qml.math.exp(1j * phi)
-        ones = qml.math.ones_like(exp_part)
-        return stack_last([ones, ones, exp_part, ones])
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-        .. seealso:: :meth:`~.CPhaseShift10.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CPhaseShift10.compute_decomposition(1.234, wires=(0,1))
-        [PauliX(wires=[1]),
-        PhaseShift(0.617, wires=[0]),
-        PhaseShift(0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(-0.617, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PauliX(wires=[1])]
-
-        """
-        decomp_ops = [
-            PauliX(wires[1]),
-            PhaseShift(phi / 2, wires=[wires[0]]),
-            PhaseShift(phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PhaseShift(-phi / 2, wires=[wires[1]]),
-            qml.CNOT(wires=wires),
-            PauliX(wires[1]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CPhaseShift10(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CPhaseShift10(self.data[0] * z, wires=self.wires)]
-
-    @property
-    def control_wires(self):
-        return self.wires[0:1]
-
-
-class CRX(Operation):
-    r"""
-    The controlled-RX operator
-
-    .. math::
-
-        \begin{align}
-            CR_x(\phi) &=
-            \begin{bmatrix}
-            & 1 & 0 & 0 & 0 \\
-            & 0 & 1 & 0 & 0\\
-            & 0 & 0 & \cos(\phi/2) & -i\sin(\phi/2)\\
-            & 0 & 0 & -i\sin(\phi/2) & \cos(\phi/2)
-            \end{bmatrix}.
-        \end{align}
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe: The controlled-RX operator satisfies a four-term parameter-shift rule
-      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
-
-      .. math::
-
-          \frac{d}{d\phi}f(CR_x(\phi)) = c_+ \left[f(CR_x(\phi+a)) - f(CR_x(\phi-a))\right] - c_- \left[f(CR_x(\phi+b)) - f(CR_x(\phi-b))\right]
-
-      where :math:`f` is an expectation value depending on :math:`CR_x(\phi)`, and
-
-      - :math:`a = \pi/2`
-      - :math:`b = 3\pi/2`
-      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    basis = "X"
-    grad_method = "A"
-    parameter_frequencies = [(0.5, 1.0)]
-
-    def generator(self):
-        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ PauliX(self.wires[1])
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label=base_label or "RX", cache=cache)
-
-    @staticmethod
-    def compute_matrix(theta):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CRX.matrix`
-
-        Args:
-            theta (tensor_like or float): rotation angle
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CRX.compute_matrix(torch.tensor(0.5))
-        tensor([[1.0+0.0j, 0.0+0.0j,    0.0+0.0j,    0.0+0.0j],
-                [0.0+0.0j, 1.0+0.0j,    0.0+0.0j,    0.0+0.0j],
-                [0.0+0.0j, 0.0+0.0j, 0.9689+0.0j, 0.0-0.2474j],
-                [0.0+0.0j, 0.0+0.0j, 0.0-0.2474j, 0.9689+0.0j]])
-        """
-        interface = qml.math.get_interface(theta)
-
-        c = qml.math.cos(theta / 2)
-        s = qml.math.sin(theta / 2)
-
-        if interface == "tensorflow":
-            c = qml.math.cast_like(c, 1j)
-            s = qml.math.cast_like(s, 1j)
-
-        # The following avoids casting an imaginary quantity to reals when backpropagating
-        c = (1 + 0j) * c
-        js = -1j * s
-        ones = qml.math.ones_like(js)
-        zeros = qml.math.zeros_like(js)
-        matrix = [
-            [ones, zeros, zeros, zeros],
-            [zeros, ones, zeros, zeros],
-            [zeros, zeros, c, js],
-            [zeros, zeros, js, c],
-        ]
-
-        return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-        .. seealso:: :meth:`~.CRot.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): the wires the operation acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CRX.compute_decomposition(1.2, wires=(0,1))
-        [RZ(1.5707963267948966, wires=[1]),
-        RY(0.6, wires=[1]),
-        CNOT(wires=[0, 1]),
-        RY(-0.6, wires=[1]),
-        CNOT(wires=[0, 1]),
-        RZ(-1.5707963267948966, wires=[1])]
-
-        """
-        pi_half = qml.math.ones_like(phi) * (np.pi / 2)
-        decomp_ops = [
-            RZ(pi_half, wires=wires[1]),
-            RY(phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            RY(-phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            RZ(-pi_half, wires=wires[1]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CRX(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CRX(self.data[0] * z, wires=self.wires)]
-
-    def simplify(self):
-        phi = self.data[0] % (4 * np.pi)
-
-        if _can_replace(phi, 0):
-            return qml.Identity(wires=self.wires[0])
-
-        return CRX(phi, wires=self.wires)
-
-    @property
-    def control_wires(self):
-        return Wires(self.wires[0])
-
-
-class CRY(Operation):
-    r"""
-    The controlled-RY operator
-
-    .. math::
-
-        \begin{align}
-            CR_y(\phi) &=
-            \begin{bmatrix}
-                1 & 0 & 0 & 0 \\
-                0 & 1 & 0 & 0\\
-                0 & 0 & \cos(\phi/2) & -\sin(\phi/2)\\
-                0 & 0 & \sin(\phi/2) & \cos(\phi/2)
-            \end{bmatrix}.
-        \end{align}
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe: The controlled-RY operator satisfies a four-term parameter-shift rule
-      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
-
-      .. math::
-
-          \frac{d}{d\phi}f(CR_y(\phi)) = c_+ \left[f(CR_y(\phi+a)) - f(CR_y(\phi-a))\right] - c_- \left[f(CR_y(\phi+b)) - f(CR_y(\phi-b))\right]
-
-      where :math:`f` is an expectation value depending on :math:`CR_y(\phi)`, and
-
-      - :math:`a = \pi/2`
-      - :math:`b = 3\pi/2`
-      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    basis = "Y"
-    grad_method = "A"
-    parameter_frequencies = [(0.5, 1.0)]
-
-    def generator(self):
-        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ PauliY(self.wires[1])
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label=base_label or "RY", cache=cache)
-
-    @staticmethod
-    def compute_matrix(theta):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CRY.matrix`
-
-
-        Args:
-            theta (tensor_like or float): rotation angle
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CRY.compute_matrix(torch.tensor(0.5))
-        tensor([[ 1.0000,  0.0000,  0.0000,  0.0000],
-                [ 0.0000,  1.0000,  0.0000,  0.0000],
-                [ 0.0000,  0.0000,  0.9689, -0.2474],
-                [ 0.0000,  0.0000,  0.2474,  0.9689]], dtype=torch.float64)
-        """
-        interface = qml.math.get_interface(theta)
-
-        c = qml.math.cos(theta / 2)
-        s = qml.math.sin(theta / 2)
-
-        if interface == "tensorflow":
-            c = qml.math.cast_like(c, 1j)
-            s = qml.math.cast_like(s, 1j)
-
-        # The following avoids casting an imaginary quantity to reals when backpropagating
-        c = (1 + 0j) * c
-        s = (1 + 0j) * s
-        ones = qml.math.ones_like(s)
-        zeros = qml.math.zeros_like(s)
-        matrix = [
-            [ones, zeros, zeros, zeros],
-            [zeros, ones, zeros, zeros],
-            [zeros, zeros, c, -s],
-            [zeros, zeros, s, c],
-        ]
-
-        return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-        .. seealso:: :meth:`~.CRY.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CRY.compute_decomposition(1.2, wires=(0,1))
-        [RY(0.6, wires=[1]),
-        CNOT(wires=[0, 1]),
-        RY(-0.6, wires=[1]),
-        CNOT(wires=[0, 1])]
-
-        """
-        decomp_ops = [
-            RY(phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            RY(-phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CRY(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CRY(self.data[0] * z, wires=self.wires)]
-
-    def simplify(self):
-        phi = self.data[0] % (4 * np.pi)
-
-        if _can_replace(phi, 0):
-            return qml.Identity(wires=self.wires[0])
-
-        return CRY(phi, wires=self.wires)
-
-    @property
-    def control_wires(self):
-        return Wires(self.wires[0])
-
-
-class CRZ(Operation):
-    r"""
-    The controlled-RZ operator
-
-    .. math::
-
-        \begin{align}
-             CR_z(\phi) &=
-             \begin{bmatrix}
-                1 & 0 & 0 & 0 \\
-                0 & 1 & 0 & 0\\
-                0 & 0 & e^{-i\phi/2} & 0\\
-                0 & 0 & 0 & e^{i\phi/2}
-            \end{bmatrix}.
-        \end{align}
-
-
-    .. note:: The subscripts of the operations in the formula refer to the wires they act on, e.g. 1 corresponds to the first element in ``wires`` that is the **control qubit**.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 1
-    * Number of dimensions per parameter: (0,)
-    * Gradient recipe: The controlled-RZ operator satisfies a four-term parameter-shift rule
-      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
-
-      .. math::
-
-          \frac{d}{d\phi}f(CR_z(\phi)) = c_+ \left[f(CR_z(\phi+a)) - f(CR_z(\phi-a))\right] - c_- \left[f(CR_z(\phi+b)) - f(CR_z(\phi-b))\right]
-
-      where :math:`f` is an expectation value depending on :math:`CR_z(\phi)`, and
-
-      - :math:`a = \pi/2`
-      - :math:`b = 3\pi/2`
-      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 1
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0,)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    basis = "Z"
-    grad_method = "A"
-    parameter_frequencies = [(0.5, 1.0)]
-
-    def generator(self):
-        return -0.5 * qml.Projector(np.array([1]), wires=self.wires[0]) @ PauliZ(self.wires[1])
-
-    def __init__(self, phi, wires, id=None):
-        super().__init__(phi, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label=base_label or "RZ", cache=cache)
-
-    @staticmethod
-    def compute_matrix(theta):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CRZ.matrix`
-
-        Args:
-            theta (tensor_like or float): rotation angle
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-        >>> qml.CRZ.compute_matrix(torch.tensor(0.5))
-        tensor([[1.0+0.0j, 0.0+0.0j,       0.0+0.0j,       0.0+0.0j],
-                [0.0+0.0j, 1.0+0.0j,       0.0+0.0j,       0.0+0.0j],
-                [0.0+0.0j, 0.0+0.0j, 0.9689-0.2474j,       0.0+0.0j],
-                [0.0+0.0j, 0.0+0.0j,       0.0+0.0j, 0.9689+0.2474j]])
-        """
-        if qml.math.get_interface(theta) == "tensorflow":
-            p = qml.math.exp(-0.5j * qml.math.cast_like(theta, 1j))
-            if qml.math.ndim(p) == 0:
-                return qml.math.diag([1, 1, p, qml.math.conj(p)])
-
-            ones = qml.math.ones_like(p)
-            diags = stack_last([ones, ones, p, qml.math.conj(p)])
-            return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
-
-        signs = qml.math.array([0, 0, 1, -1], like=theta)
-        arg = -0.5j * theta
-
-        if qml.math.ndim(arg) == 0:
-            return qml.math.diag(qml.math.exp(arg * signs))
-
-        diags = qml.math.exp(qml.math.outer(arg, signs))
-        return diags[:, :, np.newaxis] * qml.math.cast_like(qml.math.eye(4, like=diags), diags)
-
-    @staticmethod
-    def compute_eigvals(theta):  # pylint: disable=arguments-differ
-        r"""Eigenvalues of the operator in the computational basis (static method).
-
-        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
-        the operator can be reconstructed as
-
-        .. math:: O = U \Sigma U^{\dagger},
-
-        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
-
-        Otherwise, no particular order for the eigenvalues is guaranteed.
-
-        .. seealso:: :meth:`~.CRZ.eigvals`
-
-
-        Args:
-            theta (tensor_like or float): rotation angle
-
-        Returns:
-            tensor_like: eigenvalues
-
-        **Example**
-
-        >>> qml.CRZ.compute_eigvals(torch.tensor(0.5))
-        tensor([1.0000+0.0000j, 1.0000+0.0000j, 0.9689-0.2474j, 0.9689+0.2474j])
-        """
-        if qml.math.get_interface(theta) == "tensorflow":
-            phase = qml.math.exp(-0.5j * qml.math.cast_like(theta, 1j))
-            ones = qml.math.ones_like(phase)
-            return stack_last([ones, ones, phase, qml.math.conj(phase)])
-
-        prefactors = qml.math.array([0, 0, -0.5j, 0.5j], like=theta)
-        if qml.math.ndim(theta) == 0:
-            product = theta * prefactors
-        else:
-            product = qml.math.outer(theta, prefactors)
-        return qml.math.exp(product)
-
-    @staticmethod
-    def compute_decomposition(phi, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-        .. seealso:: :meth:`~.CRZ.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            wires (Iterable, Wires): wires that the operator acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.CRZ.compute_decomposition(1.2, wires=(0,1))
-        [PhaseShift(0.6, wires=[1]),
-        CNOT(wires=[0, 1]),
-        PhaseShift(-0.6, wires=[1]),
-        CNOT(wires=[0, 1])]
-
-        """
-        decomp_ops = [
-            PhaseShift(phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            PhaseShift(-phi / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        return CRZ(-self.data[0], wires=self.wires)
-
-    def pow(self, z):
-        return [CRZ(self.data[0] * z, wires=self.wires)]
-
-    def simplify(self):
-        phi = self.data[0] % (4 * np.pi)
-
-        if _can_replace(phi, 0):
-            return qml.Identity(wires=self.wires[0])
-
-        return CRZ(phi, wires=self.wires)
-
-    @property
-    def control_wires(self):
-        return Wires(self.wires[0])
-
-
-class CRot(Operation):
-    r"""
-    The controlled-Rot operator
-
-    .. math:: CR(\phi, \theta, \omega) = \begin{bmatrix}
-            1 & 0 & 0 & 0 \\
-            0 & 1 & 0 & 0\\
-            0 & 0 & e^{-i(\phi+\omega)/2}\cos(\theta/2) & -e^{i(\phi-\omega)/2}\sin(\theta/2)\\
-            0 & 0 & e^{-i(\phi-\omega)/2}\sin(\theta/2) & e^{i(\phi+\omega)/2}\cos(\theta/2)
-        \end{bmatrix}.
-
-    .. note:: The first wire provided corresponds to the **control qubit**.
-
-    **Details:**
-
-    * Number of wires: 2
-    * Number of parameters: 3
-    * Number of dimensions per parameter: (0, 0, 0)
-    * Gradient recipe: The controlled-Rot operator satisfies a four-term parameter-shift rule
-      (see Appendix F, https://doi.org/10.1088/1367-2630/ac2cb3):
-
-      .. math::
-
-          \frac{d}{d\mathbf{x}_i}f(CR(\mathbf{x}_i)) = c_+ \left[f(CR(\mathbf{x}_i+a)) - f(CR(\mathbf{x}_i-a))\right] - c_- \left[f(CR(\mathbf{x}_i+b)) - f(CR(\mathbf{x}_i-b))\right]
-
-      where :math:`f` is an expectation value depending on :math:`CR(\mathbf{x}_i)`, and
-
-      - :math:`\mathbf{x} = (\phi, \theta, \omega)` and `i` is an index to :math:`\mathbf{x}`
-      - :math:`a = \pi/2`
-      - :math:`b = 3\pi/2`
-      - :math:`c_{\pm} = (\sqrt{2} \pm 1)/{4\sqrt{2}}`
-
-    Args:
-        phi (float): rotation angle :math:`\phi`
-        theta (float): rotation angle :math:`\theta`
-        omega (float): rotation angle :math:`\omega`
-        wires (Sequence[int]): the wire the operation acts on
-        id (str or None): String representing the operation (optional)
-    """
-    num_wires = 2
-    num_params = 3
-    """int: Number of trainable parameters that the operator depends on."""
-
-    ndim_params = (0, 0, 0)
-    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
-
-    grad_method = "A"
-    parameter_frequencies = [(0.5, 1.0), (0.5, 1.0), (0.5, 1.0)]
-
-    def __init__(self, phi, theta, omega, wires, id=None):
-        super().__init__(phi, theta, omega, wires=wires, id=id)
-
-    def label(self, decimals=None, base_label=None, cache=None):
-        return super().label(decimals=decimals, base_label=base_label or "Rot", cache=cache)
-
-    @staticmethod
-    def compute_matrix(phi, theta, omega):  # pylint: disable=arguments-differ
-        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
-
-        The canonical matrix is the textbook matrix representation that does not consider wires.
-        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
-
-        .. seealso:: :meth:`~.CRot.matrix`
-
-
-        Args:
-            phi(tensor_like or float): first rotation angle
-            theta (tensor_like or float): second rotation angle
-            omega (tensor_like or float): third rotation angle
-
-        Returns:
-            tensor_like: canonical matrix
-
-        **Example**
-
-         >>> qml.CRot.compute_matrix(torch.tensor(0.1), torch.tensor(0.2), torch.tensor(0.3))
-         tensor([[ 1.0+0.0j,  0.0+0.0j,        0.0+0.0j,        0.0+0.0j],
-                [ 0.0+0.0j,  1.0+0.0j,        0.0+0.0j,        0.0+0.0j],
-                [ 0.0+0.0j,  0.0+0.0j,  0.9752-0.1977j, -0.0993+0.0100j],
-                [ 0.0+0.0j,  0.0+0.0j,  0.0993+0.0100j,  0.9752+0.1977j]])
-        """
-        # It might be that they are in different interfaces, e.g.,
-        # CRot(0.2, 0.3, tf.Variable(0.5), wires=[0, 1])
-        # So we need to make sure the matrix comes out having the right type
-        interface = qml.math.get_interface(phi, theta, omega)
-
-        c = qml.math.cos(theta / 2)
-        s = qml.math.sin(theta / 2)
-
-        # If anything is not tensorflow, it has to be casted
-        if interface == "tensorflow":
-            phi = qml.math.cast_like(qml.math.asarray(phi, like=interface), 1j)
-            omega = qml.math.cast_like(qml.math.asarray(omega, like=interface), 1j)
-            c = qml.math.cast_like(qml.math.asarray(c, like=interface), 1j)
-            s = qml.math.cast_like(qml.math.asarray(s, like=interface), 1j)
-
-        # The following variable is used to assert the all terms to be stacked have same shape
-        one = qml.math.ones_like(phi) * qml.math.ones_like(omega)
-        c = c * one
-        s = s * one
-
-        o = qml.math.ones_like(c)
-        z = qml.math.zeros_like(c)
-        mat = [
-            [o, z, z, z],
-            [z, o, z, z],
-            [
-                z,
-                z,
-                qml.math.exp(-0.5j * (phi + omega)) * c,
-                -qml.math.exp(0.5j * (phi - omega)) * s,
-            ],
-            [
-                z,
-                z,
-                qml.math.exp(-0.5j * (phi - omega)) * s,
-                qml.math.exp(0.5j * (phi + omega)) * c,
-            ],
-        ]
-
-        return qml.math.stack([stack_last(row) for row in mat], axis=-2)
-
-    @staticmethod
-    def compute_decomposition(phi, theta, omega, wires):
-        r"""Representation of the operator as a product of other operators (static method). :
-
-        .. math:: O = O_1 O_2 \dots O_n.
-
-
-        .. seealso:: :meth:`~.CRot.decomposition`.
-
-        Args:
-            phi (float): rotation angle :math:`\phi`
-            theta (float): rotation angle :math:`\theta`
-            omega (float): rotation angle :math:`\omega`
-            wires (Iterable, Wires): the wires the operation acts on
-
-        Returns:
-            list[Operator]: decomposition into lower level operations
-
-        **Example:**
-
-        >>> qml.PhaseShift.compute_decomposition(1.234, wires=0)
-        [RZ(-1.1, wires=[1]),
-        CNOT(wires=[0, 1]),
-        RZ(-2.3, wires=[1]),
-        RY(-1.15, wires=[1]),
-        CNOT(wires=[0, 1]),
-        RY(1.15, wires=[1]),
-        RZ(3.4, wires=[1])]
-
-        """
-        decomp_ops = [
-            RZ((phi - omega) / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            RZ(-(phi + omega) / 2, wires=wires[1]),
-            RY(-theta / 2, wires=wires[1]),
-            qml.CNOT(wires=wires),
-            RY(theta / 2, wires=wires[1]),
-            RZ(omega, wires=wires[1]),
-        ]
-        return decomp_ops
-
-    def adjoint(self):
-        phi, theta, omega = self.parameters
-        return CRot(-omega, -theta, -phi, wires=self.wires)
-
-    @property
-    def control_wires(self):
-        return Wires(self.wires[0])
-
-    def simplify(self):
-        """Simplifies into single controlled rotations or a controlled-Hadamard if possible.
-
-        >>> qml.CRot(np.pi / 2, 0.1, -np.pi / 2, wires=(0,1)).simplify()
-        CRX(0.1, wires=[0, 1])
-        >>> qml.CRot(0, 0.2, 0, wires=(0,1)).simplify()
-        CRY(0.2, wires=[0, 1])
-
-        """
-        target_wires = [w for w in self.wires if w not in self.control_wires]
-        wires = self.wires
-        params = self.parameters
-
-        p0, p1, p2 = [p % (4 * np.pi) for p in params]
-
-        if _can_replace(p0, 0) and _can_replace(p1, 0) and _can_replace(p2, 0):
-            return qml.Identity(wires=wires[0])
-        if _can_replace(p0, np.pi / 2) and _can_replace(p2, 7 * np.pi / 2):
-            return qml.CRX(p1, wires=wires)
-        if _can_replace(p0, 0) and _can_replace(p2, 0):
-            return qml.CRY(p1, wires=wires)
-        if _can_replace(p1, 0):
-            return qml.CRZ((p0 + p2) % (4 * np.pi), wires=wires)
-        if _can_replace(p0, np.pi) and _can_replace(p1, np.pi / 2) and _can_replace(p2, 0):
-            return qml.ctrl(Hadamard(wires=target_wires), control=self.control_wires)
-
-        return CRot(p0, p1, p2, wires=wires)
diff --git a/pennylane/ops/qubit/parametric_ops_multi_qubit.py b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
index 21a2265ed..b9a8fddc7 100644
--- a/pennylane/ops/qubit/parametric_ops_multi_qubit.py
+++ b/pennylane/ops/qubit/parametric_ops_multi_qubit.py
@@ -1516,3 +1516,541 @@ class PSWAP(Operation):
             return qml.SWAP(wires=self.wires)
 
         return PSWAP(phi, wires=self.wires)
+
+
+class CPhaseShift00(Operation):
+    r"""
+    A qubit controlled phase shift.
+
+    .. math:: CR_{00}(\phi) = \begin{bmatrix}
+                e^{i\phi} & 0 & 0 & 0 \\
+                0 & 1 & 0 & 0 \\
+                0 & 0 & 1 & 0 \\
+                0 & 0 & 0 & 1
+            \end{bmatrix}.
+
+    .. note:: The first wire provided corresponds to the **control qubit** and controls
+        on the zero state :math:`|0\rangle`.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe:
+
+    .. math::
+        \frac{d}{d \phi} CR_{00}(\phi)
+        = \frac{1}{2} \left[ CR_{00}(\phi + \pi / 2)
+            - CR_{00}(\phi - \pi / 2) \right]
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    grad_method = "A"
+    parameter_frequencies = [(1,)]
+
+    def generator(self):
+        return qml.Projector(np.array([0, 0]), wires=self.wires)
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(phi, wires=wires, id=id)
+
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label="Rϕ(00)", cache=cache)
+
+    @staticmethod
+    def compute_matrix(phi):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CPhaseShift00.matrix`
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CPhaseShift00.compute_matrix(torch.tensor(0.5))
+            tensor([[0.8776+0.4794j, 0.0+0.0j, 0.0+0.0j, 0.0+0.0j],
+                    [0.0000+0.0000j, 1.0+0.0j, 0.0+0.0j, 0.0+0.0j],
+                    [0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j, 0.0+0.0j],
+                    [0.0000+0.0000j, 0.0+0.0j, 0.0+0.0j, 1.0+0.0j]])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+
+        if qml.math.ndim(phi) > 0:
+            ones = qml.math.ones_like(exp_part)
+            zeros = qml.math.zeros_like(exp_part)
+            matrix = [
+                [exp_part, zeros, zeros, zeros],
+                [zeros, ones, zeros, zeros],
+                [zeros, zeros, ones, zeros],
+                [zeros, zeros, zeros, ones],
+            ]
+
+            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+
+        return qml.math.diag([exp_part, 1, 1, 1])
+
+    @staticmethod
+    def compute_eigvals(phi):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.CPhaseShift00.eigvals`
+
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.CPhaseShift00.compute_eigvals(torch.tensor(0.5))
+        tensor([0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j, 1.0000+0.0000j])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+        ones = qml.math.ones_like(exp_part)
+        return stack_last([exp_part, ones, ones, ones])
+
+    @staticmethod
+    def compute_decomposition(phi, wires):
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+
+
+
+        .. seealso:: :meth:`~.CPhaseShift00.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CPhaseShift00.compute_decomposition(1.234, wires=(0,1))
+        [PauliX(wires=[0]),
+        PauliX(wires=[1]),
+        PhaseShift(0.617, wires=[0]),
+        PhaseShift(0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(-0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PauliX(wires=[1]),
+        PauliX(wires=[0])]
+
+        """
+        decomp_ops = [
+            PauliX(wires[0]),
+            PauliX(wires[1]),
+            PhaseShift(phi / 2, wires=[wires[0]]),
+            PhaseShift(phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PhaseShift(-phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PauliX(wires[1]),
+            PauliX(wires[0]),
+        ]
+        return decomp_ops
+
+    def adjoint(self):
+        return CPhaseShift00(-self.data[0], wires=self.wires)
+
+    def pow(self, z):
+        return [CPhaseShift00(self.data[0] * z, wires=self.wires)]
+
+    @property
+    def control_values(self):
+        """str: The control values of the operation"""
+        return "0"
+
+    @property
+    def control_wires(self):
+        return self.wires[0:1]
+
+
+class CPhaseShift01(Operation):
+    r"""
+    A qubit controlled phase shift.
+
+    .. math:: CR_{01\phi}(\phi) = \begin{bmatrix}
+                1 & 0 & 0 & 0 \\
+                0 & e^{i\phi} & 0 & 0 \\
+                0 & 0 & 1 & 0 \\
+                0 & 0 & 0 & 1
+            \end{bmatrix}.
+
+    .. note:: The first wire provided corresponds to the **control qubit** and controls
+        on the zero state :math:`|0\rangle`.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe:
+
+    .. math::
+        \frac{d}{d \phi} CR_{01}(\phi)
+        = \frac{1}{2} \left[ CR_{01}(\phi + \pi / 2)
+            - CR_{01}(\phi - \pi / 2) \right]
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Sequence[int]): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    grad_method = "A"
+    parameter_frequencies = [(1,)]
+
+    def generator(self):
+        return qml.Projector(np.array([0, 1]), wires=self.wires)
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(phi, wires=wires, id=id)
+
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label="Rϕ(01)", cache=cache)
+
+    @staticmethod
+    def compute_matrix(phi):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CPhaseShift01.matrix`
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CPhaseShift01.compute_matrix(torch.tensor(0.5))
+            tensor([[1.0+0.0j, 0.0000+0.0000j, 0.0+0.0j, 0.0+0.0j],
+                    [0.0+0.0j, 0.8776+0.4794j, 0.0+0.0j, 0.0+0.0j],
+                    [0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j, 0.0+0.0j],
+                    [0.0+0.0j, 0.0000+0.0000j, 0.0+0.0j, 1.0+0.0j]])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+
+        if qml.math.ndim(phi) > 0:
+            ones = qml.math.ones_like(exp_part)
+            zeros = qml.math.zeros_like(exp_part)
+            matrix = [
+                [ones, zeros, zeros, zeros],
+                [zeros, exp_part, zeros, zeros],
+                [zeros, zeros, ones, zeros],
+                [zeros, zeros, zeros, ones],
+            ]
+
+            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+
+        return qml.math.diag([1, exp_part, 1, 1])
+
+    @staticmethod
+    def compute_eigvals(phi):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.CPhaseShift01.eigvals`
+
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.CPhaseShift01.compute_eigvals(torch.tensor(0.5))
+        tensor([1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j, 1.0000+0.0000j])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+        ones = qml.math.ones_like(exp_part)
+        return stack_last([ones, exp_part, ones, ones])
+
+    @staticmethod
+    def compute_decomposition(phi, wires):
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+        .. seealso:: :meth:`~.CPhaseShift01.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CPhaseShift01.compute_decomposition(1.234, wires=(0,1))
+        [PauliX(wires=[0]),
+        PhaseShift(0.617, wires=[0]),
+        PhaseShift(0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(-0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PauliX(wires=[0])]
+
+        """
+        decomp_ops = [
+            PauliX(wires[0]),
+            PhaseShift(phi / 2, wires=[wires[0]]),
+            PhaseShift(phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PhaseShift(-phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PauliX(wires[0]),
+        ]
+        return decomp_ops
+
+    def adjoint(self):
+        return CPhaseShift01(-self.data[0], wires=self.wires)
+
+    def pow(self, z):
+        return [CPhaseShift01(self.data[0] * z, wires=self.wires)]
+
+    @property
+    def control_values(self):
+        """str: The control values of the operation"""
+        return "0"
+
+    @property
+    def control_wires(self):
+        return self.wires[0:1]
+
+
+class CPhaseShift10(Operation):
+    r"""
+    A qubit controlled phase shift.
+
+    .. math:: CR_{10\phi}(\phi) = \begin{bmatrix}
+                1 & 0 & 0 & 0 \\
+                0 & 1 & 0 & 0 \\
+                0 & 0 & e^{i\phi} & 0 \\
+                0 & 0 & 0 & 1
+            \end{bmatrix}.
+
+    .. note:: The first wire provided corresponds to the **control qubit**.
+
+    **Details:**
+
+    * Number of wires: 2
+    * Number of parameters: 1
+    * Number of dimensions per parameter: (0,)
+    * Gradient recipe:
+
+    .. math::
+        \frac{d}{d \phi} CR_{10}(\phi)
+        = \frac{1}{2} \left[ CR_{10}(\phi + \pi / 2)
+            - CR_{10}(\phi - \pi / 2) \right]
+
+    Args:
+        phi (float): rotation angle :math:`\phi`
+        wires (Any, Wires): the wire the operation acts on
+        id (str or None): String representing the operation (optional)
+    """
+    num_wires = 2
+    num_params = 1
+    """int: Number of trainable parameters that the operator depends on."""
+
+    ndim_params = (0,)
+    """tuple[int]: Number of dimensions per trainable parameter that the operator depends on."""
+
+    grad_method = "A"
+    parameter_frequencies = [(1,)]
+
+    def generator(self):
+        return qml.Projector(np.array([1, 0]), wires=self.wires)
+
+    def __init__(self, phi, wires, id=None):
+        super().__init__(phi, wires=wires, id=id)
+
+    def label(self, decimals=None, base_label=None, cache=None):
+        return super().label(decimals=decimals, base_label="Rϕ(10)", cache=cache)
+
+    @staticmethod
+    def compute_matrix(phi):  # pylint: disable=arguments-differ
+        r"""Representation of the operator as a canonical matrix in the computational basis (static method).
+
+        The canonical matrix is the textbook matrix representation that does not consider wires.
+        Implicitly, this assumes that the wires of the operator correspond to the global wire order.
+
+        .. seealso:: :meth:`~.CPhaseShift10.matrix`
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: canonical matrix
+
+        **Example**
+
+        >>> qml.CPhaseShift10.compute_matrix(torch.tensor(0.5))
+            tensor([[1.0+0.0j, 0.0+0.0j, 0.0000+0.0000j, 0.0+0.0j],
+                    [0.0+0.0j, 1.0+0.0j, 0.0000+0.0000j, 0.0+0.0j],
+                    [0.0+0.0j, 0.0+0.0j, 0.8776+0.4794j, 0.0+0.0j],
+                    [0.0+0.0j, 0.0+0.0j, 0.0000+0.0000j, 1.0+0.0j]])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+
+        if qml.math.ndim(phi) > 0:
+            ones = qml.math.ones_like(exp_part)
+            zeros = qml.math.zeros_like(exp_part)
+            matrix = [
+                [ones, zeros, zeros, zeros],
+                [zeros, ones, zeros, zeros],
+                [zeros, zeros, exp_part, zeros],
+                [zeros, zeros, zeros, ones],
+            ]
+
+            return qml.math.stack([stack_last(row) for row in matrix], axis=-2)
+
+        return qml.math.diag([1, 1, exp_part, 1])
+
+    @staticmethod
+    def compute_eigvals(phi):  # pylint: disable=arguments-differ
+        r"""Eigenvalues of the operator in the computational basis (static method).
+
+        If :attr:`diagonalizing_gates` are specified and implement a unitary :math:`U^{\dagger}`,
+        the operator can be reconstructed as
+
+        .. math:: O = U \Sigma U^{\dagger},
+
+        where :math:`\Sigma` is the diagonal matrix containing the eigenvalues.
+
+        Otherwise, no particular order for the eigenvalues is guaranteed.
+
+        .. seealso:: :meth:`~.CPhaseShift10.eigvals`
+
+
+        Args:
+            phi (tensor_like or float): phase shift
+
+        Returns:
+            tensor_like: eigenvalues
+
+        **Example**
+
+        >>> qml.CPhaseShift10.compute_eigvals(torch.tensor(0.5))
+        tensor([1.0000+0.0000j, 1.0000+0.0000j, 0.8776+0.4794j, 1.0000+0.0000j])
+        """
+        if qml.math.get_interface(phi) == "tensorflow":
+            phi = qml.math.cast_like(phi, 1j)
+
+        exp_part = qml.math.exp(1j * phi)
+        ones = qml.math.ones_like(exp_part)
+        return stack_last([ones, ones, exp_part, ones])
+
+    @staticmethod
+    def compute_decomposition(phi, wires):
+        r"""Representation of the operator as a product of other operators (static method). :
+
+        .. math:: O = O_1 O_2 \dots O_n.
+        .. seealso:: :meth:`~.CPhaseShift10.decomposition`.
+
+        Args:
+            phi (float): rotation angle :math:`\phi`
+            wires (Iterable, Wires): wires that the operator acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+
+        **Example:**
+
+        >>> qml.CPhaseShift10.compute_decomposition(1.234, wires=(0,1))
+        [PauliX(wires=[1]),
+        PhaseShift(0.617, wires=[0]),
+        PhaseShift(0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PhaseShift(-0.617, wires=[1]),
+        CNOT(wires=[0, 1]),
+        PauliX(wires=[1])]
+
+        """
+        decomp_ops = [
+            PauliX(wires[1]),
+            PhaseShift(phi / 2, wires=[wires[0]]),
+            PhaseShift(phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PhaseShift(-phi / 2, wires=[wires[1]]),
+            qml.CNOT(wires=wires),
+            PauliX(wires[1]),
+        ]
+        return decomp_ops
+
+    def adjoint(self):
+        return CPhaseShift10(-self.data[0], wires=self.wires)
+
+    def pow(self, z):
+        return [CPhaseShift10(self.data[0] * z, wires=self.wires)]
+
+    @property
+    def control_wires(self):
+        return self.wires[0:1]


Test patch:
diff --git a/tests/ops/op_math/test_controlled.py b/tests/ops/op_math/test_controlled.py
index 590a87e96..54d3a3ddc 100644
--- a/tests/ops/op_math/test_controlled.py
+++ b/tests/ops/op_math/test_controlled.py
@@ -1424,10 +1424,10 @@ def test_ctrl_sanity_check():
     expanded_tape = tape.expand()
 
     expected = [
-        qml.CRX(0.123, wires=[1, 0]),
-        qml.CRY(0.456, wires=[1, 2]),
-        qml.CRX(0.789, wires=[1, 0]),
-        qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]),
+        *qml.CRX(0.123, wires=[1, 0]).decomposition(),
+        *qml.CRY(0.456, wires=[1, 2]).decomposition(),
+        *qml.CRX(0.789, wires=[1, 0]).decomposition(),
+        *qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]).decomposition(),
         qml.CNOT(wires=[1, 2]),
         *qml.CY(wires=[1, 4]).decomposition(),
         *qml.CZ(wires=[1, 0]).decomposition(),
@@ -1458,9 +1458,9 @@ def test_adjoint_of_ctrl():
 
     tape2 = QuantumScript.from_queue(q2)
     expected = [
-        qml.CRZ(4 * onp.pi - 0.456, wires=[5, 0]),
-        qml.CRY(4 * onp.pi - 0.123, wires=[5, 3]),
-        qml.CRX(4 * onp.pi - 0.789, wires=[5, 2]),
+        *qml.CRZ(4 * onp.pi - 0.456, wires=[5, 0]).decomposition(),
+        *qml.CRY(4 * onp.pi - 0.123, wires=[5, 3]).decomposition(),
+        *qml.CRX(4 * onp.pi - 0.789, wires=[5, 2]).decomposition(),
     ]
     for tape in [tape1.expand(depth=1), tape2.expand(depth=1)]:
         for op1, op2 in zip(tape, expected):
@@ -1543,9 +1543,9 @@ def test_ctrl_within_ctrl():
     tape = tape.expand(2, stop_at=lambda op: not isinstance(op, Controlled))
 
     expected = [
-        qml.CRX(0.123, wires=[2, 0]),
+        *qml.CRX(0.123, wires=[2, 0]).decomposition(),
         qml.Toffoli(wires=[2, 0, 1]),
-        qml.CRX(0.456, wires=[2, 0]),
+        *qml.CRX(0.456, wires=[2, 0]).decomposition(),
     ]
     for op1, op2 in zip(tape, expected):
         assert qml.equal(op1, op2)
@@ -1869,10 +1869,10 @@ def test_ctrl_values_sanity_check():
     tape = QuantumScript.from_queue(q_tape)
     expected = [
         qml.PauliX(wires=1),
-        qml.CRX(0.123, wires=[1, 0]),
-        qml.CRY(0.456, wires=[1, 2]),
-        qml.CRX(0.789, wires=[1, 0]),
-        qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]),
+        *qml.CRX(0.123, wires=[1, 0]).decomposition(),
+        *qml.CRY(0.456, wires=[1, 2]).decomposition(),
+        *qml.CRX(0.789, wires=[1, 0]).decomposition(),
+        *qml.CRot(0.111, 0.222, 0.333, wires=[1, 2]).decomposition(),
         qml.CNOT(wires=[1, 2]),
         *qml.CY(wires=[1, 4]).decomposition(),
         *qml.CZ(wires=[1, 0]).decomposition(),
diff --git a/tests/ops/op_math/test_controlled_ops.py b/tests/ops/op_math/test_controlled_ops.py
index e777b6885..2e46fa496 100644
--- a/tests/ops/op_math/test_controlled_ops.py
+++ b/tests/ops/op_math/test_controlled_ops.py
@@ -12,36 +12,40 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 """
-Unit tests for Operators inheriting from Controlled.
+Unit tests for Operators inheriting from ControlledOp.
 """
 
-import copy
+import functools
 
 import numpy as np
 import pytest
 from scipy.linalg import fractional_matrix_power
-from scipy.sparse import csr_matrix
 from scipy.stats import unitary_group
 
-from gate_data import (
-    CY,
-    CZ,
-)
+from gate_data import CY, CZ, CRotx, CRoty, CRotz, CRot3, ControlledPhaseShift
 
 import pennylane as qml
 from pennylane.wires import Wires
-from pennylane.operation import AnyWires
 from pennylane.ops.qubit.matrix_ops import QubitUnitary
 
-# Non-parametrized operations and their matrix representation
 NON_PARAMETRIZED_OPERATIONS = [
     (qml.CY, CY),
     (qml.CZ, CZ),
 ]
 
-SPARSE_MATRIX_SUPPORTED_OPERATIONS = (
-    (qml.CY(wires=[0, 1]), CY),
-    (qml.CZ(wires=[0, 1]), CZ),
+PARAMETRIZED_OPERATIONS = [
+    (qml.CRX, CRotx),
+    (qml.CRY, CRoty),
+    (qml.CRZ, CRotz),
+    (qml.CRot, CRot3),
+    (qml.ControlledPhaseShift, ControlledPhaseShift),
+]
+
+ALL_OPERATIONS = NON_PARAMETRIZED_OPERATIONS + PARAMETRIZED_OPERATIONS
+
+NON_PARAMETRIC_OPS_DECOMPOSITIONS = (
+    (qml.CY, [qml.CRY(np.pi, wires=[0, 1]), qml.S(0)]),
+    (qml.CZ, [qml.ControlledPhaseShift(np.pi, wires=[0, 1])]),
 )
 
 X = np.array([[0, 1], [1, 0]])
@@ -453,176 +457,495 @@ class TestControlledQubitUnitary:
             )
 
 
-class TestOperations:
-    @pytest.mark.parametrize("op_cls, _", NON_PARAMETRIZED_OPERATIONS)
-    def test_nonparametrized_op_copy(self, op_cls, _, tol):
-        """Tests that copied nonparametrized ops function as expected"""
-        op = op_cls(wires=0 if op_cls.num_wires is AnyWires else range(op_cls.num_wires))
-        copied_op = copy.copy(op)
-        assert qml.equal(copied_op, op, atol=tol, rtol=0)
-        assert copied_op is not op
-        assert qml.equal(copied_op, op, atol=tol, rtol=0)
-
-    @pytest.mark.parametrize("ops, mat", NON_PARAMETRIZED_OPERATIONS)
-    def test_matrices(self, ops, mat, tol):
-        """Test matrices of non-parametrized operations are correct"""
-        op = ops(wires=0 if ops.num_wires is AnyWires else range(ops.num_wires))
-        res_static = op.compute_matrix()
-        res_dynamic = op.matrix()
-        assert np.allclose(res_static, mat, atol=tol, rtol=0)
-        assert np.allclose(res_dynamic, mat, atol=tol, rtol=0)
-
-
-class TestDecompositions:  # pylint: disable=too-few-public-methods
-    def test_CY_decomposition(self, tol):
-        """Tests that the decomposition of the CY gate is correct"""
-        op = qml.CY(wires=[0, 1])
-        gate1, gate2 = op.decomposition()
-        decomposed_matrix = qml.matrix(op.decomposition)()
+@pytest.mark.parametrize("op_cls, _", NON_PARAMETRIZED_OPERATIONS)
+def test_map_wires_non_parametric(op_cls, _):
+    """Test get and set private wires in non-parametric controlled operations."""
 
-        assert qml.equal(gate1, qml.CRY(np.pi, wires=[0, 1]))
-        assert qml.equal(gate2, qml.S(0))
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+    op = op_cls(wires=[0, 1])
+    assert op.wires == Wires((0, 1))
 
-    def test_CZ_decomposition(self, tol):
-        """Tests that the decomposition of the CZ gate is correct"""
-        op = qml.CZ(wires=[0, 1])
-        res = op.decomposition()
+    op = op.map_wires(wire_map={0: "a", 1: "b"})
+    assert op.base.wires == Wires("b")
+    assert op.control_wires == Wires("a")
 
-        assert op.has_decomposition
-        assert len(res) == 1
-        assert qml.equal(qml.ctrl(qml.PhaseShift(np.pi, wires=1), 0), res[0], atol=tol, rtol=0)
 
-        decomposed_matrix = res[0].matrix()
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+def test_controlled_phase_shift_alias():
+    """Tests that the alias for ControlledPhaseShift works"""
+    assert qml.equal(qml.ControlledPhaseShift(0.123, wires=[0, 1]), qml.CPhase(0.123, wires=[0, 1]))
 
 
-class TestEigenval:  # pylint: disable=too-few-public-methods
-    def test_CZ_eigenval(self):
-        """Tests that the CZ eigenvalue matches the numpy eigenvalues of the CZ matrix"""
-        op = qml.CZ(wires=[0, 1])
-        exp = np.linalg.eigvals(op.matrix())
-        res = op.eigvals()
-        assert np.allclose(res, exp)
+def _arbitrary_crot(x, y, z):
+    """controlled arbitrary single qubit rotation"""
+    c = np.cos(y / 2)
+    s = np.sin(y / 2)
+    return np.array(
+        [
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, np.exp(-0.5j * (x + z)) * c, -np.exp(0.5j * (x - z)) * s],
+            [0, 0, np.exp(-0.5j * (x - z)) * s, np.exp(0.5j * (x + z)) * c],
+        ]
+    )
 
 
-period_two_ops = (
-    qml.CY(wires=(0, 1)),
-    qml.CZ(wires=(0, 1)),
-)
+def _arbitrary_crot_broadcasted(x, y, z):
+    """controlled arbitrary single qubit rotation broadcasted"""
 
+    c = np.cos(y / 2)
+    s = np.sin(y / 2)
+    return np.array(
+        [
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, np.exp(-0.5j * (_x + _z)) * _c, -np.exp(0.5j * (_x - _z)) * _s],
+                [0, 0, np.exp(-0.5j * (_x - _z)) * _s, np.exp(0.5j * (_x + _z)) * _c],
+            ]
+            for _x, _z, _c, _s in zip(x, z, c, s)
+        ]
+    )
 
-class TestPowMethod:
-    @pytest.mark.parametrize("op", period_two_ops)
-    @pytest.mark.parametrize("n", (1, 5, -1, -5))
-    def test_period_two_pow_odd(self, op, n):
-        """Test that ops with a period of 2 raised to an odd power are the same as the original op."""
-        assert qml.equal(op.pow(n)[0], op)
-        assert np.allclose(op.pow(n)[0].matrix(), op.matrix())
-        assert op.pow(n)[0].name == op.name
-
-    @pytest.mark.parametrize("op", period_two_ops)
-    @pytest.mark.parametrize("n", (2, 6, 0, -2))
-    def test_period_two_pow_even(self, op, n):
-        """Test that ops with a period of 2 raised to an even power are empty lists."""
-        assert len(op.pow(n)) == 0
-
-    @pytest.mark.parametrize("op", period_two_ops)
-    def test_period_two_noninteger_power(self, op):
-        """Test that ops with a period of 2 raised to a non-integer power raise an error."""
-        if isinstance(op, (qml.PauliZ, qml.CZ)):
-            pytest.skip("PauliZ can be raised to any power.")
-        with pytest.raises(qml.operation.PowUndefinedError):
-            op.pow(1.234)
 
-        if op.__class__ is qml.CZ:
-            pytest.skip("CZ can be raised to any power.")
-        with pytest.raises(qml.operation.PowUndefinedError):
-            op.pow(1.234)
+def _phase_shift_matrix_broadcasted(phi):
+    """Phase shift matrix broadcasted"""
 
-    @pytest.mark.parametrize("n", (0.12, -3.462, 3.693))
-    def test_cz_general_power(self, n):
-        """Check that CZ raised to an non-integer power that's not the square root
-        results in a controlled PhaseShift."""
-        op_pow = qml.CZ(wires=[0, 1]).pow(n)
+    expected = np.array([np.eye(4, dtype=complex)] * 3)
+    expected[..., 3, 3] = np.exp(1j * phi)
+    return expected
 
-        assert len(op_pow) == 1
-        assert isinstance(op_pow[0], qml.ops.ControlledOp)
-        assert isinstance(op_pow[0].base, qml.PhaseShift)
-        assert qml.math.allclose(op_pow[0].data[0], np.pi * (n % 2))
 
+EXPECTED_MATRICES = [
+    (qml.CRX, [0], np.identity(4)),
+    (
+        qml.CRX,
+        [np.pi / 2],
+        np.array(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, 1 / np.sqrt(2), -1j / np.sqrt(2)],
+                [0, 0, -1j / np.sqrt(2), 1 / np.sqrt(2)],
+            ]
+        ),
+    ),
+    (qml.CRX, [np.pi], np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1j], [0, 0, -1j, 0]])),
+    (
+        qml.CRX,
+        [np.array([np.pi / 2, np.pi])],
+        np.array(
+            [
+                [
+                    [1, 0, 0, 0],
+                    [0, 1, 0, 0],
+                    [0, 0, 1 / np.sqrt(2), -1j / np.sqrt(2)],
+                    [0, 0, -1j / np.sqrt(2), 1 / np.sqrt(2)],
+                ],
+                [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1j], [0, 0, -1j, 0]],
+            ]
+        ),
+    ),
+    (qml.CRY, [0], np.identity(4)),
+    (
+        qml.CRY,
+        [np.pi / 2],
+        np.array(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, 1 / np.sqrt(2), -1 / np.sqrt(2)],
+                [0, 0, 1 / np.sqrt(2), 1 / np.sqrt(2)],
+            ]
+        ),
+    ),
+    (qml.CRY, [np.pi], np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]])),
+    (
+        qml.CRY,
+        [np.array([np.pi / 2, np.pi])],
+        np.array(
+            [
+                [
+                    [1, 0, 0, 0],
+                    [0, 1, 0, 0],
+                    [0, 0, 1 / np.sqrt(2), -1 / np.sqrt(2)],
+                    [0, 0, 1 / np.sqrt(2), 1 / np.sqrt(2)],
+                ],
+                [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]],
+            ]
+        ),
+    ),
+    (qml.CRZ, [0], np.identity(4)),
+    (
+        qml.CRZ,
+        [np.pi / 2],
+        np.array(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, np.exp(-1j * np.pi / 4), 0],
+                [0, 0, 0, np.exp(1j * np.pi / 4)],
+            ]
+        ),
+    ),
+    (qml.CRZ, [np.pi], np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1j, 0], [0, 0, 0, 1j]])),
+    (
+        qml.CRZ,
+        [np.array([np.pi / 2, np.pi])],
+        np.array(
+            [
+                [
+                    [1, 0, 0, 0],
+                    [0, 1, 0, 0],
+                    [0, 0, np.exp(-1j * np.pi / 4), 0],
+                    [0, 0, 0, np.exp(1j * np.pi / 4)],
+                ],
+                [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1j, 0], [0, 0, 0, 1j]],
+            ]
+        ),
+    ),
+    (qml.CRot, [0, 0, 0], np.identity(4)),
+    (
+        qml.CRot,
+        [np.pi, np.pi, np.pi],
+        np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]]),
+    ),
+    (qml.CRot, [0.432, -0.152, 0.9234], _arbitrary_crot(0.432, -0.152, 0.9234)),
+    (qml.CRot, [np.zeros(5), np.zeros(5), np.zeros(5)], np.identity(4)),
+    (
+        qml.CRot,
+        [np.ones(3) * np.pi, np.ones(3) * np.pi, np.ones(3) * np.pi],
+        np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]]),
+    ),
+    (
+        qml.CRot,
+        [np.array([0.432, -0.124]), np.array([-0.152, 2.912]), np.array([0.9234, -9.2])],
+        _arbitrary_crot_broadcasted(
+            np.array([0.432, -0.124]), np.array([-0.152, 2.912]), np.array([0.9234, -9.2])
+        ),
+    ),
+    (qml.ControlledPhaseShift, [0.123], ControlledPhaseShift(0.123)),
+    (
+        qml.ControlledPhaseShift,
+        [np.array([0.2, np.pi / 2, -0.1])],
+        _phase_shift_matrix_broadcasted(np.array([0.2, np.pi / 2, -0.1])),
+    ),
+]
 
-class TestControlledMethod:  # pylint: disable=too-few-public-methods
-    """Tests for the _controlled method of non-parametric operations."""
+EXPECTED_EIGVALS = [
+    (qml.CRZ, [0], np.ones(4)),
+    (qml.CRZ, [np.pi / 2], np.array([1, 1, np.exp(-1j * np.pi / 4), np.exp(1j * np.pi / 4)])),
+    (qml.CRZ, [np.pi], np.array([1, 1, -1j, 1j])),
+    (
+        qml.CRZ,
+        [np.array([np.pi / 2, np.pi])],
+        [
+            np.array([1, 1, np.exp(-1j * np.pi / 4), np.exp(1j * np.pi / 4)]),
+            np.array([1, 1, -1j, 1j]),
+        ],
+    ),
+    (qml.ControlledPhaseShift, [0.123], np.linalg.eigvals(ControlledPhaseShift(0.123))),
+    (
+        qml.ControlledPhaseShift,
+        [np.array([0.2, np.pi / 2, -0.1])],
+        np.linalg.eigvals(_phase_shift_matrix_broadcasted(np.array([0.2, np.pi / 2, -0.1]))),
+    ),
+]
 
-    def test_CZ(self):
-        """Test the PauliZ _controlled method."""
-        out = qml.CZ(wires=[0, 1])._controlled("a")  # pylint: disable=protected-access
-        assert qml.equal(out, qml.CCZ(("a", 0, 1)))
 
+class TestComputations:
+    """Tests computing the matrix and eigenvalues of a controlled operation"""
 
-class TestSparseMatrix:  # pylint: disable=too-few-public-methods
-    @pytest.mark.parametrize("op, mat", SPARSE_MATRIX_SUPPORTED_OPERATIONS)
-    def test_sparse_matrix(self, op, mat):
-        """Tests the sparse matrix method for operations which support it."""
-        expected_sparse_mat = csr_matrix(mat)
-        sparse_mat = op.sparse_matrix()
+    @pytest.mark.parametrize("op, params, expected_matrix", EXPECTED_MATRICES)
+    def test_matrix(self, tol, op, params, expected_matrix):
+        """Tests that the correct matrix is returned"""
 
-        assert isinstance(sparse_mat, csr_matrix)
-        assert isinstance(expected_sparse_mat, csr_matrix)
-        assert all(sparse_mat.data == expected_sparse_mat.data)
-        assert all(sparse_mat.indices == expected_sparse_mat.indices)
+        assert np.allclose(op.compute_matrix(*params), expected_matrix, atol=tol, rtol=0)
+        assert np.allclose(op(*params, wires=[0, 1]).matrix(), expected_matrix, atol=tol, rtol=0)
 
+    @pytest.mark.tf
+    @pytest.mark.parametrize("op, params, expected_matrix", EXPECTED_MATRICES)
+    def test_matrix_tf(self, tol, op, params, expected_matrix):
+        """Tests that the correct matrix is returned when using TensorFlow"""
 
-label_data = [
-    (qml.CY(wires=(0, 1)), "Y"),
-    (qml.CZ(wires=(0, 1)), "Z"),
-]
+        import tensorflow as tf
 
+        variables = [tf.Variable(param) for param in params]
 
-@pytest.mark.parametrize("op, label", label_data)
-def test_label_method(op, label):
-    """Tests that the label method gives the expected result."""
-    assert op.label() == label
-    assert op.label(decimals=2) == label
+        assert np.allclose(op.compute_matrix(*variables), expected_matrix, atol=tol, rtol=0)
+        assert np.allclose(op(*variables, wires=[0, 1]).matrix(), expected_matrix, atol=tol, rtol=0)
 
+    @pytest.mark.parametrize("op, params, expected_eigvals", EXPECTED_EIGVALS)
+    def test_eigvals(self, tol, op, params, expected_eigvals):
+        """Tests that the correct eigenvalues are returned"""
 
-control_data = [
-    (qml.CY(wires=(0, 1)), Wires(0)),
-    (qml.CZ(wires=(0, 1)), Wires(0)),
-]
+        assert np.allclose(op.compute_eigvals(*params), expected_eigvals, atol=tol, rtol=0)
+        assert np.allclose(op(*params, wires=[0, 1]).eigvals(), expected_eigvals, atol=tol, rtol=0)
 
+    @pytest.mark.tf
+    @pytest.mark.parametrize("op, params, expected_eigvals", EXPECTED_EIGVALS)
+    def test_eigvals_tf(self, tol, op, params, expected_eigvals):
+        """Tests that the correct eigenvalues are returned when using TensorFlow"""
 
-@pytest.mark.parametrize("op, control_wires", control_data)
-def test_control_wires(op, control_wires):
-    """Test ``control_wires`` attribute for non-parametrized operations."""
+        import tensorflow as tf
 
-    assert op.control_wires == control_wires
+        variables = [tf.Variable(param) for param in params]
 
+        assert np.allclose(op.compute_eigvals(*variables), expected_eigvals, atol=tol, rtol=0)
+        assert np.allclose(
+            op(*variables, wires=[0, 1]).eigvals(), expected_eigvals, atol=tol, rtol=0
+        )
 
-involution_ops = [  # ops who are their own inverses
-    qml.CY((0, 1)),
-    qml.CZ(wires=(0, 1)),
-]
 
+def _dot_broadcasted(a, b):
+    return np.einsum("...ij,...jk->...ik", a, b)
 
-@pytest.mark.parametrize("op", involution_ops)
-def test_adjoint_method(op):
-    """Tests the adjoint method for operations that are their own adjoint."""
-    adj_op = copy.copy(op)
-    for _ in range(4):
-        adj_op = adj_op.adjoint()
 
-        assert qml.equal(adj_op, op)
+def _multi_dot_broadcasted(matrices):
+    return functools.reduce(_dot_broadcasted, matrices)
 
 
-@pytest.mark.parametrize("op_cls, _", NON_PARAMETRIZED_OPERATIONS)
-def test_map_wires(op_cls, _):
-    """Test that we can get and set private wires in all operations."""
+class TestDecompositions:
+    @pytest.mark.parametrize("phi", [0.432, np.array([0.1, 2.1])])
+    def test_CRX(self, phi):
+        """Test the decomposition for CRX."""
 
-    op = op_cls(wires=[0, 1])
-    assert op.wires == Wires((0, 1))
+        ops1 = qml.CRX.compute_decomposition(phi, wires=[0, 1])
+        ops2 = qml.CRX(phi, wires=(0, 1)).decomposition()
 
-    op = op.map_wires(wire_map={0: "a", 1: "b"})
-    assert op.base.wires == Wires(("b"))
-    assert op.control_wires == Wires(("a"))
+        expected_ops = [
+            qml.RZ(np.pi / 2, 1),
+            qml.RY(phi / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+            qml.RY(-phi / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+            qml.RZ(-np.pi / 2, 1),
+        ]
+
+        assert ops1 == expected_ops
+        assert ops2 == expected_ops
+
+    @pytest.mark.parametrize("phi", [0.432, np.array([2.1, 0.2])])
+    def test_CRY(self, phi):
+        """Test the decomposition for CRY."""
+
+        ops1 = qml.CRY.compute_decomposition(phi, wires=[0, 1])
+        ops2 = qml.CRY(phi, wires=(0, 1)).decomposition()
+
+        expected_ops = [
+            qml.RY(phi / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+            qml.RY(-phi / 2, 1),
+            qml.CNOT(wires=(0, 1)),
+        ]
+
+        assert ops1 == expected_ops
+        assert ops2 == expected_ops
+
+    @pytest.mark.parametrize("phi", [0.321, np.array([0.6, 2.1])])
+    def test_CRZ(self, phi):
+        """Test the decomposition for CRZ."""
+
+        ops1 = qml.CRZ.compute_decomposition(phi, wires=[1, 0])
+        ops2 = qml.CRZ(phi, wires=(1, 0)).decomposition()
+
+        expected_ops = [
+            qml.PhaseShift(phi / 2, wires=0),
+            qml.CNOT(wires=[1, 0]),
+            qml.PhaseShift(-phi / 2, wires=0),
+            qml.CNOT(wires=[1, 0]),
+        ]
+
+        assert ops1 == expected_ops
+        assert ops2 == expected_ops
+
+    @pytest.mark.parametrize("phi, theta, omega", [[0.5, 0.6, 0.7], [0.1, -0.4, 0.7], [-10, 5, -1]])
+    def test_CRot(self, tol, phi, theta, omega):
+        """Tests that the decomposition of the CRot gate is correct"""
+
+        op = qml.CRot(phi, theta, omega, wires=[0, 1])
+        res = op.decomposition()
+
+        mats = []
+        for i in reversed(res):
+            if len(i.wires) == 1:
+                mats.append(np.kron(np.eye(2), i.matrix()))
+            else:
+                mats.append(i.matrix())
+
+        decomposed_matrix = np.linalg.multi_dot(mats)
+        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+
+    @pytest.mark.parametrize(
+        "phi, theta, omega",
+        [
+            [np.array([0.1, 0.2]), np.array([-0.4, 2.19]), np.array([0.7, -0.7])],
+            [np.array([0.1, 0.2, 0.9]), -0.4, np.array([0.7, 0.0, -0.7])],
+        ],
+    )
+    def test_CRot_broadcasted(self, tol, phi, theta, omega):
+        """Tests that the decomposition of the broadcasted CRot gate is correct"""
+
+        op = qml.CRot(phi, theta, omega, wires=[0, 1])
+        res = op.decomposition()
+
+        mats = []
+        for i in reversed(res):
+            mat = i.matrix()
+            if len(i.wires) == 1:
+                I = np.eye(2)[np.newaxis] if qml.math.ndim(mat) == 3 else np.eye(2)
+                mats.append(np.kron(I, mat))
+            else:
+                mats.append(mat)
+
+        decomposed_matrix = _multi_dot_broadcasted(mats)
+        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+
+    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
+    def test_controlled_phase_shift(self, phi):
+        """Tests that the ControlledPhaseShift calculates the correct decomposition"""
+
+        op = qml.ControlledPhaseShift(phi, wires=[0, 2])
+        decomp = op.decomposition()
+
+        mats = []
+        for i in reversed(decomp):
+            if i.wires.tolist() == [0]:
+                mats.append(np.kron(i.matrix(), np.eye(4)))
+            elif i.wires.tolist() == [1]:
+                mats.append(np.kron(np.eye(2), np.kron(i.matrix(), np.eye(2))))
+            elif i.wires.tolist() == [2]:
+                mats.append(np.kron(np.eye(4), i.matrix()))
+            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
+                mats.append(np.kron(i.matrix(), np.eye(2)))
+            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
+                mats.append(
+                    np.array(
+                        [
+                            [1, 0, 0, 0, 0, 0, 0, 0],
+                            [0, 1, 0, 0, 0, 0, 0, 0],
+                            [0, 0, 1, 0, 0, 0, 0, 0],
+                            [0, 0, 0, 1, 0, 0, 0, 0],
+                            [0, 0, 0, 0, 0, 1, 0, 0],
+                            [0, 0, 0, 0, 1, 0, 0, 0],
+                            [0, 0, 0, 0, 0, 0, 0, 1],
+                            [0, 0, 0, 0, 0, 0, 1, 0],
+                        ]
+                    )
+                )
+
+        decomposed_matrix = np.linalg.multi_dot(mats)
+        lam = np.exp(1j * phi)
+        exp = np.array(
+            [
+                [1, 0, 0, 0, 0, 0, 0, 0],
+                [0, 1, 0, 0, 0, 0, 0, 0],
+                [0, 0, 1, 0, 0, 0, 0, 0],
+                [0, 0, 0, 1, 0, 0, 0, 0],
+                [0, 0, 0, 0, 1, 0, 0, 0],
+                [0, 0, 0, 0, 0, lam, 0, 0],
+                [0, 0, 0, 0, 0, 0, 1, 0],
+                [0, 0, 0, 0, 0, 0, 0, lam],
+            ]
+        )
+
+        assert np.allclose(decomposed_matrix, exp)
+
+    def test_controlled_phase_shift_broadcasted(self):
+        """Tests that the ControlledPhaseShift calculates the correct decomposition"""
+
+        phi = np.array([-0.2, 4.2, 1.8])
+        op = qml.ControlledPhaseShift(phi, wires=[0, 2])
+        decomp = op.decomposition()
+
+        mats = []
+        for i in reversed(decomp):
+            mat = i.matrix()
+            eye = np.eye(2)[np.newaxis] if np.ndim(mat) == 3 else np.eye(2)
+            if i.wires.tolist() == [0]:
+                mats.append(np.kron(mat, np.kron(eye, eye)))
+            elif i.wires.tolist() == [1]:
+                mats.append(np.kron(eye, np.kron(mat, eye)))
+            elif i.wires.tolist() == [2]:
+                mats.append(np.kron(np.kron(eye, eye), mat))
+            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
+                mats.append(np.kron(mat, eye))
+            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
+                mats.append(
+                    np.array(
+                        [
+                            [1, 0, 0, 0, 0, 0, 0, 0],
+                            [0, 1, 0, 0, 0, 0, 0, 0],
+                            [0, 0, 1, 0, 0, 0, 0, 0],
+                            [0, 0, 0, 1, 0, 0, 0, 0],
+                            [0, 0, 0, 0, 0, 1, 0, 0],
+                            [0, 0, 0, 0, 1, 0, 0, 0],
+                            [0, 0, 0, 0, 0, 0, 0, 1],
+                            [0, 0, 0, 0, 0, 0, 1, 0],
+                        ]
+                    )
+                )
+
+        decomposed_matrix = _multi_dot_broadcasted(mats)
+        lam = np.exp(1j * phi)
+        exp = np.array([np.diag([1, 1, 1, 1, 1, el, 1, el]) for el in lam])
+
+        assert np.allclose(decomposed_matrix, exp)
+
+    @pytest.mark.parametrize("ops, expected_ops", NON_PARAMETRIC_OPS_DECOMPOSITIONS)
+    def test_non_parametric_op_decompositions(self, ops, expected_ops, tol):
+        """Tests that decompositions of non-parametrized operations are correct"""
+
+        op = ops(wires=[0, 1])
+        decomps = op.decomposition()
+        decomposed_matrix = qml.matrix(op.decomposition)()
+
+        for gate, expected in zip(decomps, expected_ops):
+            assert qml.equal(gate, expected, atol=tol, rtol=0)
+
+        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
+
+
+def test_simplify_crot():
+    """Simplify CRot operations with different parameters."""
+
+    crot_x = qml.CRot(np.pi / 2, 0.1, -np.pi / 2, wires=[0, 1])
+    simplify_crot_x = crot_x.simplify()
+
+    assert simplify_crot_x.name == "CRX"
+    assert simplify_crot_x.data == (0.1,)
+    assert np.allclose(simplify_crot_x.matrix(), crot_x.matrix())
+
+    crot_y = qml.CRot(0, 0.1, 0, wires=[0, 1])
+    simplify_crot_y = crot_y.simplify()
+
+    assert simplify_crot_y.name == "CRY"
+    assert simplify_crot_y.data == (0.1,)
+    assert np.allclose(simplify_crot_y.matrix(), crot_y.matrix())
+
+    crot_z = qml.CRot(0.1, 0, 0.2, wires=[0, 1])
+    simplify_crot_z = crot_z.simplify()
+
+    assert simplify_crot_z.name == "CRZ"
+    assert np.allclose(simplify_crot_z.data, [0.3])
+    assert np.allclose(simplify_crot_z.matrix(), crot_z.matrix())
+
+    crot = qml.CRot(0.1, 0.2, 0.3, wires=[0, 1])
+    not_simplified_crot = crot.simplify()
+
+    assert not_simplified_crot.name == "CRot"
+    assert np.allclose(not_simplified_crot.matrix(), crot.matrix())
+
+
+controlled_data = [
+    (qml.RX(1.234, wires=0), qml.CRX(1.234, wires=("a", 0))),
+    (qml.RY(1.234, wires=0), qml.CRY(1.234, wires=("a", 0))),
+    (qml.RZ(1.234, wires=0), qml.CRZ(1.234, wires=("a", 0))),
+    (qml.PhaseShift(1.234, wires=0), qml.ControlledPhaseShift(1.234, wires=("a", 0))),
+    (qml.Rot(1.2, 2.3, 3.4, wires=0), qml.CRot(1.2, 2.3, 3.4, wires=("a", 0))),
+]
+
+
+@pytest.mark.parametrize("base, cbase", controlled_data)
+def test_controlled_method(base, cbase):
+    """Tests the _controlled method for parametric ops."""
+    # pylint: disable=protected-access
+    assert qml.equal(base._controlled("a"), cbase)
diff --git a/tests/ops/qubit/test_non_parametric_ops.py b/tests/ops/qubit/test_non_parametric_ops.py
index 775de28f4..56edb2717 100644
--- a/tests/ops/qubit/test_non_parametric_ops.py
+++ b/tests/ops/qubit/test_non_parametric_ops.py
@@ -38,6 +38,8 @@ from gate_data import (
     Y,
     Z,
     CH,
+    CY,
+    CZ,
 )
 from scipy.sparse import csr_matrix
 from scipy.stats import unitary_group
@@ -49,31 +51,26 @@ from pennylane.wires import Wires
 # Non-parametrized operations and their matrix representation
 NON_PARAMETRIZED_OPERATIONS = [
     (qml.Identity, I),
-    (qml.CNOT, CNOT),
     (qml.SWAP, SWAP),
     (qml.ISWAP, ISWAP),
     (qml.SISWAP, SISWAP),
-    (qml.CCZ, CCZ),
     (qml.S, S),
     (qml.T, T),
-    (qml.CSWAP, CSWAP),
-    (qml.Toffoli, Toffoli),
     (qml.ECR, ECR),
+    # Controlled operations
+    (qml.CNOT, CNOT),
     (qml.CH, CH),
+    (qml.Toffoli, Toffoli),
+    (qml.CSWAP, CSWAP),
+    (qml.CY, CY),
+    (qml.CZ, CZ),
+    (qml.CCZ, CCZ),
 ]
 
-SPARSE_MATRIX_SUPPORTED_OPERATIONS = (
-    (qml.Identity(wires=0), I),
-    (qml.Hadamard(wires=0), H),
-    (qml.PauliZ(wires=0), Z),
-    (qml.PauliX(wires=0), X),
-    (qml.PauliY(wires=0), Y),
-)
-
 
 class TestOperations:
     @pytest.mark.parametrize("op_cls, _", NON_PARAMETRIZED_OPERATIONS)
-    def test_nonparametrized_op_copy(self, op_cls, _, tol):
+    def test_op_copy(self, op_cls, _, tol):
         """Tests that copied nonparametrized ops function as expected"""
         op = op_cls(wires=0 if op_cls.num_wires is AnyWires else range(op_cls.num_wires))
         copied_op = copy.copy(op)
@@ -90,6 +87,8 @@ class TestOperations:
 
 
 class TestDecompositions:
+    """Tests that the decomposition of non-parametrized operations is correct"""
+
     def test_x_decomposition(self, tol):
         """Tests that the decomposition of the PauliX is correct"""
         op = qml.PauliX(wires=0)
@@ -537,6 +536,14 @@ class TestEigenval:
         expected = np.linalg.eigvals(qml.SX(wires=0).matrix())
         assert np.allclose(evals, expected)
 
+    def test_cz_eigenval(self):
+        """Tests that the CZ eigenvalue matches the numpy eigenvalues of the CZ matrix"""
+
+        op = qml.CZ(wires=[0, 1])
+        exp = np.linalg.eigvals(op.matrix())
+        res = op.eigvals()
+        assert np.allclose(res, exp)
+
 
 class TestMultiControlledX:
     """Tests for the MultiControlledX"""
@@ -937,11 +944,14 @@ period_two_ops = (
     qml.PauliY(0),
     qml.PauliZ(0),
     qml.Hadamard("a"),
-    qml.CNOT(wires=(0, 1)),
-    qml.CH(wires=(0, 1)),
     qml.SWAP(wires=(0, 1)),
     qml.ISWAP(wires=(0, 1)),
     qml.ECR(wires=(0, 1)),
+    # Controlled operations
+    qml.CNOT(wires=(0, 1)),
+    qml.CY(wires=(0, 1)),
+    qml.CZ(wires=(0, 1)),
+    qml.CH(wires=(0, 1)),
     qml.CCZ(wires=(0, 1, 2)),
     qml.CSWAP(wires=(0, 1, 2)),
     qml.Toffoli(wires=(0, 1, 2)),
@@ -965,11 +975,22 @@ class TestPowMethod:
     @pytest.mark.parametrize("op", period_two_ops)
     def test_period_two_noninteger_power(self, op):
         """Test that ops with a period of 2 raised to a non-integer power raise an error."""
-        if op.__class__ is qml.PauliZ:
+        if op.__class__ in [qml.PauliZ, qml.CZ]:
             pytest.skip("PauliZ can be raised to any power.")
         with pytest.raises(qml.operation.PowUndefinedError):
             op.pow(1.234)
 
+    @pytest.mark.parametrize("n", (0.12, -3.462, 3.693))
+    def test_cz_general_power(self, n):
+        """Check that CZ raised to a non-integer power that's not the square root
+        results in a controlled PhaseShift."""
+        op_pow = qml.CZ(wires=[0, 1]).pow(n)
+
+        assert len(op_pow) == 1
+        assert isinstance(op_pow[0], qml.ops.ControlledOp)
+        assert isinstance(op_pow[0].base, qml.PhaseShift)
+        assert qml.math.allclose(op_pow[0].data[0], np.pi * (n % 2))
+
     @pytest.mark.parametrize("n", (0.5, 2.5, -1.5))
     def test_paulix_squareroot(self, n):
         """Check that the square root of PauliX is SX"""
@@ -1130,16 +1151,31 @@ class TestControlledMethod:
         out = original._controlled("a")
         assert qml.equal(original, out)
 
+    def test_CZ(self):
+        """Test the PauliZ _controlled method."""
+        out = qml.CZ(wires=[0, 1])._controlled("a")
+        assert qml.equal(out, qml.CCZ(("a", 0, 1)))
 
-class TestSparseMatrix:
-    @pytest.mark.parametrize("op, mat", SPARSE_MATRIX_SUPPORTED_OPERATIONS)
-    def test_sparse_matrix(self, op, mat):
-        expected_sparse_mat = csr_matrix(mat)
-        sparse_mat = op.sparse_matrix()
 
-        assert isinstance(sparse_mat, type(expected_sparse_mat))
-        assert all(sparse_mat.data == expected_sparse_mat.data)
-        assert all(sparse_mat.indices == expected_sparse_mat.indices)
+SPARSE_MATRIX_SUPPORTED_OPERATIONS = (
+    (qml.Identity(wires=0), I),
+    (qml.Hadamard(wires=0), H),
+    (qml.PauliZ(wires=0), Z),
+    (qml.PauliX(wires=0), X),
+    (qml.PauliY(wires=0), Y),
+    (qml.CY(wires=[0, 1]), CY),
+    (qml.CZ(wires=[0, 1]), CZ),
+)
+
+
+@pytest.mark.parametrize("op, mat", SPARSE_MATRIX_SUPPORTED_OPERATIONS)
+def test_sparse_matrix(op, mat):
+    expected_sparse_mat = csr_matrix(mat)
+    sparse_mat = op.sparse_matrix()
+
+    assert isinstance(sparse_mat, type(expected_sparse_mat))
+    assert all(sparse_mat.data == expected_sparse_mat.data)
+    assert all(sparse_mat.indices == expected_sparse_mat.indices)
 
 
 label_data = [
@@ -1151,19 +1187,22 @@ label_data = [
     (qml.S(wires=0), "S"),
     (qml.T(wires=0), "T"),
     (qml.SX(wires=0), "SX"),
-    (qml.CNOT(wires=(0, 1)), "X"),
-    (qml.CH(wires=(0, 1)), "H"),
     (qml.SWAP(wires=(0, 1)), "SWAP"),
     (qml.ISWAP(wires=(0, 1)), "ISWAP"),
     (qml.ECR(wires=(0, 1)), "ECR"),
     (qml.SISWAP(wires=(0, 1)), "SISWAP"),
     (qml.SQISW(wires=(0, 1)), "SISWAP"),
+    (qml.Barrier(0), "||"),
+    (qml.WireCut(wires=0), "//"),
+    # Controlled operations
+    (qml.CY(wires=(0, 1)), "Y"),
+    (qml.CZ(wires=(0, 1)), "Z"),
+    (qml.CNOT(wires=(0, 1)), "X"),
+    (qml.CH(wires=(0, 1)), "H"),
     (qml.CCZ(wires=(0, 1, 2)), "Z"),
     (qml.CSWAP(wires=(0, 1, 2)), "SWAP"),
     (qml.Toffoli(wires=(0, 1, 2)), "X"),
     (qml.MultiControlledX(wires=(0, 1, 2, 3)), "X"),
-    (qml.Barrier(0), "||"),
-    (qml.WireCut(wires=0), "//"),
 ]
 
 
@@ -1185,6 +1224,9 @@ control_data = [
     (qml.ISWAP(wires=(0, 1)), Wires([])),
     (qml.SISWAP(wires=(0, 1)), Wires([])),
     (qml.ECR(wires=(0, 1)), Wires([])),
+    # Controlled operations
+    (qml.CY(wires=(0, 1)), Wires(0)),
+    (qml.CZ(wires=(0, 1)), Wires(0)),
     (qml.CNOT(wires=(0, 1)), Wires(0)),
     (qml.CH(wires=(0, 1)), Wires(0)),
     (qml.CSWAP(wires=(0, 1, 2)), Wires([0])),
@@ -1197,7 +1239,6 @@ control_data = [
 @pytest.mark.parametrize("op, control_wires", control_data)
 def test_control_wires(op, control_wires):
     """Test ``control_wires`` attribute for non-parametrized operations."""
-
     assert op.control_wires == control_wires
 
 
@@ -1207,21 +1248,24 @@ involution_ops = [  # ops who are their own inverses
     qml.PauliX(0),
     qml.PauliY(0),
     qml.PauliZ(0),
-    qml.CNOT((0, 1)),
-    qml.CH((0, 1)),
     qml.SWAP((0, 1)),
     qml.ECR((0, 1)),
+    qml.Barrier(0),
+    qml.WireCut(0),
+    # Controlled operations
+    qml.CNOT((0, 1)),
+    qml.CH((0, 1)),
+    qml.CY((0, 1)),
+    qml.CZ(wires=(0, 1)),
     qml.CSWAP((0, 1, 2)),
     qml.CCZ((0, 1, 2)),
     qml.Toffoli((0, 1, 2)),
     qml.MultiControlledX(wires=(0, 1, 2, 3)),
-    qml.Barrier(0),
-    qml.WireCut(0),
 ]
 
 
 @pytest.mark.parametrize("op", involution_ops)
-def test_adjoint_method(op):
+def test_involution_operators(op):
     adj_op = copy.copy(op)
     for _ in range(4):
         adj_op = adj_op.adjoint()
diff --git a/tests/ops/qubit/test_parametric_ops.py b/tests/ops/qubit/test_parametric_ops.py
index 854805de8..cfddda99d 100644
--- a/tests/ops/qubit/test_parametric_ops.py
+++ b/tests/ops/qubit/test_parametric_ops.py
@@ -20,14 +20,13 @@ from functools import reduce
 
 import numpy as np
 import pytest
-from gate_data import ControlledPhaseShift, CPhaseShift00, CPhaseShift01, CPhaseShift10, Z
+from gate_data import CPhaseShift00, CPhaseShift01, CPhaseShift10, Z
 from scipy import sparse
 
 import pennylane as qml
 from pennylane import numpy as npp
 from pennylane.ops.qubit import (
     RX as old_loc_RX,
-    ControlledPhaseShift as old_loc_ControlledPhaseShift,
     MultiRZ as old_loc_MultiRZ,
 )
 from pennylane.wires import Wires
@@ -109,7 +108,6 @@ BROADCASTED_OPERATIONS = [
     qml.DiagonalQubitUnitary(np.array([[1.0, 1.0j], [1.0j, 1.0j]]), wires=1),
 ]
 
-
 NON_PARAMETRIZED_OPERATIONS = [
     qml.Identity(0),
     qml.S(wires=0),
@@ -319,150 +317,6 @@ class TestDecompositions:
                 assert isinstance(op, c)
                 assert op.parameters == p
 
-    def test_CRX_decomposition(self):
-        """Test the decomposition for CRX."""
-        phi = 0.432
-
-        ops1 = qml.CRX.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRX(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.RZ, qml.RY, qml.CNOT, qml.RY, qml.CNOT, qml.RZ]
-        params = [[np.pi / 2], [phi / 2], [], [-phi / 2], [], [-np.pi / 2]]
-        wires = [Wires(1), Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1)), Wires(1)]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert op.parameters == p
-                assert op.wires == w
-
-    def test_CRX_decomposition_broadcasted(self):
-        """Test the decomposition for broadcasted CRX."""
-        phi = np.array([0.1, 2.1])
-
-        ops1 = qml.CRX.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRX(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.RZ, qml.RY, qml.CNOT, qml.RY, qml.CNOT, qml.RZ]
-        params = [[np.pi / 2], [phi / 2], [], [-phi / 2], [], [-np.pi / 2]]
-        wires = [Wires(1), Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1)), Wires(1)]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert qml.math.allclose(op.parameters, p)
-                assert op.wires == w
-
-    def test_CRY_decomposition(self):
-        """Test the decomposition for CRY."""
-        phi = 0.432
-
-        ops1 = qml.CRY.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRY(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.RY, qml.CNOT, qml.RY, qml.CNOT]
-        params = [[phi / 2], [], [-phi / 2], []]
-        wires = [Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1))]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert np.allclose(op.parameters, p)
-                assert op.wires == w
-
-    def test_CRY_decomposition_broadcasted(self):
-        """Test the decomposition for broadcastedCRY."""
-        phi = np.array([2.1, 0.2])
-
-        ops1 = qml.CRY.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRY(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.RY, qml.CNOT, qml.RY, qml.CNOT]
-        params = [[phi / 2], [], [-phi / 2], []]
-        wires = [Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1))]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert np.allclose(op.parameters, p)
-                assert op.wires == w
-
-    def test_CRZ_decomposition(self):
-        """Test the decomposition for CRZ."""
-        phi = 0.321
-
-        ops1 = qml.CRZ.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRZ(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.PhaseShift, qml.CNOT, qml.PhaseShift, qml.CNOT]
-        params = [[phi / 2], [], [-phi / 2], []]
-        wires = [Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1))]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert np.allclose(op.parameters, p)
-                assert op.wires == w
-
-    def test_CRZ_decomposition_broadcasted(self):
-        """Test the decomposition for broadcasted CRZ."""
-        phi = np.array([0.6, 2.1])
-
-        ops1 = qml.CRZ.compute_decomposition(phi, wires=[0, 1])
-        ops2 = qml.CRZ(phi, wires=(0, 1)).decomposition()
-
-        classes = [qml.PhaseShift, qml.CNOT, qml.PhaseShift, qml.CNOT]
-        params = [[phi / 2], [], [-phi / 2], []]
-        wires = [Wires(1), Wires((0, 1)), Wires(1), Wires((0, 1))]
-
-        for ops in [ops1, ops2]:
-            for op, c, p, w in zip(ops, classes, params, wires):
-                assert isinstance(op, c)
-                assert np.allclose(op.parameters, p)
-                assert op.wires == w
-
-    @pytest.mark.parametrize("phi, theta, omega", [[0.5, 0.6, 0.7], [0.1, -0.4, 0.7], [-10, 5, -1]])
-    def test_CRot_decomposition(self, tol, phi, theta, omega):
-        """Tests that the decomposition of the CRot gate is correct"""
-        op = qml.CRot(phi, theta, omega, wires=[0, 1])
-        res = op.decomposition()
-
-        mats = []
-        for i in reversed(res):
-            if len(i.wires) == 1:
-                mats.append(np.kron(np.eye(2), i.matrix()))
-            else:
-                mats.append(i.matrix())
-
-        decomposed_matrix = np.linalg.multi_dot(mats)
-
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
-
-    @pytest.mark.parametrize(
-        "phi, theta, omega",
-        [
-            [np.array([0.1, 0.2]), np.array([-0.4, 2.19]), np.array([0.7, -0.7])],
-            [np.array([0.1, 0.2, 0.9]), -0.4, np.array([0.7, 0.0, -0.7])],
-        ],
-    )
-    def test_CRot_decomposition_broadcasted(self, tol, phi, theta, omega):
-        """Tests that the decomposition of the broadcasted CRot gate is correct"""
-        op = qml.CRot(phi, theta, omega, wires=[0, 1])
-        res = op.decomposition()
-
-        mats = []
-        for i in reversed(res):
-            mat = i.matrix()
-            if len(i.wires) == 1:
-                I = np.eye(2)[np.newaxis] if qml.math.ndim(mat) == 3 else np.eye(2)
-                mats.append(np.kron(I, mat))
-            else:
-                mats.append(mat)
-
-        decomposed_matrix = multi_dot_broadcasted(mats)
-
-        assert np.allclose(decomposed_matrix, op.matrix(), atol=tol, rtol=0)
-
     def test_U1_decomposition(self):
         """Test the decomposition for U1."""
         phi = 0.432
@@ -815,99 +669,6 @@ class TestDecompositions:
         for expected_mat, decomp_mat in zip(expected_mats, decomp_mats):
             assert np.allclose(expected_mat, decomp_mat)
 
-    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_decomp(self, phi, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation
-        calculates the correct decomposition"""
-        op = cphase_op(phi, wires=[0, 2])
-        decomp = op.decomposition()
-
-        mats = []
-        for i in reversed(decomp):
-            if i.wires.tolist() == [0]:
-                mats.append(np.kron(i.matrix(), np.eye(4)))
-            elif i.wires.tolist() == [1]:
-                mats.append(np.kron(np.eye(2), np.kron(i.matrix(), np.eye(2))))
-            elif i.wires.tolist() == [2]:
-                mats.append(np.kron(np.eye(4), i.matrix()))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
-                mats.append(np.kron(i.matrix(), np.eye(2)))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
-                mats.append(
-                    np.array(
-                        [
-                            [1, 0, 0, 0, 0, 0, 0, 0],
-                            [0, 1, 0, 0, 0, 0, 0, 0],
-                            [0, 0, 1, 0, 0, 0, 0, 0],
-                            [0, 0, 0, 1, 0, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 1, 0, 0],
-                            [0, 0, 0, 0, 1, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 0, 0, 1],
-                            [0, 0, 0, 0, 0, 0, 1, 0],
-                        ]
-                    )
-                )
-
-        decomposed_matrix = np.linalg.multi_dot(mats)
-        lam = np.exp(1j * phi)
-        exp = np.array(
-            [
-                [1, 0, 0, 0, 0, 0, 0, 0],
-                [0, 1, 0, 0, 0, 0, 0, 0],
-                [0, 0, 1, 0, 0, 0, 0, 0],
-                [0, 0, 0, 1, 0, 0, 0, 0],
-                [0, 0, 0, 0, 1, 0, 0, 0],
-                [0, 0, 0, 0, 0, lam, 0, 0],
-                [0, 0, 0, 0, 0, 0, 1, 0],
-                [0, 0, 0, 0, 0, 0, 0, lam],
-            ]
-        )
-
-        assert np.allclose(decomposed_matrix, exp)
-
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_decomp_broadcasted(self, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation
-        calculates the correct decomposition"""
-        phi = np.array([-0.2, 4.2, 1.8])
-        op = cphase_op(phi, wires=[0, 2])
-        decomp = op.decomposition()
-
-        mats = []
-        for i in reversed(decomp):
-            mat = i.matrix()
-            eye = np.eye(2)[np.newaxis] if np.ndim(mat) == 3 else np.eye(2)
-            if i.wires.tolist() == [0]:
-                mats.append(np.kron(mat, np.kron(eye, eye)))
-            elif i.wires.tolist() == [1]:
-                mats.append(np.kron(eye, np.kron(mat, eye)))
-            elif i.wires.tolist() == [2]:
-                mats.append(np.kron(np.kron(eye, eye), mat))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 1]:
-                mats.append(np.kron(mat, eye))
-            elif isinstance(i, qml.CNOT) and i.wires.tolist() == [0, 2]:
-                mats.append(
-                    np.array(
-                        [
-                            [1, 0, 0, 0, 0, 0, 0, 0],
-                            [0, 1, 0, 0, 0, 0, 0, 0],
-                            [0, 0, 1, 0, 0, 0, 0, 0],
-                            [0, 0, 0, 1, 0, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 1, 0, 0],
-                            [0, 0, 0, 0, 1, 0, 0, 0],
-                            [0, 0, 0, 0, 0, 0, 0, 1],
-                            [0, 0, 0, 0, 0, 0, 1, 0],
-                        ]
-                    )
-                )
-
-        decomposed_matrix = multi_dot_broadcasted(mats)
-        lam = np.exp(1j * phi)
-        exp = np.array([np.diag([1, 1, 1, 1, 1, el, 1, el]) for el in lam])
-
-        assert np.allclose(decomposed_matrix, exp)
-
     @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
     @pytest.mark.parametrize(
         "cphase_op,lam_pos",
@@ -1667,222 +1428,6 @@ class TestMatrix:
             qml.Rot(a, b, c, wires=0).matrix(), arbitrary_rotation(a, b, c), atol=tol, rtol=0
         )
 
-    def test_CRX(self, tol):
-        """Test controlled x rotation is correct"""
-
-        # test identity for theta=0
-        assert np.allclose(qml.CRX.compute_matrix(0), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRX(0, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array(
-            [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 0, 1 / np.sqrt(2), -1j / np.sqrt(2)],
-                [0, 0, -1j / np.sqrt(2), 1 / np.sqrt(2)],
-            ]
-        )
-        assert np.allclose(qml.CRX.compute_matrix(np.pi / 2), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRX(np.pi / 2, wires=[0, 1]).matrix(), expected_pi_half, atol=tol, rtol=0
-        )
-
-        # test identity for theta=pi
-        expected_pi = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1j], [0, 0, -1j, 0]])
-        assert np.allclose(qml.CRX.compute_matrix(np.pi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRX(np.pi, wires=[0, 1]).matrix(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        expected = [expected_pi_half, expected_pi]
-        assert np.allclose(qml.CRX.compute_matrix(param), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRX(param, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0)
-
-    def test_CRY(self, tol):
-        """Test controlled y rotation is correct"""
-
-        # test identity for theta=0
-        assert np.allclose(qml.CRY.compute_matrix(0), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRY(0, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array(
-            [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 0, 1 / np.sqrt(2), -1 / np.sqrt(2)],
-                [0, 0, 1 / np.sqrt(2), 1 / np.sqrt(2)],
-            ]
-        )
-        assert np.allclose(qml.CRY.compute_matrix(np.pi / 2), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRY(np.pi / 2, wires=[0, 1]).matrix(), expected_pi_half, atol=tol, rtol=0
-        )
-
-        # test identity for theta=pi
-        expected_pi = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]])
-        assert np.allclose(qml.CRY.compute_matrix(np.pi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRY(np.pi, wires=[0, 1]).matrix(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        expected = [expected_pi_half, expected_pi]
-        assert np.allclose(qml.CRY.compute_matrix(param), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRY(param, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0)
-
-    def test_CRZ(self, tol):
-        """Test controlled z rotation is correct"""
-
-        # test identity for theta=0
-        assert np.allclose(qml.CRZ.compute_matrix(0), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(0, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array(
-            [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 0, np.exp(-1j * np.pi / 4), 0],
-                [0, 0, 0, np.exp(1j * np.pi / 4)],
-            ]
-        )
-        assert np.allclose(qml.CRZ.compute_matrix(np.pi / 2), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRZ(np.pi / 2, wires=[0, 1]).matrix(), expected_pi_half, atol=tol, rtol=0
-        )
-
-        # test identity for theta=pi
-        expected_pi = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1j, 0], [0, 0, 0, 1j]])
-        assert np.allclose(qml.CRZ.compute_matrix(np.pi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(np.pi, wires=[0, 1]).matrix(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        expected = [expected_pi_half, expected_pi]
-        assert np.allclose(qml.CRZ.compute_matrix(param), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(param, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0)
-
-    @pytest.mark.tf
-    def test_CRZ_tf(self, tol):
-        """Test controlled z rotation is correct when used with Tensorflow,
-        because the code differs in that case."""
-        import tensorflow as tf
-
-        # test identity for theta=0
-        z = tf.Variable(0.0)
-        assert np.allclose(qml.CRZ.compute_matrix(z), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(z, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array(
-            [
-                [1, 0, 0, 0],
-                [0, 1, 0, 0],
-                [0, 0, np.exp(-1j * np.pi / 4), 0],
-                [0, 0, 0, np.exp(1j * np.pi / 4)],
-            ]
-        )
-        phi = tf.Variable(np.pi / 2)
-        assert np.allclose(qml.CRZ.compute_matrix(phi), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(phi, wires=[0, 1]).matrix(), expected_pi_half, atol=tol, rtol=0)
-
-        # test identity for theta=pi
-        expected_pi = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, -1j, 0], [0, 0, 0, 1j]])
-        phi = tf.Variable(np.pi)
-        assert np.allclose(qml.CRZ.compute_matrix(phi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(phi, wires=[0, 1]).matrix(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        expected = [expected_pi_half, expected_pi]
-        param_tf = tf.Variable(param)
-        assert np.allclose(qml.CRZ.compute_matrix(param_tf), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(param_tf, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0)
-
-    def test_CRot(self, tol):
-        """Test controlled arbitrary rotation is correct"""
-
-        # test identity for phi,theta,omega=0
-        assert np.allclose(qml.CRot.compute_matrix(0, 0, 0), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRot(0, 0, 0, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0
-        )
-
-        # test identity for phi,theta,omega=pi
-        expected = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]])
-        assert np.allclose(qml.CRot.compute_matrix(np.pi, np.pi, np.pi), expected, atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRot(np.pi, np.pi, np.pi, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0
-        )
-
-        def arbitrary_Crotation(x, y, z):
-            """controlled arbitrary single qubit rotation"""
-            c = np.cos(y / 2)
-            s = np.sin(y / 2)
-            return np.array(
-                [
-                    [1, 0, 0, 0],
-                    [0, 1, 0, 0],
-                    [0, 0, np.exp(-0.5j * (x + z)) * c, -np.exp(0.5j * (x - z)) * s],
-                    [0, 0, np.exp(-0.5j * (x - z)) * s, np.exp(0.5j * (x + z)) * c],
-                ]
-            )
-
-        a, b, c = 0.432, -0.152, 0.9234
-        assert np.allclose(
-            qml.CRot.compute_matrix(a, b, c), arbitrary_Crotation(a, b, c), atol=tol, rtol=0
-        )
-        assert np.allclose(
-            qml.CRot(a, b, c, wires=[0, 1]).matrix(),
-            arbitrary_Crotation(a, b, c),
-            atol=tol,
-            rtol=0,
-        )
-
-    def test_CRot_broadcasted(self, tol):
-        """Test broadcasted controlled arbitrary rotation is correct"""
-
-        # test identity for phi,theta,omega=0
-        z = np.zeros(5)
-        assert np.allclose(qml.CRot.compute_matrix(z, z, z), np.identity(4), atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRot(z, z, z, wires=[0, 1]).matrix(), np.identity(4), atol=tol, rtol=0
-        )
-
-        # test -i*CY for phi,theta,omega=pi
-        expected = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, -1], [0, 0, 1, 0]])
-        pi = np.ones(3) * np.pi
-        assert np.allclose(qml.CRot.compute_matrix(pi, pi, pi), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRot(pi, pi, pi, wires=[0, 1]).matrix(), expected, atol=tol, rtol=0)
-
-        def arbitrary_Crotation(x, y, z):
-            """controlled arbitrary single qubit rotation"""
-            c = np.cos(y / 2)
-            s = np.sin(y / 2)
-            return np.array(
-                [
-                    [
-                        [1, 0, 0, 0],
-                        [0, 1, 0, 0],
-                        [0, 0, np.exp(-0.5j * (_x + _z)) * _c, -np.exp(0.5j * (_x - _z)) * _s],
-                        [0, 0, np.exp(-0.5j * (_x - _z)) * _s, np.exp(0.5j * (_x + _z)) * _c],
-                    ]
-                    for _x, _z, _c, _s in zip(x, z, c, s)
-                ]
-            )
-
-        a, b, c = np.array([0.432, -0.124]), np.array([-0.152, 2.912]), np.array([0.9234, -9.2])
-        assert np.allclose(
-            qml.CRot.compute_matrix(a, b, c), arbitrary_Crotation(a, b, c), atol=tol, rtol=0
-        )
-        assert np.allclose(
-            qml.CRot(a, b, c, wires=[0, 1]).matrix(),
-            arbitrary_Crotation(a, b, c),
-            atol=tol,
-            rtol=0,
-        )
-
     def test_U2_gate(self, tol):
         """Test U2 gate matrix matches the documentation"""
         phi = 0.432
@@ -1961,73 +1506,6 @@ class TestMatrix:
         assert np.allclose(qml.U3.compute_matrix(theta, phi, lam), expected, atol=tol, rtol=0)
         assert np.allclose(qml.U3(theta, phi, lam, wires=[0]).matrix(), expected, atol=tol, rtol=0)
 
-    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_matrix_and_eigvals(self, phi, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation calculate the correct
-        matrix and eigenvalues"""
-        op = cphase_op(phi, wires=[0, 1])
-        res = op.matrix()
-        exp = ControlledPhaseShift(phi)
-        assert np.allclose(res, exp)
-
-        res = op.eigvals()
-        assert np.allclose(np.diag(res), exp)
-
-    @pytest.mark.tf
-    @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_matrix_and_eigvals_tf(self, phi, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation calculate the correct
-        matrix and eigenvalues for the Tensorflow interface, because the code differs
-        in that case."""
-        import tensorflow as tf
-
-        op = cphase_op(tf.Variable(phi), wires=[0, 1])
-        res = op.matrix()
-        exp = ControlledPhaseShift(phi)
-        assert np.allclose(res, exp)
-
-        res = op.eigvals()
-        assert np.allclose(np.diag(res), exp)
-
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_matrix_and_eigvals_broadcasted(self, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation calculate the
-        correct matrix and eigenvalues for broadcasted parameters"""
-        phi = np.array([0.2, np.pi / 2, -0.1])
-        op = cphase_op(phi, wires=[0, 1])
-        res = op.matrix()
-        expected = np.array([np.eye(4, dtype=complex)] * 3)
-        expected[..., 3, 3] = np.exp(1j * phi)
-        assert np.allclose(res, expected)
-
-        res = op.eigvals()
-        exp_eigvals = np.ones((3, 4), dtype=complex)
-        exp_eigvals[..., 3] = np.exp(1j * phi)
-        assert np.allclose(res, exp_eigvals)
-
-    @pytest.mark.tf
-    @pytest.mark.parametrize("cphase_op", [qml.ControlledPhaseShift, qml.CPhase])
-    def test_controlled_phase_shift_matrix_and_eigvals_broadcasted_tf(self, cphase_op):
-        """Tests that the ControlledPhaseShift and CPhase operation calculate the
-        correct matrix and eigenvalues for broadcasted parameters and Tensorflow,
-        because the code differs for that interface."""
-        import tensorflow as tf
-
-        phi = np.array([0.2, np.pi / 2, -0.1])
-        phi_tf = tf.Variable(phi)
-        op = cphase_op(phi_tf, wires=[0, 1])
-        res = op.matrix()
-        expected = np.array([np.eye(4, dtype=complex)] * 3)
-        expected[..., 3, 3] = np.exp(1j * phi)
-        assert np.allclose(res, expected)
-
-        res = op.eigvals()
-        exp_eigvals = np.ones((3, 4), dtype=complex)
-        exp_eigvals[..., 3] = np.exp(1j * phi)
-        assert np.allclose(res, exp_eigvals)
-
     @pytest.mark.parametrize("phi", [-0.1, 0.2, 0.5])
     @pytest.mark.parametrize(
         "cphase_op,gate_data_mat",
@@ -2207,61 +1685,6 @@ class TestEigvals:
         )
         assert np.allclose(op.eigvals(), expected)
 
-    def test_crz_eigvals(self, tol):
-        """Test controlled z rotation eigvals are correct"""
-
-        # test identity for theta=0
-        assert np.allclose(qml.CRZ.compute_eigvals(0), np.ones(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(0, wires=[0, 1]).eigvals(), np.ones(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array([1, 1, np.exp(-1j * np.pi / 4), np.exp(1j * np.pi / 4)])
-        assert np.allclose(qml.CRZ.compute_eigvals(np.pi / 2), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(
-            qml.CRZ(np.pi / 2, wires=[0, 1]).eigvals(), expected_pi_half, atol=tol, rtol=0
-        )
-
-        # test identity for theta=pi
-        expected_pi = np.array([1, 1, -1j, 1j])
-        assert np.allclose(qml.CRZ.compute_eigvals(np.pi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(np.pi, wires=[0, 1]).eigvals(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        expected = [expected_pi_half, expected_pi]
-        assert np.allclose(qml.CRZ.compute_eigvals(param), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(param, wires=[0, 1]).eigvals(), expected, atol=tol, rtol=0)
-
-    @pytest.mark.tf
-    def test_crz_eigvals_tf(self, tol):
-        """Test controlled z rotation eigvals are correct with Tensorflow, because the
-        code differs for that interface."""
-        import tensorflow as tf
-
-        # test identity for theta=0
-        z = tf.Variable(0)
-        assert np.allclose(qml.CRZ.compute_eigvals(z), np.ones(4), atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(z, wires=[0, 1]).eigvals(), np.ones(4), atol=tol, rtol=0)
-
-        # test identity for theta=pi/2
-        expected_pi_half = np.array([1, 1, np.exp(-1j * np.pi / 4), np.exp(1j * np.pi / 4)])
-        phi = tf.Variable(np.pi / 2)
-        assert np.allclose(qml.CRZ.compute_eigvals(phi), expected_pi_half, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(phi, wires=[0, 1]).eigvals(), expected_pi_half, atol=tol, rtol=0)
-
-        # test identity for theta=pi
-        expected_pi = np.array([1, 1, -1j, 1j])
-        phi = tf.Variable(np.pi)
-        assert np.allclose(qml.CRZ.compute_eigvals(phi), expected_pi, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(phi, wires=[0, 1]).eigvals(), expected_pi, atol=tol, rtol=0)
-
-        # test broadcasting
-        param = np.array([np.pi / 2, np.pi])
-        param_tf = tf.Variable(param)
-        expected = [expected_pi_half, expected_pi]
-        assert np.allclose(qml.CRZ.compute_eigvals(param_tf), expected, atol=tol, rtol=0)
-        assert np.allclose(qml.CRZ(param_tf, wires=[0, 1]).eigvals(), expected, atol=tol, rtol=0)
-
     def test_global_phase_eigvals(self):
         """Test GlobalPhase eigenvalues are correct"""
 
@@ -4090,7 +3513,7 @@ class TestSimplify:
         simplified_op = qml.simplify(unsimplified_op)
 
         if op != qml.PSWAP:
-            assert qml.equal(simplified_op, qml.Identity(0))
+            assert isinstance(simplified_op, qml.Identity)
         else:
             # PSWAP reduces to SWAP when the angle is 0
             assert qml.equal(simplified_op, qml.SWAP(wires=[0, 1]))
@@ -4131,36 +3554,6 @@ class TestSimplify:
         assert not_simplified_rot.name == "Rot"
         assert np.allclose(not_simplified_rot.matrix(), rot.matrix())
 
-    def test_simplify_crot(self):
-        """Simplify CRot operations with different parameters."""
-
-        crot_x = qml.CRot(np.pi / 2, 0.1, -np.pi / 2, wires=[0, 1])
-        simplify_crot_x = crot_x.simplify()
-
-        assert simplify_crot_x.name == "CRX"
-        assert simplify_crot_x.data == (0.1,)
-        assert np.allclose(simplify_crot_x.matrix(), crot_x.matrix())
-
-        crot_y = qml.CRot(0, 0.1, 0, wires=[0, 1])
-        simplify_crot_y = crot_y.simplify()
-
-        assert simplify_crot_y.name == "CRY"
-        assert simplify_crot_y.data == (0.1,)
-        assert np.allclose(simplify_crot_y.matrix(), crot_y.matrix())
-
-        crot_z = qml.CRot(0.1, 0, 0.2, wires=[0, 1])
-        simplify_crot_z = crot_z.simplify()
-
-        assert simplify_crot_z.name == "CRZ"
-        assert np.allclose(simplify_crot_z.data, [0.3])
-        assert np.allclose(simplify_crot_z.matrix(), crot_z.matrix())
-
-        crot = qml.CRot(0.1, 0.2, 0.3, wires=[0, 1])
-        not_simplified_crot = crot.simplify()
-
-        assert not_simplified_crot.name == "CRot"
-        assert np.allclose(not_simplified_crot.matrix(), crot.matrix())
-
     def test_simplify_u2(self):
         """Simplify u2 operations with different parameters."""
 
@@ -4217,22 +3610,6 @@ class TestSimplify:
         assert np.allclose(u3_not_simplified.matrix(), u3.matrix())
 
 
-controlled_data = [
-    (qml.RX(1.234, wires=0), qml.CRX(1.234, wires=("a", 0))),
-    (qml.RY(1.234, wires=0), qml.CRY(1.234, wires=("a", 0))),
-    (qml.RZ(1.234, wires=0), qml.CRZ(1.234, wires=("a", 0))),
-    (qml.PhaseShift(1.234, wires=0), qml.ControlledPhaseShift(1.234, wires=("a", 0))),
-    (qml.Rot(1.2, 2.3, 3.4, wires=0), qml.CRot(1.2, 2.3, 3.4, wires=("a", 0))),
-]
-
-
-@pytest.mark.parametrize("base, cbase", controlled_data)
-def test_controlled_method(base, cbase):
-    """Tests the _controlled method for parametric ops."""
-    # pylint: disable=protected-access
-    assert qml.equal(base._controlled("a"), cbase)
-
-
 label_data = [
     (
         qml.Rot(1.23456, 2.3456, 3.45678, wires=0),
@@ -4262,12 +3639,6 @@ label_data = [
         "∏_ϕ\n(1.23)",
         "∏_ϕ\n(1)",
     ),
-    (
-        qml.ControlledPhaseShift(1.2345, wires=(0, 1)),
-        "Rϕ",
-        "Rϕ\n(1.23)",
-        "Rϕ\n(1)",
-    ),
     (
         qml.CPhaseShift00(1.2345, wires=(0, 1)),
         "Rϕ(00)",
@@ -4286,15 +3657,6 @@ label_data = [
         "Rϕ(10)\n(1.23)",
         "Rϕ(10)\n(1)",
     ),
-    (qml.CRX(1.234, wires=(0, 1)), "RX", "RX\n(1.23)", "RX\n(1)"),
-    (qml.CRY(1.234, wires=(0, 1)), "RY", "RY\n(1.23)", "RY\n(1)"),
-    (qml.CRZ(1.234, wires=(0, 1)), "RZ", "RZ\n(1.23)", "RZ\n(1)"),
-    (
-        qml.CRot(1.234, 2.3456, 3.456, wires=(0, 1)),
-        "Rot",
-        "Rot\n(1.23,\n2.35,\n3.46)",
-        "Rot\n(1,\n2,\n3)",
-    ),
     (qml.U1(1.2345, wires=0), "U1", "U1\n(1.23)", "U1\n(1)"),
     (qml.U2(1.2345, 2.3456, wires=0), "U2", "U2\n(1.23,\n2.35)", "U2\n(1,\n2)"),
     (
@@ -4321,6 +3683,22 @@ label_data = [
         "IsingZZ\n(1.23)",
         "IsingZZ\n(1)",
     ),
+    # Controlled operations
+    (qml.CRX(1.234, wires=(0, 1)), "RX", "RX\n(1.23)", "RX\n(1)"),
+    (qml.CRY(1.234, wires=(0, 1)), "RY", "RY\n(1.23)", "RY\n(1)"),
+    (qml.CRZ(1.234, wires=(0, 1)), "RZ", "RZ\n(1.23)", "RZ\n(1)"),
+    (
+        qml.CRot(1.234, 2.3456, 3.456, wires=(0, 1)),
+        "Rot",
+        "Rot\n(1.23,\n2.35,\n3.46)",
+        "Rot\n(1,\n2,\n3)",
+    ),
+    (
+        qml.ControlledPhaseShift(1.2345, wires=(0, 1)),
+        "Rϕ",
+        "Rϕ\n(1.23)",
+        "Rϕ\n(1)",
+    ),
 ]
 
 # labels with broadcasted parameters are not implemented properly yet, the parameters are truncated
@@ -4533,7 +3911,7 @@ control_data = [
     (qml.IsingYY(np.array([-5.1, 0.219]), wires=(0, 1)), Wires([])),
     (qml.IsingZZ(1.234, wires=(0, 1)), Wires([])),
     (qml.PSWAP(1.234, wires=(0, 1)), Wires([])),
-    ### Controlled Ops
+    # Controlled Ops
     (qml.ControlledPhaseShift(1.234, wires=(0, 1)), Wires(0)),
     (qml.CPhaseShift00(1.234, wires=(0, 1)), Wires(0)),
     (qml.CPhaseShift01(1.234, wires=(0, 1)), Wires(0)),
@@ -4566,8 +3944,5 @@ def test_control_values(op, control_values):
 
 def test_op_aliases_are_valid():
     """Tests that ops in new files can still be accessed from the old parametric_ops module."""
-    assert (
-        qml.ops.qubit.parametric_ops_controlled.ControlledPhaseShift is old_loc_ControlledPhaseShift
-    )
     assert qml.ops.qubit.parametric_ops_multi_qubit.MultiRZ is old_loc_MultiRZ
     assert qml.ops.qubit.parametric_ops_single_qubit.RX is old_loc_RX
diff --git a/tests/templates/test_subroutines/test_controlled_sequence.py b/tests/templates/test_subroutines/test_controlled_sequence.py
index 1f556fde9..9c937fc8b 100644
--- a/tests/templates/test_subroutines/test_controlled_sequence.py
+++ b/tests/templates/test_subroutines/test_controlled_sequence.py
@@ -165,9 +165,9 @@ class TestMethods:
 
         decomp = op.compute_decomposition(base=op.base, control_wires=op.control, lazy=False)
         expected_decomp = [
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 4, wires=3), control_wires="a"),
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 2, wires=3), control_wires=1),
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 1, wires=3), control_wires="blue"),
+            qml.CRX(0.25 * 4, wires=["a", 3]),
+            qml.CRX(0.25 * 2, wires=[1, 3]),
+            qml.CRX(0.25 * 1, wires=["blue", 3]),
         ]
 
         for op1, op2 in zip(decomp, expected_decomp):
@@ -178,9 +178,9 @@ class TestMethods:
 
         decomp = op.decomposition()
         expected_decomp = [
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 4, wires=3), control_wires="a"),
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 2, wires=3), control_wires=1),
-            qml.ops.op_math.Controlled(qml.RX(0.25 * 1, wires=3), control_wires="blue"),
+            qml.CRX(0.25 * 4, wires=["a", 3]),
+            qml.CRX(0.25 * 2, wires=[1, 3]),
+            qml.CRX(0.25 * 1, wires=["blue", 3]),
         ]
 
         for op1, op2 in zip(decomp, expected_decomp):
diff --git a/tests/transforms/test_defer_measurements.py b/tests/transforms/test_defer_measurements.py
index 5b533e3cf..b8104f0ff 100644
--- a/tests/transforms/test_defer_measurements.py
+++ b/tests/transforms/test_defer_measurements.py
@@ -1372,7 +1372,7 @@ class TestQubitReuseAndReset:
             qml.Hadamard(0),
             qml.CNOT([0, 2]),
             qml.CNOT([2, 0]),
-            qml.ops.Controlled(qml.RX(0.123, 1), 2),
+            qml.CRX(0.123, wires=[2, 1]),
             qml.expval(qml.PauliZ(1)),
         ]