We give an algorithmic approach to the approximative solution of operator Lyapunov equations for controllability. Motivated by the successfully applied alternating direction implicit (ADI) iteration for matrix Lyapunov equations, we consider this method for the determination of Gramian operators of infinite-dimensional control systems. In the case where the input space is finite-dimensional, this method provides approximative solutions of finite rank. Under the assumption of infinite-time admissibility and boundedness of the semigroup, we analyze convergence in several operator norms. We show that under a mild assumption on the shift parameters, convergence to the Gramian is obtained. Particular emphasis is placed on systems governed by a heat equation with boundary control. We present that ADI iteration for the heat equation consists of solving a sequence of Helmholtz equations. Two numerical examples are presented; the first showing the benefit of adaptive finite elements and the second illustrating convergence to something other than the Gramian in a case where our condition on the shift parameters is not satisfied.