The finite element analysis of high frequency vibrations of quartz crystal plates is a necessary process required in the design of quartz crystal resonators of precision types for applications in filters and sensors.The anisotropic materials and extremely high frequency in radiofrequency range of resonators determine that vibration frequency spectra are complicated with strong couplings of large number of different vibration modes representing deformations which do not appear in usual structural problems.For instance,the higher-order thickness-shear vibrations usually representing the sharp deformation of thin plates in the thickness direction,expecting the analysis is to be done with refined meshing schemes along the relatively small thickness and consequently the large plane area.To be able to represent the precise vibration mode shapes,a very large number of elements are needed in the finite element analysis with either the three-dimensional theory or the higher-order plate theory,although considerable reduction of numbers of degree-of-freedom(DOF)are expected for the two-dimensional analysis without scarifying the accuracy.We have successfully implemented the Mindlin plate theory for the analysis of quartz crystal resonators with the finite element method based on both linear and nonlinear formulations for vibration frequencies and mode shapes in the vicinity of the thickness-shear frequency of the fundamental and third-order overtone modes.The large number of DOF of the linear system resulted from the finite element formulation has challenged the procedure of eigenvalue computation in a specified interval,again which is not encountered often in structural vibrations.As part of the software development,many different libraries and functions have been utilized in the transformation and eventual evaluation of the eigenvalue problem and the efficiency and accuracy have to be taken into consideration in the solution process.In this paper,we reviewed the software architecture for the analysis and demonstrated the evaluation and tuning of parameters for the improvement of the analysis with problems of elements with a large number of DOF in each node,or a problem with unusually large bandwidth of the banded stiffness and mass matrices in comparison with conventional finite element formulation.Such a problem can be used as an example for the optimization and tuning of problems from multi-physics analysis which are increasingly important in applications with excessive large number of DOF and bandwidth in engineering.