## 2 Matrix Sign Function-Based Solvers  
The sign function method [2] is an efficient numerical tool to solve Lyapunov, Sylvester and Riccati equations, with dense coefficient matrices, on parallel message-passing computers [3; 4], as well as on hardware accelerators [10]. The convenience of the sign function method is based on two properties: First, it is composed of well-known basic linear algebra operations that exhibit a high degree of concurrency. Moreover, high performance implementations for parallel architectures of these operations are included in linear algebra libraries, like BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage), and their extension for GPUs, CUBLAS (CUDA BLAS) from NVIDIA and for message-passing platforms, e.g., ScaLAPACK (Scalable LAPACK). Second, it is an iterative algorithm that presents a fast convergence rate that is asymptotically quadratic. Several schemes have been proposed in the literature to compute the sign function. Among them, the Newton iteration, illustrated in Algorithm 1--GESINE below, exhibits a remarkable simplicity and efficiency. ```
\(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{k+1}\leftarrow\left(A_{k}+A_{k}^{-1}\right)/2\) \(k\gets k+1\) until\(\sqrt{\|A_{k}-A_{k-1}\|_{1}}<\tau_{s}\|A_{k}\|_{1}\)
```  
**Algorithm 1**GESINE**:  
```
\(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\gets 0\) repeat \(A_{0}\gets A\) \(k\The most time-consuming operation in Algorithm 1--GESINE is the inversion of \(A_{k}\). Given a square matrix, \(A\), of order \(n\), this operation renders the cost of the algorithm as \(2n^{3}\) flops per iteration. The cubic computational cost in the matrix dimension is inherited by all solvers described next. To avoid stagnation of the iteration, we set \(\tau_{s}=n\cdot\sqrt{\varepsilon}\), where \(\varepsilon\) stands for the machine precision, and perform one additional iteration step after the stopping criterion is satisfied. Due to the asymptotic quadratic convergence rate of the Newton iteration, this is usually enough to reach the attainable accuracy. ### Solution of Lyapunov Equations

Algorithm 2--Geclnc, presented below, illustrates a variant of the sign function method for the solution of a Lyapunov equation of the form:

\[AX+XA^{T}\ =\ -BB^{T} \tag{1}\]

where \(A\in\mathbb{R}^{n\times n}\) is c-stable (_i.e._, all its eigenvalues have a negative real part), \(B\in\mathbb{R}^{n\times m}\) and \(X\in\mathbb{R}^{n\times n}\) is the desired solution. ```
\(A_{0}\gets A\), \(\tilde{S}_{0}\gets B^{T}\) \(k\gets 0\) repeat  Compute the rank-revealing QR (rank-revealing QR (RRQR)) decomposition \(\frac{1}{\sqrt{2c_{k}}}\left[ \begin{array}{cc}\tilde{S}_{k},&c_{k}\tilde{S} _{k}A_{k}^{-T}\end{array} \right]=Q_{s}\left[ \begin{array}{c}U_{s}\\ 0\end{array} \right]\Pi_{s}\) \(\tilde{S}_{k+1}\gets U_{s}\Pi_{s}\) \(A_{k+1}\leftarrow\frac{1}{\sqrt{2}}\left(A_{k}/c_{k}+c_{k}A_{k}^{-1}\right)\) \(k\gets k+1\) until\(\sqrt{\left\|A_{k}-I\right\|}_{1}<\tau_{l}\)
```
**Algorithm 2****Geclnc:**

On convergence, after \(\tilde{k}\) iterations, \(\tilde{S}=\frac{1}{\sqrt{2}}\tilde{S}_{k}\) is a full (row-)rank approximation of \(S\), so that \(X=S^{T}S\approx\tilde{S}^{T}\tilde{S}\). In practice, the scaling factor, \(c_{k}\), is used to accelerate the convergence rate of the algorithm, (_i.e._, to reduce the number of required iterations). In our case, we set:

\[c_{k}=\|A_{k}\|/\|A_{k}^{-1}\|\]

Furthermore, we choose \(\tau_{l}=\tau_{s}\) and perform an extra step after the convergence criterion is satisfied. Note that the number of columns of \(\tilde{S}_{k}\) is doubled at each iteration and, in consequence, the computational and storage costs associated with its update increase with each iteration.