Model order reduction (MOR) is common in simulation, control and optimization of complex dynamical systems arising in modeling of physical processes and in the spatial discretization of parabolic partial differential equations (PDEs) in two or more dimensions. Typically, after a semidiscretization of the differential operator by the finite element method (FEM) or by the boundary element method, we have a large state-space dimension n = O(104). It is assumed that the number of inputs and outputs is equal and much smaller than n. We show how to compute an approximate reduced-order system of order r n with a balancing-related model reduction method. The method is based on the computation of the cross-Gramian (CG) X, which is the solution of one Sylvester equation. As standard algorithms for the solution of Sylvester equations are of limited use for largescale systems, we investigate approaches based on the sign function method. To make this iterative method applicable in the large-scale ...