When constructing multivariate Pad approximants, highly structured linear systems arise in almost all existing definitions . Until now little or no attention has been paid to fast algorithms for the computation of multivariate Pad approximants, with the exception of . In this paper we show that a suitable arrangement of the unknowns and equations, for the multivariate definitions of Pad approximant under consideration, leads to a Toeplitz-block linear system with coefficient matrix of low displacement rank. Moreover, the matrix is very sparse, especially in higher dimensions. In Section 2 we discuss this for the so-called equation lattice definition and in Section 3 for the homogeneous definition of the multivariate Pad approximant. We do not discuss definitions based on multivariate generalizations of continued fractions [12, 25], or approaches that require some symbolic computations [6, 18]. In Section 4 we present an explicit formula for the factorization of the matrix that results from applying the displacement operator to the Toeplitz-block coefficient matrix. We then generalize the well-known fast Gaussian elimination procedure with partial pivoting developed in [14, 19], to deal with a rectangular block structure where the number and size of the blocks vary. We do not aim for a superfast solver because of the higher risk for instability. Instead we show how the developed technique can be combined with an easy interval arithmetic verification step. Numerical results illustrate the technique in Section 5.