fix(jacobian): size fwd-mode outputs from the range of f

The fwd jacobian sized both outputs from the domain before f was ever
applied, so the range dimension m was simply assumed equal to n:

    J.resize(x_fvar.size(), x.size());
    fx.resize(x_fvar.size());
    ...
    Matrix<fvar<T>, Dynamic, 1> fx_fvar = f(x_fvar);
    J.col(0) = fx_fvar.d();

For m != n that assigns an m-vector into an n-row column. Move the
J.resize below the call to f and take its rows from fx_fvar.size(),
which mirrors what the rev implementation already does.

The fx.resize is dropped rather than moved: fx = fx_fvar.val() is a
full-matrix assignment that resizes on its own, which is why fx.size()
came out correct even in the aborting runs. Moving it would only have
restated that.

Failure modes before the fix, both directions, with a functor that is
not square:

  assertions on (math's own build never defines NDEBUG) - abort.

  -DNDEBUG, m > n - Eigen's assignment loop bounds on the destination
  block, so the copy truncates: J stays n x n and holds a partial
  Jacobian. Silent wrong answer, no out-of-bounds access.

  -DNDEBUG, m < n - the same destination bound reads past the source
  temporary. ASan reports a heap-buffer-overflow READ at jacobian.hpp:25,
  8 bytes past fx_fvar.

Note that this narrows the issue's description: there is no write past
J's storage in either direction, so no heap corruption.

No signature change, so overload resolution is untouched and the eight
in-tree callers (all in rev/functor, all resolving to the rev overload,
and all square) are unaffected. Verified with ASan clean in both
directions, agreement with the rev overload on the same functor, the
higher-order T = fvar<double> instantiation, the fwd/functor suite,
mix/functor/autodiff_test, and rev/functor/algebra_solver_fp_test.

Closes #3385

Author: SteveBronder

stan/math/fwd/functor/jacobian.hpp

@@ -14,14 +14,14 @@ inline void jacobian(const F& f, const Eigen::Matrix<T, Eigen::Dynamic, 1>& x,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   Matrix<fvar<T>, Dynamic, 1> x_fvar(x.size());
-  J.resize(x_fvar.size(), x.size());
-  fx.resize(x_fvar.size());
   for (int k = 0; k < x.size(); ++k) {
     x_fvar(k) = fvar<T>(x(k), 0);
   }
   x_fvar(0) = fvar<T>(x(0), 1);
   Matrix<fvar<T>, Dynamic, 1> fx_fvar = f(x_fvar);
+  // size the outputs from the range of f, which is only known once f is applied
   fx = fx_fvar.val();
+  J.resize(fx_fvar.size(), x.size());
   J.col(0) = fx_fvar.d();
   const fvar<T> switch_fvar(0, 1);  // flips the tangents on and off
   for (int i = 1; i < x.size(); ++i) {