Solving Linear System with BlockMatrix and BlockVector

Question

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 PETSc

N     = 128
alpha = 1e-3
mesh = 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))*dx
bc = 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 = z
Ap.zero()
Am.zero()
Ap_old = Ap
Am_old = Am

gamma = 1

As = as_backend_type(A)
As.zero()
dvec  = (Ap+Am)
dvec*=(gamma)
As.set_diagonal(dvec)

H  = BlockMatrix(2,2)
H[0, 0] = A
H[1, 0] = As
H[0, 1] = -1*M
H[1, 1] = AT

g = BlockVector(2)
g[0] = alpha*gamma*(Ap-Am)
g[1] = -1*z

x = BlockVector(2)
solve(H, x, g)  # Evidently, this is where it throws an error
(y,p) = split(x)  

Answer 1

Consider using cbc.block for this:
https://bitbucket.org/fenics-apps/cbc.block

Comments

I am not sure if there is a confusion. cbc.block seems to be more targeted towards a block preconditioning. In my sample code, I attempt to combine two linear systems, out of which one arises from a variational form, the other is by construction in a single Hx=g form giving a blocked structure for H. Neither of BlockMatrix and block_matrix objects are conducive to be used inside a solve() command which desires a GenericLinearOperator object for H. I am looking for ways to get such an assembled block system as a GenericLinearOperator object or to that of its derived classes.

---

Consider using the backend PETSc to handle the solve rutine and
avoid using the GenericLinearOperator. PETSc has tools for block
matrices and block vectors. I don't think the BlockMatrix <-> GenericLinearOperator
has been much used and it might take some work.