Are you interested in a direct frequency response or a modal frequency response analysis ? If you have several excitation frequencies and large matrices the direct approach is ruled out.
Anyway, you have to solve a complex algebraic system
(-\Omega_i^2 M + 1j \Omega D + K) u = f
for each excitation frequency \Omega_i i = 1,...,n