  We present a new direct logarithmically optimal in theory and fast in
practice algorithm to implement the high order finite element method on
multi-dimensional rectangular parallelepipeds for solving PDEs of the Poisson
kind. The key points are the fast direct and inverse FFT-based algorithms for
decomposition in eigenvectors of the 1D eigenvalue problems for the high order
FEM. The algorithm can further be used for numerous applications, in
particular, to implement the high order finite element methods for various
time-dependent PDEs.
