The integration of surface normals for the purpose of computing the shape of a surface in 3D space is a classic problem in computer vision. However, even nowadays it is still a challenging task to devise a method that is flexible enough to work on non-trivial computational domains with high accuracy, robustness, and computational efficiency. By uniting a classic approach for surface normal integration with modern computational techniques, we construct a solver that fulfils these requirements. Building upon the Poisson integration model, we use an iterative Krylov subspace solver as a core step in tackling the task. While such a method can be very efficient, it may only show its full potential when combined with suitable numerical preconditioning and problem-specific initialisation. We perform a thorough numerical study in order to identify an appropriate preconditioner for this purpose. To provide suitable initialisation, we compute this initial state using a recently developed fast marching integrator. Detailed numerical experiments illustrate the benefits of this novel combination. In addition, we show on real-world photometric stereo datasets that the developed numerical framework is flexible enough to tackle modern computer vision applications.