In analogy with the multigrid method, the authors present a solution method for problems based on the \(p\)-version of the finite element method whereby degrees of freedom associated with approximation order \(p\) are visited at once and separately from those of approximation order \(p-1\) and \(p+1\). The authors present the method in an algebraic setting, independently of the finite element formulation from which the system arises. They present three variations on the theme: a straightforward ``V-cycle,'' a nested variation in which lower order levels are visited more frequently than higher order levels, and a combined approach with the nested variation used initially, followed by the V-cycle. Multi-\(p\) methods can be used to precondition the conjugate gradient iterative method. Indeed, the authors remark that an arbitrary preconditioner matrix \(M\) for the conjugate gradient method can be related to a smoothing matrix in a multi-\(p\) method, and the resulting method can be used as a preconditioner for the conjugate gradient method. They go on to demonstrate that this preconditioner reduces the condition number of the iteration and improves its convergence rate. One commercially available computer program employs a particular choice of hierarchical shape functions based on quadrilateral finite elements, solving the resulting system with a direct solver. An alternative approach to solution involves eliminating the so-called internal degrees of freedom using the Schur complement and solving the resulting system by preconditioned conjugate gradients. In this paper, the authors compare the direct solution method, preconditioning by the SSOR method as implemented in ITPACK, and preconditioning by the multi-\(p\) method presented above. Numerical results show multi-\(p\) preconditioning to be superior to the other two methods. The numerical results also indicate that the number of iterations required for solution does not increase as \(p\) increases, in contrast to other available methods.