In this work, we present the implementation of a variational density fitting methodology that uses iterative linear algebra for solving the associated system of linear equations. It is well known that most difficulties with this system arise from the fact that the coefficient matrix is in general ill-conditioned and, due to finite precision round-off errors, it may not be positive definite. The dimensionality, given by the number of auxiliary functions, also poses a challenge in terms of memory and time demand since the coefficient matrix is dense. The methodology presented is based on a preconditioned Krylov subspace method able to deal with indefinite ill-conditioned equation systems. To assess its potential, it has been combined with double asymptotic electron repulsion integral expansions as implemented in the deMon2k package. A numerical study on a set of problems with up to 130,000 auxiliary functions shows its effectiveness to alleviate the abovementioned problematic. A comparison with the default methodology used in deMon2k based on a truncated eigenvalue decomposition of the coefficient matrix indicates that the proposed method exhibits excellent robustness and scalability when implemented in a parallel setting.