  To help reveal the complete picture of linear kinetic drift modes, four
independent numerical approaches, based on integral equation, Euler initial
value simulation, Euler matrix eigenvalue solution and Lagrangian particle
simulation, respectively, are used to solve the linear gyrokinetic
electrostatic drift modes equation in Z-pinch with slab simplification and in
tokamak with ballooning space coordinate. We identify that these approaches can
yield the same solution with the difference smaller than 1\%, and the
discrepancies mainly come from the numerical convergence, which is the first
detailed benchmark of four independent numerical approaches for gyrokinetic
linear drift modes. Using these approaches, we find that the entropy mode and
interchange mode are on the same branch in Z-pinch, and the entropy mode can
have both electron and ion branches. And, at strong gradient, more than one
eigenstate of the ion temperature gradient mode (ITG) can be unstable and the
most unstable one can be on non-ground eigenstates. The propagation of ITGs
from ion to electron diamagnetic direction at strong gradient is also observed,
which implies that the propagation direction is not a decisive criterion for
the experimental diagnosis of turbulent mode at the edge plasmas.
