Making two correlated normals from independent ones

You have a routine that outputs independent standard normal draws. Using it, how would you construct a pair of $N(0,1)$ random variables whose correlation equals a given $\rho \in [-1,1]$?

Take $X=Z_1$ and $Y=\rho Z_1+\sqrt{1-\rho^2}\,Z_2$ with $Z_1,Z_2$ independent standard normals. Then $\operatorname{Var}(Y)=\rho^2+(1-\rho^2)=1$ and $\operatorname{Cov}(X,Y)=\rho$. This is the $2\times 2$ Cholesky factor of the correlation matrix; the same idea extends to $n$ variables via $X=LZ$ with $LL^\top=\Sigma$.