Implement a Solve function that can take an initial estimate as input

[tvercaut]

As per the lsmr documentation:

 * Note that x is not an input parameter.
 * If some initial estimate x0 is known and if damp = 0,
 * one could proceed as follows:
 *
 * 1. Compute a residual vector     r0 = b - A*x0.
 * 2. Use LSMR to solve the system  A*dx = r0.
 * 3. Add the correction dx to obtain a final solution x = x0 + dx.

[tvercaut]

Quick and dirty python implementation:

def lsmr_with_init(A,b,x0):
    r0 = b - scipy.sparse.linalg.aslinearoperator(A).matvec(x0)
    deltax_pack = scipy.sparse.linalg.lsmr(A,r0)
    return x0 + deltax_pack[0]

Note however that it makes it more difficult to handle both initialisation and damping.