In your code above, you dont need to wrap it in a custom operator, as your operation is reshaping the 3 scaler into a matrix, the easiest way to do is something like:
with pm.Model() as m:
c11 = ...
c12 = ...
c22 = ...
x = tt.zeros((2, 2))
x = tt.inc_subtensor(x[0, 0], c11)
x = tt.inc_subtensor(x[0, 1], c12)
x = tt.inc_subtensor(x[1, 0], c12)
x = tt.inc_subtensor(x[1, 1], c22)
And theano will take care of the gradient for each scaler.
I will reply to you in Multivariate Gaussian with Custom Operator Mean shortly.