Given a large square real matrix A and a rectangular tall matrix Q, many application problems require the approximation of the operation exp(A)Q\exp(A)Q. Under certain hypotheses on A, the matrix exp(A)Q\exp(A)Q preserves the orthogonality characteristics of Q; this property is particularly attractive when the associated application problem requires some geometric constraints to be satisfied. For small size problems numerical methods have been devised to approximate exp(A)Q\exp(A)Q while maintaining the structure properties. On the other hand, no algorithm for large A has been derived with similar preservation properties. In this paper we show that an appropriate use of the block Lanczos method allows one to obtain a structure preserving approximation to exp(A)Q\exp(A)Q when A is skew-symmetric or skew-symmetric and Hamiltonian. Moreover, for A Hamiltonian we derive a new variant of the block Lanczos method that again preserves the geometric properties of the exact scheme. Numerical results are reported to support our theoretical findings, with particular attention to the numerical solution of linear dynamical systems by means of structure preserving integrators.