## 1 Introduction  
Eigenvalue analysis is a fundamental component of electric power system Small-Signal Stability Analysis (SSSA). It is utilized in order to determine the stability of the operating points of the system [1; 2], design controllers [3], determine machine clusters and reduced equivalent models [4], and evaluate the properties of the transmission grid [5]. On the other hand, there is rich literature on numerical algorithms that determine the full or a partial solution of a given eigenvalue problem. Relevant monographs on the topic are, for example [6; 7]. However, not all available algorithms are suitable for the SSSA of power systems. The vast majority of numerical algorithms, in fact, solve exclusively symmetric, Hermitian, or tridiagonal Hermitian eigenvalue problems, which arise in many engineering, physics, and computer applications. Such algorithms are, for example, the Davidson [8], the implicitly restarted Lanczos [9], and the locally optimal block preconditioned conjugate gradient [10] methods. However, the matrices that describe a linearized power system model are typically non-symmetric. When compared to symmetric problems, non-symmetric eigenvalue problems are more difficult and computationally demanding to solve. This is also due to the fact that symmetric eigenvalue problems are generally well-conditioned, while non-symmetric eigenvalue problems are not. Consequently, the available free and open-source software libraries that solve non-symmetric eigenvalue problems are a small subset of all existing libraries. The scalability of the numerical solution of eigenvalue problems is also very important, since real-world electrical power networks are large-scale dynamic systems. Unfortunately, the most reliable methods for finding the full spectrum of an eigenvalue problem are dense-matrix methods, and their computational complexity and memory requirements increase more than quadratically (in some cases even cubically) as the size of the matrix increases. As a matter of fact, the largest ever eigenvalue analysis with a dense algorithm to date was the solution of a \(10^{6}\times 10^{6}\) problem in about 1 h, and it was carried out in 2014 by the Japanese K computer in Riken. To be able to obtain this result, the K computer includes 88,000 processors that draw a peak power of 12.6 MW, while its operation costs annually US$10 million. The solution of large-scale power system eigenvalue problems is challenging even when using sparse matrices and limiting the search to a subset of the spectrum. A coarse taxonomy of existing algorithms for the solution of non-symmetric eigenvalue problems is as follows:  
* Vector iteration methods, which, in turn, are separated to single and simultaneous vector iteration methods. Single vector iteration methods include the power method [11] and its variants, such as the inverse power and Rayleigh quotient iteration. Simultaneous vector iteration methods include the subspace iteration method [12] and its variants, such as the inverse subspace method. Regarding the application of vector iteration methods in power systems, the inverse power iteration was first discussed in [13]. A Rayleigh quotient iteration based algorithm was proposed in [14; 15], while recent studies have provided subspace accelerated and deflated versions of the same algorithm [16; 17; 18]. Finally, simultaneous vector iteration methods were studied in [19; 20]. * Schur decomposition methods, which mainly include the QR algorithm [21], the QZ algorithm [22], and their variants, such as the QR algorithm with shifts. Schur decomposition based methods have been the standard methods employed for the eigenvalue analysis of small to medium size power systems [23; 24]. * Krylov subspace methods, which basically include the Arnoldi iteration [25] and its variants, such as the implicitly restarted Arnoldi [26] and the Krylov-Schur method [27].