The matrix is composed of submatrices of the form $\begin{bmatrix}u & -v \\ v & u\end{bmatrix}$, for example $\begin{bmatrix}a & -c \\ c & a\end{bmatrix}$ is in the first and the third rows and columns. These matrices all commute and can be simultaneously diagonalized: we have $\begin{bmatrix}u & -v \\ v & u\end{bmatrix}\begin{bmatrix}1\\i\end{bmatrix} = (u - iv) \begin{bmatrix}1\\i\end{bmatrix}$, $\begin{bmatrix}u & -v \\ v & u\end{bmatrix}\begin{bmatrix}1\\-i\end{bmatrix} = (u + iv) \begin{bmatrix}1\\-i\end{bmatrix}$, so after the basis change to $\begin{bmatrix}1\\i\end{bmatrix}, \begin{bmatrix}1\\-i\end{bmatrix}$ the matrix is $\begin{bmatrix}u - iv & 0\\0 & u + iv\end{bmatrix}$. For computational reasons it's better to use the basis $\frac12\begin{bmatrix}1\\i\end{bmatrix}, \frac12\begin{bmatrix}1\\-i\end{bmatrix}$, so that a vector $\begin{bmatrix}x\\y\end{bmatrix}$ is transformed to $\begin{bmatrix}x'\\y'\end{bmatrix} = \begin{bmatrix}x -iy \\ x + iy\end{bmatrix}$. The diagonalized matrix is of course still the same.
In our matrix we have $4$ submatrices like this, so the equation $$\begin{bmatrix}a & -b & -c & d \\ b & e & -d & -f \\ c & -d & a & -b \\ d & f & b & e\end{bmatrix}\begin{bmatrix}x\\y\\z\\t\end{bmatrix} = \begin{bmatrix}p\\q\\r\\s\end{bmatrix}$$ after the basis change to $\frac12\begin{bmatrix}1\\0\\i\\0\end{bmatrix}, \frac12\begin{bmatrix}0\\1\\0\\i\end{bmatrix}, \frac12\begin{bmatrix}1\\0\\-i\\0\end{bmatrix}, \frac12\begin{bmatrix}0\\1\\0\\-i\end{bmatrix}$ translates to a block diagonal system $$\begin{bmatrix}a - ic & -(b - id) & 0 & 0 \\ b - id & e - if & 0 & 0 \\ 0 & 0 & a + ic & -(b + id) \\ 0 & 0 & b + id & e + if\end{bmatrix}\begin{bmatrix}x - iz \\ y - it \\ x + iz \\ y + it\end{bmatrix} = \begin{bmatrix}p - ir \\ q - is \\ p + ir \\ q + is\end{bmatrix},$$ we only need to solve two $2 \times 2$ systems and then find $x,y,z,t$, for example, $x = \frac{(x + iz) + (x - iz)}{2}$.
