## 3 Parallel Solution and Preconditioning of Triangular Toeplitz Linear Systems
### Preconditioning for (Block) Triangular Toeplitz Systems  
The algorithm described in the previous subsection can be applied in a very straightforward way to the preconditioning of linear systems with (block) triangular Toeplitz structure, allowing for an extremely efficient computation of the solution. On the other hand, special attention must be paid to the issue of stability, which strongly relies on the parameter \(\omega\). We consider an example to present the approach, but we emphasize that the same basic ideas pertain to a variety of situations. The example we use is taken from Ref. [21]. Here, [PERSON] and [PERSON] seek the numerical solution of the linear wave equation model  
\[ \begin{cases}y_{tt}-\Delta y=f,&\text{in }\Omega\times(0,T),\\ y=0,&\text{on }\partial\Omega\times(0,T),\\ y(\cdot,0)=\psi_{0},\ y_{t}(\cdot,0)=\psi_{1}&\text{in }\Omega,\end{cases}  \tag{4}\]  
where \(\Omega\subset\mathbb{R}^{d}\) with \(d>1\) is a bounded and open domain with Lipschitz boundary, \(\psi_{0}\), \(\psi_{1}\) are the initial conditions and \(f\) is a given source term. The equations are discretized all-at-once in time by means of implicit finite-difference schemes (see Ref. [21] for details). In the two-dimensional case with \(\Omega\) a rectangle, the attained matrix shows the following block lower triangular Toeplitz structure  
\[\mathcal{K}=\frac{1}{\tau^{2}} \begin{bmatrix}L&&&&\\ -2I_{m}&L&&&\\ L&-2I_{m}&L&&\\ &\ddots&\ddots&\ddots&\\ &&&L&-2I_{m}&L\end{bmatrix} ,\qquad L=I_{m}-\frac{\tau^{2}}{2}\Delta_{ht},\]  
in which  
* \(m\) is the number of grid nodes in the spatial mesh and \(h=\frac{1}{m+1}\) is the spatial step size;
* \(\tau\) is the time mesh step size, obtained as \(\frac{T}{n}\), where \(n\) is the number of time steps;
* \(I_{m}\) is the identity matrix of size \(m\);
* \(\Delta_{h}\) is the matrix obtained by discretizing the Laplacian operator \(\Delta\) in (4) with the central finite difference method. The associated linear system encompasses all the time steps at once and its solution corresponds to the solution to (4) at each time step simultaneously. In other words, if \(y_{j}\in\mathbb{R}^{m}\) denotes the solution to (4) at the \(j\)-th time step, the system has the form  
\[\mathcal{K}\mathbf{y}=\mathbf{b},\quad\mathbf{y}= \begin{bmatrix}y_{1}\\ \vdots\\ y_{n}\end{bmatrix} ,\ \mathbf{b}= \begin{bmatrix}b_{1}\\ \vdots\\ b_{n}\end{bmatrix} ,\ b_{j}\in\mathbb{R}^{m},\ j=1,\ldots,n\]and therefore the vectors \(y_{j}\), \(j=1,\ldots,n\) are computed in parallel as the solution to the all-at-once linear system with the coefficient matrix \(\mathcal{K}\). To solve this linear system, the authors adopt a preconditioned GMRES method and construct a class of block \(\omega\)-circulant preconditioners as follows. Defining \(B_{1}\) and \(B_{2}\) as  
\[B_{1}= \begin{bmatrix}1\\ 0&1\\ 1&0&1\\ &\ddots&\ddots&\ddots\\ &&1&0&1\end{bmatrix} _{n\times n},\qquad B_{2}= \begin{bmatrix}0\\ 1&0\\ &1&0\\ &&\ddots&\ddots\\ &&1&0\end{bmatrix} _{n\times n},\]  
as described in Ref. [22], \(\mathcal{K}\) can be expressed in the form  
\[\mathcal{K}=\frac{1}{\tau^{2}}(B_{1}\otimes L-B_{2}\otimes 2I_{m}).\]  
Then, given the matrix \(\Pi_{n,\omega}\) as in Section 2 and  
\[\Psi_{n,\omega}= \begin{bmatrix}1&&\omega&0\\ 0&1&&\omega\\ 1&0&1&\\ &\ddots&\ddots&\ddots\\ &&1&0&1\end{bmatrix}