Several modifications to the Davidson algorithm are systematically explored to establish their performance for an assortment of configuration interaction (CI) computations. The combination of a generalized Davidson method, a periodic two-vector subspace collapse, and a blocked Davidson approach for multiple roots is determined to retain the convergence characteristics of the full subspace method. This approach permits the efficient computation of wave functions for large-scale CI matrices by eliminating the need to ever store more than three expansion vectors (bi) and associated matrix-vector products (sigmai), thereby dramatically reducing the I/O requirements relative to the full subspace scheme. The minimal-storage, single-vector method of Olsen is found to be a reasonable alternative for obtaining energies of well-behaved systems to within microhartree accuracy, although it typically requires around 50% more iterations and at times is too inefficient to yield high accuracy (ca. 10-10 hartree) for very large CI problems. Several approximations to the diagonal elements of the CI Hamiltonian matrix are found to allow simple on-the-fly computation of the preconditioning matrix, to maintain spin symmetry of the determinant-based wave function, and to preserve the convergence characteristics of the diagonalization procedure.