How to compute gradient of Softmax in backward pass of a neural network?

I am trying to implement the backward pass of a Softmax layer.

As an input to my backward-function, I receive the gradient from the next upper layer and I have to pass the calculated gradient of the Softmax layer to the layer "beneath" it.

This is my current code (100% wrong):

def backward(self, error_tensor):
Sz = self.last_forward_value
D = np.zeros(error_tensor.shape)
for i in range(D.shape[0]):
for j in range(D.shape[1]):
if i == j:
D[i, j] = Sz[i, j] * (1 - Sz[i, j])
else:
D[i, j] = -(Sz[i, j] * Sz[i, j])
return D


self.last_forward_value is the softmax 2D-array that was calculated at a previous call to forward.

I have troubles understanding which dimensions the gradient (in my example it is called D) has (or should have). error_tensor is a 2D array (as well as self.last_forward_value) - does this mean the gradient should be 3D?

I tried to google this problem but it seems everyone is only working with 1D input arrays. I saw some articles mentioning that the gradient of a 1D input is 2D (is this correct)?

Updated code (still wrong I guess):

J = np.zeros((self.last_forward_value.shape[0], self.last_forward_value.shape[1], self.last_forward_value.shape[1]))
for k in range(error_tensor.shape[0]):
row = self.last_forward_value[k, :]
for i in range(row.shape[0]):
for j in range(row.shape[0]):
i_val = row[i]
j_val = row[j]
if i == j:
J[k, i, j] = i_val * (1 - i_val)
else:
J[k, i, j] = -i_val * j_val
return np.sum(J, axis=1)