We develop a data-sparse and accurate approximation of the normalised hyperbolic operator sine family generated by a strongly P-positive elliptic operator defined in [Gavrilyuk, J. Math. Anal. Appl. 236: 327–349, 1999, Gavrilyuk and Makarov, Numer. Funct. Anal. Optimization 20: 695–717, 1999]. In the preceding papers [Hackbusch, Computing 62: 89–108, 1999–Hackbusch, Khoromskij, and Sauter, On ℋ 2 -matrices, Springer-Verlag, 2000], a class of ℋ-matrices has been analysed which is data-sparse and allow an approximate matrix arithmetic with almost linear complexity. An ℋ-matrix approximation to the operator exponent with a strongly P-positive operator was proposed in [Gavrilyuk, Hackbusch, and Khoromskij, ℋ-matrix approximation for the operator exponential with applications, 2000]. In the present paper, we apply the ℋ-matrix techniques to approximate the elliptic solution operator on cylinder-like domains Ω × [ a, b ] associated with the elliptic operator , x ∈ [ a, b ]. It is explicitly presented by the normalised hyperbolic sine of an elliptic operator ℒ defined in Ω. The approach is then applied to elliptic partial differential equations in domains composed of rectangles or cylinders. In particular, we consider the ℋ-matrix approximation to the interface Poincaré-Steklov operators with application in the Schur-complement domain decomposition method. Starting with the Dunford–Cauchy representation for the hyperbolic sine operator, we then discretise the integral by the exponentially convergent quadrature rule involving a short sum of resolvents. The latter are approximated by the ℋ-matrix techniques. Our algorithm inherits a two-level parallelism with respect to both the computation of resolvents and the treatment of different values of the spatial variable x ∈ [ a, b ].