In this article, we propose a fast method to solve large-scale three-dimensional topology optimization problems subject to buckling constraints. Buckling analysis entails the solution of a generalized eigenvalue problem. For problems with large degrees of freedom, the current numerical methods tend to be memory-hungry, leading to high computational costs. First, a low-memory assembly-free linear buckling analysis method is proposed. Specifically, this method is based on the voxelization model, an assembly-free version of the deflated conjugate gradient is used to accelerate the iteration solution of linear systems of equations, where neither the stiffness matrix nor the deflation matrix is assembled, and the parallelization of matrix–vector multiplication is achieved by the congruency of voxels. Due to the particularity of the stress stiffness matrix in the buckling analysis, the inverse iteration is used to solve the general eigenvalue problem, which can reduce the operations of stress stiffness matrix considerably. Based on the efficient buckling analysis method, we extend the level-set method for buckling constraints in a semi-analytical manner. Several numerical experiments demonstrate that the proposed method can solve large-scale three-dimensional buckling analysis and topology optimization against buckling constraints effectively.