Numerical methods based on transformation of variables are developed to improve the computational efficiency of the variational nodal method (VNM). Reordering and orthogonal transformations of the nodal unknowns are found to reduce the coefficient matrices of VNM into block-diagonal forms. These forms make it possible to reduce greatly the number of floating-point operations in matrix manipulations and hence to reduce the computational times. The red-black response matrix acceleration by transformation of interface partial-current variables has been extended to three-dimensional geometries and higher orders of spatial and angular approximations. These combined methods are incorporated within the algorithms currently used in the variational nodal code VARIANT at Argonne National Laboratory. All primary algorithms ranging from the generation of response matrices to the iterative solution method for the response matrix equations are modified to implement the new formulation. The efficiency of the new methods is tested on eigenvalue problems by comparing the computation times of the new and existing methods. Three-dimensional calculations are performed in hexagonal and Cartesian geometry for various spatial and angular approximations. The test results show that very significant gains can be obtained especially for the coupling coefficient calculations in higher angular approximations. More than an order of magnitude reduction of the total computing time is achieved in the best case.