$ \newcommand{\BAhat}{\hat{\boldsymbol{A}}} \newcommand{\Buhat}{\hat{\boldsymbol{u}}} \newcommand{\Bbhat}{\hat{\boldsymbol{b}}} $

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

∑k=1N∑l=1MAijkl ukl=bij\begin{equation} \suml{k=1}{N} \suml{l=1}{M} A_{ijkl} \, u_{kl} = b_{ij} \htmlId{eq:1}{} \tag{1}\end{equation}

for i=1,…,Ni=1,\dots,N and j=1,…,Mj=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 A^\BAhat and vector of vectors u^\Buhat and b^\Bbhat:

[A1111⋯A111MA11N1⋯A11NM⋮⋱⋮⋯⋮⋱⋮A1M11⋯A1M1MA1MN1⋯A1MNM⋮⋱⋮AN111⋯AN11MAN1N1⋯AN1NM⋮⋱⋮⋯⋮⋱⋮ANM11⋯ANM1MANMN1⋯ANMNM]⏟A^[u11⋮u1M⋮uN1⋮uNM]⏟u^=[b11⋮b1M⋮bN1⋮bNM]⏟b^\begin{equation} \underbrace{ \left[ \begin{array}{ccc|c|ccc} A_{1111} & \cdots & A_{111M} & & A_{11N1} & \cdots & A_{11NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{1M11} & \cdots & A_{1M1M} & & A_{1MN1} & \cdots & A_{1MNM} \\[1ex] \hline & \vdots & &\ddots & & \vdots & \\ \hline &&&&&& \\[-1.5ex] A_{N111} & \cdots & A_{N11M} & & A_{N1N1} & \cdots & A_{N1NM} \\ \vdots & \ddots & \vdots & \cdots & \vdots & \ddots & \vdots \\ A_{NM11} & \cdots & A_{NM1M} & & A_{NMN1} & \cdots & A_{NMNM} \end{array} \right] }_{\BAhat} \underbrace{ \left[ \begin{array}{c} u_{11}\\ \vdots \\ u_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] u_{N1} \\ \vdots \\ u_{NM} \end{array} \right] }_{\Buhat} = \underbrace{ \left[ \begin{array}{c} b_{11}\\ \vdots \\ b_{1M} \\[1ex] \hline \vdots\\ \hline \\[-1.5ex] b_{N1} \\ \vdots \\ b_{NM} \end{array} \right]}_{\Bbhat} \tag{2}\end{equation}

This allows us to express the linear system as

A^u^=b^\begin{equation} \BAhat \Buhat = \Bbhat \htmlId{eq:3}{} \tag{3}\end{equation}

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

id:={[1,N]×[1,M]→[1,NM](i,j)↦M(i−1)+j\begin{equation} i_d := \left\{ \begin{array}{rl} [1,N]\times[1,M] & \to [1,NM]\\[1ex] (i,j) & \mapsto M(i-1)+j\\ \end{array} \right. \tag{4}\end{equation}

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

α:=id(i,j)=M(i−1)+jβ:=id(k,l)=M(k−1)+l\begin{equation} \boxed{ \begin{alignedat}{3} \alpha &:= i_d(i,j) &&= M(i-1) + j \\ \beta &:= i_d(k,l) &&= M(k-1) + l \\ \end{alignedat} } \tag{5}\end{equation}

and write

A^αβ=Aijkl,u^β=uklandb^α=bij\begin{equation} \hat{A}_{\alpha\beta} = A_{ijkl} \quad,\quad \hat{u}_{\beta} = u_{kl} \eqand \hat{b}_{\alpha} = b_{ij} \tag{6}\end{equation}

For reference, the inverse of the index mapping reads

id−1:={[1,NM]→[1,N]×[1,M]α↦(1+(α−mod(α,M))/M  ,  mod(α,M))\begin{equation} i_d^{-1} := \left\{ \begin{array}{rl} [1,NM] & \to [1,N]\times[1,M] \\[1ex] \alpha & \mapsto (1+(\alpha-\modop(\alpha,M))/M\;,\; \modop(\alpha,M)) \end{array} \right. \tag{7}\end{equation}

Thus, we have for our reshaped indices,

j=mod(α,M)i=1+(α−j)/Mandl=mod(β,M)k=1+(β−l)/M\begin{equation} \begin{aligned} j &= \modop(\alpha,M) \\ i &= 1+(\alpha-j)/M \end{aligned} \quad\eqand\quad \begin{aligned} l &= \modop(\beta,M) \\ k &= 1+(\beta-l)/M \end{aligned} \tag{8}\end{equation}

Expressed as a regular linear system (3), the higher-order system (1) can be solved with a linear solver such as LAPACK.