A variational method is developed to determine configuration interaction wave functions. The method is straightforward and is applied to a pairing Hamiltonian with constant matrix elements, for which exact eigenvalues are available. Comparisons are made with the exact results. Calculations can be carried out to any desired degree of accuracy. The method is also applied to a Hamiltonian that has neutron-neutron, proton-proton, and neutron-proton pairing. No difficulties are found in extending the method to this Hamiltonian that has many collective modes. In practice, the method scales linearly with , where is the number of variational basis states.