Primary instability of the lid-driven flow in a cube is studied by a linear stability approach.Two cases, in which the lid moves parallel to the cube sidewall or parallel to the diagonal plane, are considered. It is shown that Krylov vectors required for application of the Newton and Arnoldi iteration methods can be evaluated by the SIMPLE procedure. The finite volume grid is gradually refined from 100 3 to 256 3 nodes. The computations result in grid converging values of the critical Reynolds number and oscillation frequency that allow for Richardson extrapolation to the zero grid size. Three-dimensional flow and most unstable perturbations are visualized by a recently proposed approach that allows for a better insight in the flow patterns and appearance of the instability. New arguments regarding the assumption that the centrifugal mechanism triggers the instability are given for both cases.Keywords: lid-driven cavity flow, Newton method, Arnoldi method, linear stability, Krylovsubspace iteration, SIMPLE steady-oscillatory transition was studied recently in [22], again by the straightforward integration in time.Our primary goal is to apply the eigenvalue analysis for calculation of the critical Reynolds numbers, at which the steady flow bifurcates to an oscillatory state. The computations are carried out using gradually refined stretched grids, containing 100 3 , 150 3 , 200 3 , and 256 3 finite volumes. In spite that pseudo-spectral and collocation methods may yield better accuracy for model problems in simple geometries, as those considered here, they are not applicable to most of the applied problems that involve more complicated flow regions and liquid-liquid interfaces. Besides a wide set of curved-boundary-fitted numerical methods, such problems can be treated by the fixed grid finite volume approach used below with the immersed boundary technique, as it was done recently in [23].The critical Reynolds numbers are calculated together with the critical oscillation frequencies. These are defined by the imaginary part of the leading eigenvalue, whose real part crosses the imaginary axis at the stability limit. The most unstable infinitesimally small perturbation is defined by the eigenvector corresponding to the leading eigenvalue. The eigenvector spatial pattern allows us to compare with the results of previous studies that applied the straightforward integration in time, to visualize spatio-temporal behavior of the most unstable disturbance, and basing on this to speculate further about the reasons that cause the instability onset.In the following we formulate the problem and briefly describe the numerical method, as well as our method of visualization of three-dimensional divergence-free velocity field [24,25]. The present numerical approach includes Newton iteration for calculation of the steady states and Arnoldi iteration for computation of leading eigenvalues. The Newton corrections are computed by the Krylov-subspace-iteration techniques BiCGstab(2) or GMRES. The SIMPLE procedure [26] is reformulated ...