code: def symm_block_tridiag_matmul(H_diag, H_upper_diag, v):
      
      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)
      assert v.shape == (T, D)
  
      out = np.matmul(H_diag, v[:, :, None])[:, :, 0]
      out[:-1] += np.matmul(H_upper_diag, v[1:][:, :, None])[:, :, 0]
      out[1:] += np.matmul(np.swapaxes(H_upper_diag, -2, -1), v[:-1][:, :, None])[:, :, 0]
      return out
docs: Compute matrix-vector product with a symmetric block
      tridiagonal matrix H and vector v.
  
      :param H_diag:          block diagonal terms of H
      :param H_upper_diag:    upper block diagonal terms of H
      :param v:               vector to multiple
      :return:                H * v