We propose a fast and versatile algorithm to calculate local and transport properties such as conductance, shot noise, local density of state or local currents in mesoscopic quantum systems. Within the non equilibrium Green function formalism, we generalize the recursive Green function technique to tackle multiterminal devices with arbitrary geometries. We apply our method to analyze two recent experiments: an electronic Mach-Zehnder interferometer in a 2D gas and a Hall bar made of graphene nanoribbons in quantum Hall regime. In the latter case, we find that the Landau edge state pinned to the Dirac point gets diluted upon increasing carrier density.The field of quantum transport at the nanometer scale now includes a large number of systems involving very different physics. Examples include for instance mesoscopic devices in two dimensional heterostructures [1], graphene nanoribbons [2], superconducting weak links [3], molecular electronic devices [4] or ferromagnetic multilayers nanopillars [5,6]. Although those systems have different structures and geometries, they are all quantum systems connected to the macroscopic world through electrodes and, consequently, formalisms developed to describe one of them can often be adapted to the others. This is in particular the case of the widely used Landauer-Büttiker formalism [7] which focuses on the scattering properties of the system. The formalism is very intuitive and general. However, it is not well suited when one is interested in what happens inside the sample or for performing a microscopic calculation for a given device. An alternative mathematically equivalent approach is referred as the NEGF (non equilibrium Green function) formalism [8,9]. NEGF, which is derived from the Keldysh formalism [10], provides a simple route to compute the physical observables from a microscopic model. It is now an extremely popular numerical approach to a very wide class of physical problems (see references in [11]). For instance, all the references mentioned in examples above correspond to calculation done with this technique either from ab-initio or from phenomenological models [1,2,3,4,5,6,12,13]. At the core of NEGF is the calculation of the retarded Green function G of the mesoscopic region in presence of the (semiinfinite) electrodes. A straightforward method consists of direct inversion of the Hamiltonian H. However, when doing so, one is restricted to rather small systems of a few thousand sites: for a system of size L in dimension d, the computing time scales as L 3d while the needed memory scales as L2d . An alternative algorithm, known as the recursive Green function technique [1,5,14], takes advantage of the structure of H to reduce drastically the computing time down to L 3d−2 , putting systems of a few million sites within reach. In its original version [14,15,16], only the transport properties of the device could be computed, but recent progresses made it possible to get access to observables inside the sample (like local electronic density or local currents [1,17]...