The classical Pade theory is viewed in an abstract algebraic framework with respect to power series over an arbitrary integral domain. A symbolic manipulation algorithm is developed to compute Pade approximants for power series with polynomial coefficients. The algorithm is based on a new fraction-free elimination algorithm for symmetric indefinite systems of linear equations and it exploits the block structure of the Pade table in case of singularity. The algorithm is presented in a precise algorithmic notation that is suitable for translation into any of the languages for symbolic computation.