Constraint multipliers

I changed solve() so that it returns the Lagrange multipliers associated with the constraint functions (previously, only the multipliers associated with the upper and lower bounds on variables were returned.) 

Because this changes the number of values returned by solve(), it will break existing code. If you want to avoid that I could rewrite this to make the return of the constraint multipliers optional.