An efficient direct integral-driven algorithm for the random-phase approximation (RPA) is introduced using the equation of motion on a transition density matrix representing a "quasi particle". In the algorithm, several roots are obtained at the same time by solving a set of coupled equations that are projected on a space spanned by a set of error vectors representing the quasi particles. The most time-consuming RPA operation on the vectors is accomplished by a single call of an integral-generation routine, and the time per iteration is comparable to that for a direct SCF cycle. The algorithm is implemented using a new integral package based on accompanying coordinate expansion (ACE), as well as traditional integral routines from GAMESS. The example applications indicate good convergence of the iterative scheme. In some computationally intensive cases, the RPA computation for several excited states is completed in less time than the SCF computation.