User defined gradients

This has been come up a few times, and we discussed this yesterday

If a user defined function includes iterative algorithm with a fixed point the implict gradient is 0, but autodiff will diff through that iterative algorithm which may significant cost per each iteration. In such cases it would be useful if the user could define the gradient computation. If that gradient computation could be part of the user defined function (e.g. returning a tuple of value and gradient), then the value and gradient could be computed using the same intermediate values. This would be useful in cases, where the user defined function is something which would not be likely to be added to Stan math, but where avoiding the autodiff through iterative part could make the gradient computation order of magnitude faster.

It's not clear how complicated this would be to implement, but it was suggested that I create this issue and tag @SteveBronder

The issue does not name specific files, tests, or entry points. Start by locating the existing user-defined function and autodiff mechanisms, then determine how a function could provide both its value and gradient while reusing intermediate values. Done means a documented API and tests showing that iterative algorithms avoid unnecessary autodiff work.