In vectorial problems, we end up with linear systems of higher order, such as

for $i=1,\dots,N$ and $j=1,\dots,M$.

These systems cannot be solved readily with existing software. In order to be able to solve them with existing software, we need to reshape them by defining a matrix of matrices $\BAhat$ and vector of vectors $\Buhat$ and $\Bbhat$:

This allows us to express the linear system as

Here, we reshape the system by defining a map $i_d$ that maps original indices to the reshaped indices

where we used 1-based indexing of the arrays. We set

and write

For reference, the inverse of the index mapping reads

Thus, we have for our reshaped indices,

Expressed as a regular linear system \eqref{eq:3}, the higher-order system \eqref{eq:1} can be solved with a linear solver such as LAPACK.