A parallel finite element implementation on tetrahedral grids of the nonlinear three-dimensional full-Stokes model for the dynamics of ice-sheets is presented. Effective solvers by using preconditioning techniques for the saddle-point system resulting from the discretization are discussed and implemented. We validated our finite element full-Stokes model through the use of well-known ice-sheet benchmark experiments, and the solver is shown to be efficient, robust, and scalable.