A linear algebraic theory called the 'multiple Arnoldi method' is presented and realizes large-scale (order-N) electronic structure calculations with generalized eigenvalue equations. A set of linear equations, in the form of (zS - H)x = b, are solved simultaneously with multiple Krylov subspaces. The method is implemented in a simulation package ELSES (www.elses.jp) with tight-binding-form Hamiltonians. A finite-temperature molecular dynamics simulation is carried out for metallic and insulating materials. A calculation with 10(7) atoms was realized by a workstation. The parallel efficiency is shown up to 1024 CPU cores.