For a saddle point problem you can use direct solvers as umfpack or also the build-in sparsecholesky solver if you regularize the problem at the lower-right zero block. E.g. for Stokes involving velocity u and pressure p you can add

For vector space, I can use V=VectorL2(mesh, order = ), then in the following
MatrixL2 = L2(mesh, order=order, dim=4)
u,v = MatrixL2.TnT()
u.dims = (2,2)
v.dims = (2,2)
gf = GridFunction(MatrixL2)
gf.dims = (2,2)

which one is like the above space V written by me? It seems that MatrixL2 = L2(mesh, order=order, dim=4) is different from the space V.