WIP: Provide linear solver with error estimates

Related issue: #4286 

Currently wraps LAPACK's `xGESVX` routines.