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.