A perturbation method based on the variational nodal method for solving the neutron transport equation is developed for multidimensional geometries. The method utilizes the solution of the corresponding adjoint transport equation to calculate changes in the critical eigenvalue due to crosssection changes. Both first-order and exact perturbation theory expressions are derived. The adjoint solution algorithm has been formulated and incorporated into the variational nodal option of the Argonne National Laboratory DIF3D production code. To demonstrate the efficacy of the methods, perturbation calculations are performed on the three-dimensional Takeda benchmark problems in both Cartesian and hexagonal geometries. The resulting changes in eigenvalue are also obtained by direct calculation with the variational nodal method and compared with the change approximated by the first-order and exact theory expressions from the perturbation method. Exact perturbation results are in excellent agreement with the actual eigenvalue differences calculated in VARIANT. First-order theory holds well for sufficiently small perturbations. The times required for the perturbation calculations are small compared with those expended for the base-forward and adjoint calculations.