I have set up a linear system using a BlockMatrix and BlockVector, where blocks arise from assembled systems and self introduced linear systems. But I see that the solve doesn't take these entities as valid inputs. How can solve such a system which has been constructed using these objects? Also will it be possible to separate out parts which correspond to different blocks (something akin to the split() method used for MixedFunctionSpaces ) ?Below is a minimal code for what I seek to do inside an iteration eventually. I'd appreciate any comments regarding the quality of the code so far: please suggest improvements if you think some things can be done more elegantly. For now I am struggling to get the code working but this input will also go a long way.from dolfin import *from mshr import *from petsc4py import PETScN     = 128alpha = 1e-3mesh = UnitSquareMesh(N,N)V = FunctionSpace(mesh, 'CG', 1)u = TrialFunction(V)v = TestFunction(V)yd = Expression('exp(-50*(pow(x[0]-0.5,2) + pow(x[1]-0.5,2)))')a  = inner(grad(u),grad(v))*dxbc = DirichletBC(V, Constant(0.0), on_boundary)z  = yd*v*dx      A, z  = assemble_system(a,z,bc)            AT = assemble(adjoint(a)) bc.apply(AT)     M = assemble(u*v*dx)       bc.apply(M)n = V.dim()                Ap = Am = zAp.zero()Am.zero()Ap_old = ApAm_old = Amgamma = 1As = as_backend_type(A)As.zero()dvec  = (Ap+Am)dvec*=(gamma)As.set_diagonal(dvec)H  = BlockMatrix(2,2)H[0, 0] = AH[1, 0] = AsH[0, 1] = -1*MH[1, 1] = ATg = BlockVector(2)g[0] = alpha*gamma*(Ap-Am)g[1] = -1*zx = BlockVector(2)solve(H, x, g)  # Evidently, this is where it throws an error(y,p) = split(x)  
Fenics: Solving Linear System with BlockMatrix and BlockVector
petsc;fenics