Abstract: A matrix-forming algorithm has been developed for the massively parallel solution of multi-dimensional non-symmetric eigenvalue and singular value decomposition problems governing global linear fluid flow instability. In the high supersonic and hypersonic regimes of interest, the EVP is described by dense matrices, while the SVD destroys any potentially existing sparsity pattern. The ScaLAPACK library, available as part of vendor-optimized libraries on most modern supercomputers, is used to perform a full distributed LU-decomposition, as required by shift-and-invert strategies which promote physically-interesting and suppress spurious eigenmodes. The iterative generation of Krylov subspaces and the computation of the Ritz vectors are performed by linear algebra operations distributed over the available processors, while the serial computation of the Ritz values represents a negligible part of the overall computing effort. The algorithm is validated by solving two- and three-dimensional eigenvalue problems, results of which are known analytically, as well as BiGlobal eigenvalue problems in incompressible and hypersonic flow configurations. The proposed methodology sets new standards regarding the size of problems that can now be addressed by global modal or non-modal linear stability analysis techniques.