Accept Function not Dat in _la_solve

Additionally, document arguments better.
Add Matrix class

A thin wrapper around a PyOP2 Mat that remembers the bilinear form it
was assembled from.
Use Matrix not Mat for assembled bilinear forms
Apply boundary conditions correctly in _la_solve case

Now that the assembled bilinear form remembers where it came from, we
can correctly apply strong boundary conditions.
Add test of preassembled solve with boundary conditions
Implement lazy bilinear form assembly

When pre-assembling systems that contain boundary conditions it is legal
to do:

    A = assemble(a)
    for bc in bcs: bc.apply(A)
    solve(A, ...)

However, the way we apply boundary conditions to the bilinear form
requires that they are available when we assemble it.  To support the
above, we need to postpone assembling the form until the solve
call (when any boundary conditions must be known).

To do this, when assembling a bilinear form we now return an unassembled
Matrix object, with a assembly callback stashed on it.  When we need it,
we must now call Matrix.assemble().  If the Matrix has already been
assembled, this is a no-op.
Add more tests of preassembled Matrix boundary conditions
Make Matrix.M property force assembly

If we're looking at the op2.Mat through the public property, we probably
want some values in it, so force the assembly of the Matrix.
DTRT if a preassembled Matrix is reused with different bcs

If we preassemble a Matrix and then use it with some set of boundary
conditions, if we later reuse it with different boundary conditions,
force reassembly with the new boundary conditions.  This way, we get the
speed of preassembled things with the DWIMness of just using forms
everywhere.
Add __str__ and __repr__ methods on Matrix
Add more tests of Matrix behaviour

Check that assemble is actually doing things lazily.
Make remaining tests executable as python scripts

If run as python test_foo.py, import pytest and run the file.
Add test that solve with preassembly doesn't modify RHS
Ensure user-provided RHS is not overwritten in la_solve

Previously, we were putting the corrected (BC-applied) RHS into the
user-provided input RHS.  This was bad, because solve(A, u, b) would
modify b unexpectedly.  Fix this, by using the newly assembled Function
instead.
Replace BCs if their subdomain matches an existing one

Rather than adding boundary conditions on the same subdomain to the set
of BCs on the Matrix, if the new boundary condition is on the same
subdomain as an existing one, replace the old one.

Additionally, make checks for when a matrix needs reassembling smarter.
If the boundary conditions have changed, but the subdomains haven't, we
don't need to reassemble.
Add more tests of matrices with boundary conditions