In this paper we propose an iterative mean algorithm involving arithmetic and geometric means of n positive definite matrices which generalizes the 3-dimensional algorithm of positive reals discovered by Carlson (1970) [10]. We show that the iterative mean algorithm is convergent and the common limit satisfies multidimensional versions of all properties (permutation symmetry, concavity, monotonicity, homogeneity, congruence invariancy, duality, mean inequalities) that one would expect for the Carlson mean of positive reals. Convergence and perturbation analysis with numerical experiments are presented in terms of the Thompson metric and the spectral norm.