I would like to generate $S$ samples of a $T \cdot M$ dimensional vector, where $T$ is the number of time steps and $M$ the number of locations, i.e., the vector is a stack with $T$ values for location 1, then $T$ values for location 2, and so on. For small $T$ and $M$, it is convenient to simulate the $S$ samples using a $T \cdot M \times T \cdot M$ covariance matrix and Cholesky decomposition.
Say that I would like to generate $T = 30\cdot24 = 720$ hourly timesteps, i.e., a full month. Already for small M, the above methodology becomes computationally infeasible as $720 \cdot M$ becomes large. However, the temporal correlation in my data is weak and it can be assumed to be zero after 12 hours, for example. So, in practice, I could simulate 30 independent instances of a vector of the length $24 \cdot M$, i.e. dividing the month to days. However, to obtain a valid sample of $720 \cdot M$, I cannot simply join the 30 vectors of length $24 \cdot M$ as there will be discontinuities between the days.
So, my question is that how to make the above simulation computationally feasible? Is it possible to implement some kind of conditional sequential sampling, where a new day would be conditioned on the previous day so that the temporal correlation would be preserved around the change of the day? If yes, what is this technique called and how could it be implemented in Matlab or R?