  The five-dimensional quadrupole collective model based on energy density
functionals (EDF) has often been employed to treat long-range correlations
associated with shape fluctuations in nuclei. Our goal is to derive the
collective inertial functions in the collective Hamiltonian by the local
quasiparticle random phase approximation (QRPA) that correctly takes into
account time-odd mean-field effects. Currently, practical framework to perform
the QRPA calculation with the modern EDFs on the $(\beta,\gamma)$ deformation
space is not available. Toward this goal, we develop an efficient numerical
method to perform the QRPA calculation on the $(\beta,\gamma)$ deformation
space based on the Skyrme EDF. We use the finite amplitude method (FAM) for
efficient calculation of QRPA strength functions for multipole external fields.
We construct a computational code of FAM-QRPA in the three-dimensional
Cartesian coordinate space to handle triaxially deformed superfluid nuclei. We
validate our new code by comparing our results with former QRPA calculations
for axially symmetric nuclei. Isoscalar quadrupole strength functions in
triaxial superfluid nuclei, ${}^{110}$Ru and ${}^{190}$Pt, are obtained within
a reasonable computational cost. QRPA calculations for triaxially deformed
superfluid nuclei based on the Skyrme EDF are achieved with the help of FAM.
This is an important step toward the microscopic calculation of collective
inertial functions of the local QRPA.
