If \(G\) is the transfer function of a stable linear time-invariant system and suppose \(G(s)G^T(-s)=W^T(-s)W(s)\), then a reduced order model \(G_r\) can be obtained by truncating a balanced realization of \(G\). For balancing, one needs the controllability grammian \(P\) of \(G\) and the observability grammian \(Q\) of \(W\). This \(P\) is the solution of a Lyapunov equation (LE), and \(Q\) is the stabilizing solution of an algebraic Riccati equation (ARE). Because of the numerical conditioning, it is better to work with square roots of \(P\) and \(Q\) instead of with \(P\) and \(Q\) themselves. Square roots can be obtained for example as Cholesky factors. However, for large scale systems, these factors tend to be of relatively low rank, so that it is more efficient to work with rectangular square roots that are of full rank. This paper gives a detailed description of efficient, reliable and parallelizable numerical implementations of the different steps of such an algorithm for large scale model reduction that is based on truncated balanced realization. The LEs are solved by a sign-function technique, and the same technique is used in the inner loop of the Newton iteration to solve the ARE. Numerical examples show the effectiveness of the method, which allows problems with a couple of thousands of state variables to be solved.