The method of direct elimination can be used to solve the constrained least squares problem \begin{equation} \min_{\mathbf{x}}\left\Vert \mathbf{Ax}-\mathbf{b}\right\Vert _{2} \end{equation} \begin{equation} \textrm{such that } \mathbf{Cx = d} \end{equation} where $\mathbf{A} \in \mathbb{R}^{m \times n}$, $\mathbf{C} \in \mathbb{R}^{p \times n}$ and $\mathbf{x} \in \mathbb{R}^{n}$. A solution to the above is unique iff $\operatorname{rank}(\mathbf{B}) = p$ and $\operatorname{rank}([\mathbf{A \; B}]^T)=n$.
This solution methodology (see Björck 1996 page 188) involves deriving an equivalent unconstrained least squares problem of lower dimension. To obtain the solution, we begin by computing the QR decomposition of $\mathbf{C}$, \begin{equation} \mathbf{C}\varPi_{\mathbf{C}}=\mathbf{Q}_{C}^{T}\left(\begin{array}{cc} \mathbf{R}_{11} & \mathbf{R}_{12}\\ 0 & 0 \end{array}\right), \end{equation} where $\varPi_{\mathbf{C}}$ is the permutation matrix. The number of rows in $\mathbf{R}_{11}$ is $r$. Using the permutation and orthogonal matrices, let us define \begin{equation} \mathbf{A}\varPi_{\mathbf{C}}=\left(\begin{array}{cc} \bar{\mathbf{A}}_{1} & \bar{\mathbf{A}}_{2}\end{array}\right),\; \varPi_{\mathbf{C}}^{T}\mathbf{x}=\left(\begin{array}{c} \bar{\mathbf{x}}_{1}\\ \bar{\mathbf{x}}_{2} \end{array}\right),\; \mathbf{Q}_{\mathbf{C}}^{T}\mathbf{d}=\left(\begin{array}{c} \bar{\mathbf{d}}_{1}\\ \bar{\mathbf{d}}_{2} \end{array}\right). \end{equation} Then the solution to the constrained least squares problem is given by \begin{equation} \min_{\bar{\mathbf{x}}} \left\Vert \hat{\mathbf{A}}_{2}\bar{\mathbf{x}}_{2}-\hat{\mathbf{b}}\right\Vert_{2}, \end{equation} where \begin{equation} \hat{\mathbf{A}}_{2}=\bar{\mathbf{A}}_{2}-\bar{\mathbf{A}}_{1}\mathbf{R}_{11}^{-1}\mathbf{R}_{12},\; \; \hat{\mathbf{b}}=\mathbf{b}-\bar{\mathbf{A}}_{1}\mathbf{R}_{11}^{-1}\bar{\mathbf{d}}_1. \end{equation} Now, how can we incorporate multiple linear constraints like $\mathbf{Ex=f}$--in addition to $\mathbf{Cx=d}$--using a variation of this method?
There's nothing particularly special about "multiple linear constraints". In fact, if $C$ already has more than one row, you already have multiple linear constraints. What you need to do here is simply stack the matrices on top of each other. That is: $$\begin{array}{ll} \text{minimize} & \|Ax-b\|_2 \\ \text{subject to} & C x = d \\ & E x = f \end{array}$$ is equivalent to $$\begin{array}{ll} \text{minimize} & \|Ax-b\|_2 \\ \text{subject to} & \bar{C} x = \bar{d} \end{array}$$ where we have defined $$\bar{C} \triangleq \begin{bmatrix} C \\ E \end{bmatrix} \qquad \bar{d} \triangleq \begin{bmatrix} d \\ f \end{bmatrix}$$ Now you are free to apply all of the same derivations as before.