Global stability modes of flows provide significant insight into their dynamics. Direct methods to obtain these modes are restricted by the daunting sizes and complexity of Jacobians encountered in general three-dimensional flows. Jacobian-free approaches have greatly alleviated the required computational burden. Particularly, the most common Arnoldi-based methods obtain the desired subset of the eigenmodes by considering Jacobian-vector products to create a smaller iterative subspace, instead of working with the Jacobian itself. However, operations such as orthonormalization and shift-and-invert transformation of matrices with appropriate shift guesses can introduce computational and parameter-dependent costs that inhibit their routine application to general three-dimensional flowfields. Further, in time-stepper type approaches, where the linearized perturbation snapshots directly obtained from an aerodynamic code are treated as Jacobian-vector products, the inversion operation necessitates use of approximate iterative linear solvers with several parameters. The present work addresses these limitations by proposing and implementing a robust, generalizable approach to extract the principal global modes, suited for curvilinear coordinates as well as the effects of compressibility. Accurate linearized perturbation snapshots are obtained using high-order schemes by leveraging the same non-linear Navier-Stokes code as used to obtain the basic state by appropriately constraining the equations using a body-force. The forcing includes not only that required to obtain the products, but also to ensure that the basic state does not drift.It is shown that with random impulse forcing, dynamic mode decomposition (DMD) of the subspace formed by these products yields the desired physically meaningful modes, when appropriately scaled. The leading eigenmodes are thus obtained without spurious modes or the need for an iterative procedure. Further since orthonormalization is not required, large subspaces can be processed to capture converged low frequency or stationary modes. The validity and versatility of the method is demonstrated with numerous examples encompassing essential elements expected in realistic flows, such as compressibility effects and complicated domains requiring general curvilinear meshes. Favorable qualitative as well as quantitative comparisons with Arnoldi-based method, complemented with substantial savings in computational resources show the potential of current approach for relatively complex flows.