  We study the average of the product of the central values of two
$L$-functions of modular forms $f$ and $g$ twisted by Dirichlet characters to a
large prime modulus $q$. As our principal tools, we use spectral theory to
develop bounds on averages of shifted convolution sums with differences ranging
over multiples of $q$, and we use the theory of Deligne and Katz to estimate
certain complete exponential sums in several variables and prove new bounds on
bilinear forms in Kloosterman sums with power savings when both variables are
near the square root of $q$. When at least one of the forms $f$ and $g$ is
non-cuspidal, we obtain an asymptotic formula for the mixed second moment of
twisted $L$-functions with a power saving error term. In particular, when both
are non-cuspidal, this gives a significant improvement on M.~Young's asymptotic
evaluation of the fourth moment of Dirichlet $L$-functions. In the general
case, the asymptotic formula with a power saving is proved under a conjectural
estimate for certain bilinear forms in Kloosterman sums.
