## 3 Parallel Solution and Preconditioning of Triangular Toeplitz Linear Systems
### Parallel Solution of (Block) Triangular Toeplitz Systems  
Since a diagonalization of \(\omega\)-circulant matrices through the Fourier matrix is available, it can be used to solve a linear system of the form \(\mathbf{C}_{\omega}\mathbf{x}_{\omega}=\mathbf{r}\), with \(\mathbf{C}_{\omega}\) as in Definition 2, \(\mathbf{x}_{\omega},\mathbf{r}\in\mathbb{R}^{n}\) and \(\omega\
eq 0\), in an efficient way. More precisely, supposing that \(\mathbf{C}_{\omega}\) is invertible and deriving from Proposition 2 the spectral decomposition  
\[\mathbf{C}_{\omega}^{-1}=D_{\omega}^{-1}F_{n}^{*}L_{n,\omega}^{-1}F_{n}D_{ \omega},\]  
\(\mathbf{x}_{\omega}\) can be calculated as the matrix-vector product \(\mathbf{C}_{\omega}^{-1}\mathbf{r}\) via the following algorithm  
1. Compute \(\mathbf{y}=F_{n}D_{\omega}\mathbf{r}\) via the FFT, with sequential cost \(\mathcal{O}(n\log(n))\);
2. Compute \(\mathbf{z}=L_{n,\omega}^{-1}\mathbf{y}\), with sequential cost \(\mathcal{O}(n)\);
3. Compute \(\mathbf{x}_{\omega}=D_{\omega}^{-1}F_{n}^{*}\mathbf{z}\) via the FFT, with sequential cost \(\mathcal{O}(n\log(n))\).  
This algorithm has a total sequential cost of \(\mathcal{O}(n\log(n))\), but, since all the steps can be fully parallelized, it can be lowered to \(\mathcal{O}(\log(n))\) in a parallel model of computation.  
Now, the key idea to compute the solution of a (lower) triangular Toeplitz linear system \(A\mathbf{x}=\mathbf{r}\), with \(\mathbf{x},\mathbf{r}\in\mathbb{R}^{n}\) and  
\[A=\begin{bmatrix}a_{0}&&\\ a_{1}&a_{0}&&\\ \vdots&\ddots&\ddots&\\ a_{n-1}&\cdots&a_{1}&a_{0}\end{bmatrix},\qquad a_{0},a_{1},\ldots,a_{n-1}\in \mathbb{C},\]  
consists in completing \(A\) to the following \(\omega\)-circulant  
\[\mathbf{C}_{\omega}=\begin{bmatrix}a_{0}&\omega a_{n-1}&\cdots&\omega a_{1} \\ a_{1}&a_{0}&\ddots&\vdots\\ \vdots&\ddots&\ddots&\omega a_{n-1}\\ a_{n-1}&\cdots&a_{1}&a_{0}\end{bmatrix}=A+\omega\begin{bmatrix}0&a_{n-1}& \cdots&a_{1}\\ 0&\ddots&\vdots\\ &&\ddots&a_{n-1}\\ &&&0\end{bmatrix}\]  
and observing that  
\[\lim_{\omega\to 0}\mathbf{C}_{\omega}^{-1}=A^{-1},\quad\lim_{\omega\to 0} \mathbf{x}_{\omega}^{-1}=\mathbf{x},\quad\text{where }\mathbf{x}_{\omega}=\mathbf{C}_{\omega}^{-1}\mathbf{r}.\]Hence, applying the algorithm described above, the solution of interest \(\mathbf{x}\) can be approximated by computing \(\mathbf{x}_{\omega}\), in theory up to an arbitrarily small error for \(\omega\to 0\). Nevertheless, in practice an arbitrarily small error may not be attainable, since the condition number of the transform \(F_{n}D_{\omega}\) corresponds to the condition number of \(D_{\omega}\), which tends to \(\max\{|\omega|,|\omega|^{-1}\}\) as \(n\to\infty\) and thus tends to infinity as \(\omega\to 0\), causing stability issues. See the next subsection for a discussion in the setting of preconditioning.  
The cost of this computation amounts to \(\mathcal{O}(n\log(n))\) in a sequential implementation and to \(\mathcal{O}(\log(n))\) in a parallel implementation. We refer the reader to Ref. [8] for details and two strategies that allow for the retrieval of the exact solution.  
We stress that, although for the sake of clarity the above discussion takes into account the lower triangular case only, the technique can be easily extended to block triangular and upper triangular Toeplitz matrices.