code: def cholesky(spatial_cov, cross_corr):
      
      M, N = spatial_cov.shape[:2]
      L = numpy.array([numpy.linalg.cholesky(spatial_cov[i]) for i in range(M)])
      LLT = []
      for i in range(M):
          row = [numpy.dot(L[i], L[j].T) * cross_corr[i, j] for j in range(M)]
          for j in range(N):
              singlerow = numpy.zeros(M * N)
              for i in range(M):
                  singlerow[i * N:(i + 1) * N] = row[i][j]
              LLT.append(singlerow)
      return numpy.linalg.cholesky(numpy.array(LLT))
docs: Decompose the spatial covariance and cross correlation matrices.
  
      :param spatial_cov: array of shape (M, N, N)
      :param cross_corr: array of shape (M, M)
      :returns: a triangular matrix of shape (M * N, M * N)