Compute correct operator sizes in parallel

In parallel, we cannot look at Function.dof_count to determine the
number of local degrees of freedom, since it counts all dofs (including
those in the halo).  Instead use Function.dof_dset.size where
appropriate.  Furthermore, correct instantiation of PETSc Vecs and Mats
such that we always pass the local size and PETSC_DECIDE for the global
size.  Solver now runs in parallel, but possibly does far too much
work.