The calculation of dynamic response of multiple-degree-of-freedom viscoelastic linear systems is considered. Viscoelastic forces depend on the past history of motion via convolution integrals over exponentially decaying kernel functions. Exact closed-form expressions for the dynamic response due to general forces and initial conditions are derived in terms of the eigensolutions of the system in the original space. Eigensolutions of the viscoelastic system in turn are obtained approximately as functions of the elastic eigensolutions. This enables one to approximately calculate the dynamic response of complex viscoelastic systems by simple post-processing of the elastic (undamped) eigensolutions. Suitable examples are given to illustrate the derived results.