## 3 Parallel Solution and Preconditioning of Triangular Toeplitz Linear Systems
### Preconditioning for (Block) Triangular Toeplitz Systems
,\]  
the generalized preconditioner is defined for \(\omega\in(0,1]\) as  
\[\mathcal{P}_{\omega}=\frac{1}{\tau^{2}}(\Psi_{n,\omega}\otimes L-\Pi_{n, \omega}\otimes 2I_{m}).\]  
By exploiting the simultaneous diagonalization of the \(\omega\)-circulant matrices \(\Pi_{n,\omega}\) and \(\Gamma_{n,\omega}\) and the properties of the Kronecker product, \(\mathcal{P}_{\omega}\) can be written as  
\[\mathcal{P}_{\omega}=\frac{1}{\tau^{2}}(\Gamma_{\omega}^{-1}F_{n}^{*}\otimes I _{m})(D_{1}\otimes L-2D_{2}\otimes I_{m})(F_{n}\Gamma_{\omega}\otimes I_{m}),\]  
where \(F_{n}\) is the Fourier matrix, \(\Gamma_{\omega}=\text{diag}\Big{(}1,\omega^{\frac{1}{2}},\ldots,\omega^{\frac {n-1}{n}}\Big{)}\), \(D_{1}=\text{diag}\Big{(}\sqrt{n}F_{n}\Gamma_{\omega}\Psi_{n,\omega}^{(1)} \Big{)}\) and \(D_{2}=\text{diag}\Big{(}\sqrt{n}F_{n}\Gamma_{\omega}\Pi_{n,\omega}^{(1)} \Big{)}\), with \(\Psi_{n,\omega}^{(1)}\), \(\Pi_{n,\omega}^{(1)}\) being the first column of \(\Psi_{n,\omega}\) and \(\Pi_{n,\omega}\). Therefore, given a vector \(\mathbf{r}\), the proposed algorithm for computing \(\mathbf{z}=\mathcal{P}_{\omega}^{-1}\mathbf{r}\) is the following  
1. Compute \(\mathbf{s}_{1}=(F_{n}\Gamma_{\omega}\otimes I_{m})\mathbf{r}\) via the FFT;
2. Compute \(\mathbf{s}_{2}^{(k)}=\tau^{2}\Big{(}D_{1}(k,k)L-2D_{2}(k,k)I_{m}\Big{)}^{-1} \mathbf{s}_{1}^{(k)},\quad k=1,\ldots,n\), where \(\mathbf{s}_{1}^{(k)}\), \(\mathbf{s}_{2}^{(k)}\) denote the \(k\)-th block of dimension \(m\) of \(\mathbf{s}_{1}\), \(\mathbf{s}_{2}\);
3. Compute \(\mathbf{z}=(\Gamma_{\omega}^{-1}F_{n}^{*}\otimes I_{m})\mathbf{s}_{2}\) via the FFT;  
Which in its essence is identical to the one proposed by [PERSON], although we must observe that Ref. [8] is not referenced in the work by [PERSON] and [PERSON] [21]. We conclude this section by discussing the matter of stability when dealing with \(\omega\)-circulants. While theoretically the aforementioned algorithm may be applied whenever \(\omega\
eq 0\), in practical implementation exceedingly small values of the parameter will cause stability issues, which should be considered, when stating the global precision of a solution method. In fact, as we mentioned in the previous section, the condition number of the transform \(F_{n}\Gamma_{\omega}\) is equal to that of \(\Gamma_{\omega}\), which tends to \(\max\{|\omega|,|\omega|^{-1}\}\) as \(n\to\infty\) and therefore tends to infinity as \(\omega\to 0\). We tested the stability of the procedure described above by fixing a vector \(\mathbf{x}\), computing \(\mathbf{b}=\mathcal{P}_{\omega}\mathbf{x}\) and retrieving the starting vector as \(\mathcal{P}_{\omega}^{-1}\mathbf{b}\) with the spectral decomposition of \(\mathcal{P}_{\omega}\), using the algorithm described above and in Ref. [21]. We denote the retrieved vector with \(\bar{\mathbf{x}}\). As a first example, we set \(m=2^{8}\), \(n=2^{8}\), \(\mathbf{x}=\left[\sin(\frac{j\tau}{nm}+3)\right]_{j=1,\ldots,mm}\) and show the relative errors between \(\mathbf{x}\) and \(\bar{\mathbf{x}}\) as the parameter \(\omega\) decreases in Table 1. Clearly, as \(\omega\) gets closer to zero, the error grows and \(\bar{\mathbf{x}}\) ceases to be an accurate approximation of \(\mathbf{x}\), especially in the right part in accordance with the accumulation errors in the solution of a lower triangular system. In particular the error for very small values of \(\omega\) becomes even worse than that obtained when using a direct inherently sequential algorithm working for triangular systems (\(\omega=0\)). Figure 1 shows the plots of the two vectors for \(\omega=10^{-16}\). Our second example, where \(\mathbf{x}\) is this time chosen at random, grants similar results, collected in Table 2. We observe that this matter is not addressed in Ref. [21]. Because of these instability issues, the possibility of using preconditioned MINRES recently emerged, where the preconditioner is taken from the \(\tau\) algebra whose eigenvector matrix can be chosen as a real symmetric orthogonal and very stable transform (see Ref. [22] and the references therein).