Add diagonalization option to ph.x

Running series of phonon calculations at different geometries of the same crystal is usual in the calculation of thermodynamic properties (e.g. Gibbs free energy) in the quasiharmonic approximation. When many phases and are involved, invariably some of them will fail due to some cdiaghg error or another. This patch implements the "diagonalization" option in ph.x, which allows changing the diagonalization method for the non-scf calculations.