## 2 Matrix Sign Function-Based Solvers
```
\(t_{0}\gets a\) for\(k\gets 0\) to \(\lceil\frac{b-a}{h}\rceil\) \(\tilde{A}_{k}\gets A-B(Z_{k}(Z_{k}^{T}B))^{T}-\frac{1}{2h}I\) \(N_{k}\leftarrow[\,C^{T}\ Z_{k}(Z_{k}^{T}B)\ \sqrt{h^{-1}Z_{k}}\,]\)  Apply Algorithm 2--GECINC to obtain \(Z_{k+1}\), a low-rank approximation to  the solution of \(\tilde{A}_{k}^{T}X_{k+1}+X_{k+1}\tilde{A}_{k}=-N_{k}N_{k}^{T}\) \(t_{k+1}\gets t_{k}+h\) endfor
```
**Algorithm 4****ROS1****:**

Note that, since we obtain the low rank factor of the solution of Equation (8), the cost associated with the updates of both \(\tilde{A}_{k}\) and \(\tilde{N}_{k}\) is drastically decreased, and thus, the most time-consuming operation in Algorithm 4--ROS1 is the solution of Equation (8), _i.e._, the execution of Algorithm 2--GECINC. ## 3 High Performance Implementations

In this section, we briefly review several high performance implementations for the solution of the matrix equations in the study. In particular, we describe Lyapunov, ARE and DRE solvers based on Algorithm 2--GECINC, Algorithm 3--GECRSG and Algorithm 4--ROS1, respectively. These algorithms require the computation of the sign function and, therefore, inversion of a large dense matrix, which in all cases represents the most expensive operation from the computational point of view. In the next two subsections, we describe an efficient and reliable numerical method for matrix inversion and the parallelization techniques applied in each platform: multi-core, hybrid CPU-GPUsystems and hybrid multi-GPU platforms. In Section 3.3, we then discuss some additional details on the parallel solution of matrix equations. ### Matrix Inversion

Traditional matrix inversion via Gaussian elimination requires the computation of the LU factorization of the matrix and the solution of two triangular linear systems. This method presents two drawbacks when implemented in parallel architectures: the process is partitioned into three consecutive stages, and the solution of the linear systems involves triangular matrices (impairing load-balance). ```
Partition \(A\rightarrow\left( \begin{array}{c|c}A_{TL}&A_{TR}\\ \hline A_{BL}&A_{BR}\end{array} \right),\) where \(A_{TL}\) is \(0\times 0\) and \(A_{BR}\) is \(n\times n\) while\(m(A_{TL})<m(A)\)do Determine block size \(b\) Repartition \(\left( \begin{array}{c|c}A_{TL}&A_{TR}\\ \hline A_{BL}&A_{BR}\end{array} \right)\rightarrow\left( \begin{array}{c|c}A_{0 0}&A_{01}&A_{02}\\ \hline A_{10}&A_{11}&A_{12}\\ \hline A_{20}&A_{21}&A_{22}\end{array} \right),\) where \(A_{11}\) is \(b\times b\) \(\left[ \begin{array}{c|c}A_{01}\\ A_{11}\\ A_{21}\end{array} \right]\) \(\leftarrow\) GJE\({}_{\text{UNB}}\left(\left[ \begin{array}{c}A_{01}\\ A_{11}\\ A_{21}\end{array} \right]\right)\) \(A_{00}\) \(\leftarrow\)\(A_{00}+A_{01}A_{10}\) \(A_{20}\) \(\leftarrow\)\(A_{20}+A_{21}A_{10}\) \(A_{10}\) \(\leftarrow\)\(A_{11}A_{10}\) \(A_{02}\) \(\leftarrow\)\(A_{02}+A_{01}A_{12}\) \(A_{22}\) \(\leftarrow\)\(A_{22}+A_{21}A_{12}\) \(A_{12}\) \(\leftarrow\)\(A_{11}A_{12}\) Continue with \(\left( \begin{array}{c|c}A_{TL}&A_{TR}\\ \hline A_{BL}&A_{BR}\end{array} \right)\leftarrow\left( \begin{array}{c|c}A_{0 0}&A_{01}&A_{02}\\ \hline A_{10}&A_{11}&A_{12}\\ \hline A_{20}&A_{21}&A_{22}\end{array} \right)\) endwhile
```

**Algorithm 5**GJEBLK':

On the other hand, matrix inversion via the Gauss-Jordan elimination method (GJE) is mathematically equivalent and presents the same arithmetic cost to/as the Gaussian-based counterpart, but computes the inverse in a single stage, does not require to operate with triangular matrices and casts the bulk of the computations in terms of coarse-grain highly parallel operations (matrix-matrix products). Algorithm 5--GJE\({}_{\text{BLK}}\), presented above, illustrates a blocked version of the matrix inversion GJE-based algorithm using a FLAME (Formal Linear Algebra Methods Environment) notation [16]. There, \(m(A)\) stands for the number of rows of the matrix. We believe the rest of the notation to be intuitive; for further details, see [16].