BlockSeries produced of wrong shape

Running block_diagonalize in a 1x1 AA subspace produces BlockSeries of wrong shapes. For example,

import numpy as np
from pymablock import block_diagonalize

N = 10
H_0 = np.diag(np.arange(N))
H_1 = np.random.rand(N, N) + 1j * np.random.rand(N, N)
H_1 = H_1 + H_1.T.conj()

# First, a series only in H_1
N_A = 1
subspace_indices = N_A * [0] + (N - N_A) * [1]
_, U, _ = block_diagonalize([H_0, H_1], subspace_indices=subspace_indices)

print(U[0, 0, 3].shape)
print(U[1, 1, 3].shape)
print(U[0, 1, 3].shape)

prints (9, 9) for all three blocks of U. These should be (1, 1), (9, 9) and (1, 9), respectively. The code works correctly if N_A > 1. The code works correctly with sympy matrices instead of numpy arrays, so this probably has to do with numpy slicing of a 1-sized array. We should write a test for it too.