An efficient and numerically stable eigensolution method for structures with multiple natural frequencies is presented. The proposed method is developed by improving the well-known subspace iteration method with shift. A major difficulty of the subspace iteration method with shift is that because of singularity problem, a shift close to an eigenvalue can not be used, resulting in slower convergence. In this paper, the above singularity problem has been solved by introducing side conditions without sacrifice of convergence. The proposed method is always nonsingular even if a shift is on a distinct eigenvalue or multiple ones. This is one of the significant characteristics of the proposed method. The nonsingularity is proved analytically. The convergence of the proposed method is at least equal to that of the subspace iteration method with shift, and the operation counts of above two methods are almost the same when a large number of eigenpairs are required. To show the effectiveness of the proposed method, two numerical examples are considered.