code: def convert_block_tridiag_to_banded(H_diag, H_upper_diag, lower=True):
      
      T, D, _ = H_diag.shape
      assert H_diag.ndim == 3 and H_diag.shape[2] == D
      assert H_upper_diag.shape == (T - 1, D, D)
      H_lower_diag = np.swapaxes(H_upper_diag, -2, -1)
  
      ab = np.zeros((2 * D, T * D))
  
      # Fill in blocks along the diagonal
      for d in range(D):
          # Get indices of (-d)-th diagonal of H_diag
          i = np.arange(d, D)
          j = np.arange(0, D - d)
          h = np.column_stack((H_diag[:, i, j], np.zeros((T, d))))
          ab[d] = h.ravel()
  
      # Fill in lower left corner of blocks below the diagonal
      for d in range(0, D):
          # Get indices of (-d)-th diagonal of H_diag
          i = np.arange(d, D)
          j = np.arange(0, D - d)
          h = np.column_stack((H_lower_diag[:, i, j], np.zeros((T - 1, d))))
          ab[D + d, :D * (T - 1)] = h.ravel()
  
      # Fill in upper corner of blocks below the diagonal
      for d in range(1, D):
          # Get indices of (+d)-th diagonal of H_lower_diag
          i = np.arange(0, D - d)
          j = np.arange(d, D)
          h = np.column_stack((np.zeros((T - 1, d)), H_lower_diag[:, i, j]))
          ab[D - d, :D * (T - 1)] += h.ravel()
  
      return ab if lower else transpose_lower_banded_matrix(ab)
docs: convert blocks to banded matrix representation required for scipy.
      we are using the "lower form."
      see https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solveh_banded.html