A select-divide-and-conquer variational method to approximate configuration interaction ͑CI͒ is presented. Given an orthonormal set made up of occupied orbitals ͑Hartree-Fock or similar͒ and suitable correlation orbitals ͑natural or localized orbitals͒, a large N-electron target space S is split into subspaces S 0 , S 1 , S 2 , ... ,S R . S 0 , of dimension d 0 , contains all configurations K with attributes ͑energy contributions, etc.͒ above thresholds T 0 ϵ͕T 0 egy , T 0 etc.͖; the CI coefficients in S 0 remain always free to vary. S 1 accommodates Ks with attributes above T 1 ഛ T 0 . An eigenproblem of dimension d 0 + d 1 for S 0 + S 1 is solved first, after which the last d 1 rows and columns are contracted into a single row and column, thus freezing the last d 1 CI coefficients hereinafter. The process is repeated with successive S j ͑j ജ 2͒ chosen so that corresponding CI matrices fit random access memory ͑RAM͒. Davidson's eigensolver is used R times. The final energy eigenvalue ͑lowest or excited one͒ is always above the corresponding exact eigenvalue in S. Threshold values ͕T j ; j =0,1,2, ... ,R͖ regulate accuracy; for large-dimensional S, high accuracy requires S 0 + S 1 to be solved outside RAM. From there on, however, usually a few Davidson iterations in RAM are needed for each step, so that Hamiltonian matrix-element evaluation becomes rate determining. One hartree accuracy is achieved for an eigenproblem of order 24ϫ 10 6 , involving 1.2ϫ 10 12 nonzero matrix elements, and 8.4ϫ 10 9 Slater determinants.