We present an efficient implementation of an approximate balanced truncation model reduction technique for general large-scale RLC systems, described by a state-space model where the "C" matrix in the time-domain modified nodal analysis (MNA) circuit equation "C\dot{x}=-Gx+Bu" is not necessarily invertible. The large sizes of the models that we consider make most implementations of the balance-and-truncate method impractical from the points of view of computational load and numerical conditioning. This motivates our use of Krylov subspace methods to directly compute approximate low-rank square roots of the Gramians of the original system. The approximate low-order general balanced and truncated model can then be constructed directly from these square roots. We demonstrate using three practical circuit examples that our new approach effectively gives approximate balanced and reduced order coordinates with little truncation error.