Unable to set given (local) rows to identity matrix.

Question

I have simple code for solving Navier-Stokes equation in 3D. Unfortunately I'm getting weird error, It looks that newton solver does one iteration and than the error occurs for what ever reason. I tested the code on friend's computer and it outputs the same error(same OS and dolfin version though).

Why I'm getting this error?

---

edit: I was able to get rid of the problem. By changing line

 v, p = w.split()

to line

 v, p = split(w)

But why does this help?

---

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.900e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Traceback (most recent call last):
  File "myStaticNavierStokes.py", line 35, in <module>
    solve( F == 0, w, bcs )
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 297, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 345, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to set given (local) rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.6.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

.

from dolfin import *

mesh = UnitCubeMesh(10,10,10)

bndry = FacetFunction("size_t", mesh)
for f in facets(mesh):
    mp = f.midpoint()
    if near(mp[0], 0.0): # inflow
        bndry[f] = 3
    elif near(mp[0], 1.0): # outflow
        bndry[f] = 2
    elif near(mp[1], 0.0) or near(mp[1], 1.0) or near(mp[2], 0.0) or near(mp[2], 1.0): # walls
        bndry[f] = 1

V = VectorFunctionSpace(mesh, "CG", 2)
P = FunctionSpace(mesh,"CG",1)
W = MixedFunctionSpace([V,P])

noslip = Constant((0,0,0))
vb     = Constant((0,0,1))
bcv_walls = DirichletBC(W.sub(0), noslip, bndry, 1)
bcv_cap   = DirichletBC(W.sub(0), vb, bndry, 3)

bcs = [bcv_cap,bcv_walls]

v_, p_ = TestFunctions(W)
w = Function(W)
v, p = w.split()

Re = Constant(0.1)
I = Identity(mesh.geometry().dim())
T = -p*I + 2.0/Re*sym(grad(v))
F = (inner(T, grad(v_)) - p_*div(v) + inner(grad(v)*v, v_))*dx

solve( F == 0, w, bcs )

plot(w.sub(0),interactive=True)