Refactor matmul code with mul!

[blegat]

The new implementation catches mul! calls for which the eltype of the result is a JuMP affine or quadratic expressions (resolving #1990).
This allows the generic codes using mul! to benefit from our specialized multiplication methods.
For SparseArrays, the situation is less ideal as some calls are redirected to mul! and some have custom implementation that do not work for JuMP as they start by promoting the arrays to the same eltype...
This PR redirects some of them so that all tests pass but more are needed (more tests and more redirections) but this is left for future work. I don't think there is any regression as we already had this issue: #1276

Closes #1990
Closes #1716

[codecov]

Codecov Report

Merging #1995 into master will increase coverage by 0.13%.
The diff coverage is 82.02%.

@@            Coverage Diff             @@
##           master    #1995      +/-   ##
==========================================
+ Coverage   88.86%   88.99%   +0.13%     
==========================================
  Files          33       33              
  Lines        4282     4306      +24     
==========================================
+ Hits         3805     3832      +27     
+ Misses        477      474       -3

Impacted Files	Coverage Δ
src/quad_expr.jl	89.03% <0%> (-1.1%)	⬇️
src/aff_expr.jl	85.82% <66.66%> (-0.44%)	⬇️
src/operators.jl	86.36% <85.54%> (+3.03%)	⬆️
src/optimizer_interface.jl	80.39% <0%> (+3.46%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 36dc086...6af7e98. Read the comment docs.

[mlubin]

Could you add a test to cover #1990?

Does this close #1716?

[mlubin]

LGTM assuming tests pass.