Sampling from a Normal distribution with sparse covariance matrix

Question

To sample from a gaussian distribution with mean zero and covariance matrix S, we can do the following:

from scipy import sparse
from numpy import np
S = sparse.diags([np.full(100,1),0.1*np.ones(99),0.1*np.ones(99)],[0,-1,1])
S.A
array([[1. , 0.1, 0. , ..., 0. , 0. , 0. ],
       [0.1, 1. , 0.1, ..., 0. , 0. , 0. ],
       [0. , 0.1, 1. , ..., 0. , 0. , 0. ],
       ...,
       [0. , 0. , 0. , ..., 1. , 0.1, 0. ],
       [0. , 0. , 0. , ..., 0.1, 1. , 0.1],
       [0. , 0. , 0. , ..., 0. , 0.1, 1. ]])
np.random.multivariate_normal(np.zeros(100), S.A, 1)

Obviously this completely ignores the fact that the covariance matrix is sparse. I am wondering if there is a better way to more efficiently sample from a Gaussian with a sparse matrix without having to use the full 100x100 matrix. In my problem, I am dealing with a 1Mx1M matrix which exhausts my memory, but the matrix is extremely sparse so there should hopefully be a better way to do this, ideally in Python.