We present a general variational method to calculate vibrational energy levels of polyatomic molecules without dynamical approximation. The method is based on a Lanczos algorithm, which does not require storage of the Hamiltonian matrix. The rate-determining step of each Lanczos iteration is the evaluation of the product of the matrix and a trial vector. We use simple product basis functions and write the Hamiltonian as a sum of factorizable terms. With n one-dimensional functions in each of f dimensions, the matrix-vector product requires no more than cnf+1 multiplications for a single term involving c coordinates. Choosing a (potential optimized) discrete variable representation (DVR) in each dimension, the potential energy matrix is diagonal. The rate-determining step is now the multiplication of a vector by the kinetic energy matrix and c is effectively (with rare exceptions) at most two. The nf+1 scaling holds for both diagonal and mixed second derivative operators. The method is directly applicable to any three-atom and any nonlinear four-atom molecule. We use a variety of coordinate systems (Jacobi, Radau, a hybrid of the two, and bond), for which the total number of factorizable terms in the exact kinetic energy operator is never large, to calculate very well-converged band origins of H2O up to 22 000 cm−1, of H+3 up to 18 000 cm−1, and of CH2O up to 5700 cm−1; and low-lying levels of H2O2. The results for CH2O are new, and those for H+3 clarify the causes of discrepancies in published work. The product basis results in very large matrices (up to 500 000×500 000 for four atoms), but the cost is within an order of magnitude of that of contracted-basis approaches using explicit diagonalization. While contracted basis approaches are molecule and Hamiltonian specific, it was possible to apply the DVR-Lanczos method to all the examples presented here with a single computer program. The principal advantage of our method is thus its generality, and in this context it is efficient, with the cost scaling slowly with basis size. It is also easily parallelized.