I want to generate correlated random variables with a given correlation matrix, means, and variances. Does the Cholesky decomposition only work when the initial random variables are iids with the same mean and variance?
Let Z be uncorrelated random variables normally distributed with mean 0 and variance 1. This means $$Z\sim N\left(0,I\right)$$ If you make an affine transformation $$X\equiv A+BZ$$ then $X$ has a distribution $$X\sim N\left(A,BB'\right)$$ In our case, we want $BB'=\Sigma$, so applying the cholesky decomposition to $\Sigma$ is one way to find a suitable $B$. Thus, to simulate from $X\sim N\left(\mu,\Sigma\right)$, you would simulate $Z$, set $A=\mu$ and $B=chol(\Sigma)$, and apply the transformation above.
So to answer your question, uncorrelated variables of mean 0 and variance 1 can be transformed to generic multivariate normal distributions through the use of affine transformations, depending on the mean vector and cholesky decomposition of the covariance matrix. This is how multivariate normal random number generators generally work.
To address your other point in the comments, suppose you simulate $Z$ and make a transformation $$Y\equiv chol(C)Z$$ so that $Y\sim N(0, C)$. You can still get to $X$ from $Y$. As I discussed in the comments, you can easily do this by multiplying the univariate distributions by the respective standard deviations and adding the respective mean.
More formally, you can consider this another affine transformation $$X\equiv A+SY$$ However, since $Y$ is not uncorrelated, the affine transformation would give $$X\sim N\left(A,SCS'\right)$$ In this case, to get $SCS'=\Sigma$, you would need to set $S$ equal to a matrix with the standard deviations on the diagonal and zeros elsewhere. This is equivalent to applying the transformations on a univariate basis.
We need to decompose the convariance matrix Σ into orthogonal matrix product as mentionned above.
There are different ways to get the orthogonal matrix B, Cholesky is one of them. We can also use SVD (Singular Value Decomposition) to get B, equivalent to calculate the eigenvalues/vectors for PDM.
If you want to use the decomposition to simulate brownians over time, you can also use Brownian Bridge with Haar transformation.