Using the help provided here I have written user functions in C++ that are used in my likelihood calculation, and for calculating the derivatives of the likelihood.

I got an example working where the return type was univariate (double in my C++ source code).

I am now trying to implement another example where the return type is VectorXd.

I have managed to get my code to compile, but when I run the executable, there is some kind of run-time error:

```
> ./fhn.exe sample data file=fhn.d
method = sample (Default)
sample
num_samples = 1000 (Default)
num_warmup = 1000 (Default)
...
Aborted (core dumped)
```

Outside of Stan I have written a simple test program that just evaluates my user functions. This seems to work fine.

The next level up from my source code are wrappers that can take `var`

arguments. Here are my wrapper definitions,

```
namespace stan{
namespace math{
template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, -1,1>
spec_fun(const Eigen::Matrix<T0__, -1,1>& p,
const Eigen::Matrix<T1__, -1,1>& theta, std::ostream* pstream__) {
throw std::logic_error("not implemented"); // this should never be called
}
// this is a template specialization for <double, double>
template<>
Eigen::Matrix<double, -1,1> spec_fun<double, double>(const Eigen::Matrix<double, -1,1>& p,
const Eigen::Matrix<double,-1,1>& x,
std::ostream* pstream__) {
return linear_sde::spec_fun(p, x);
}
// this is a template specialization for <double, var>
template<>
Eigen::Matrix<var, -1, 1> spec_fun<double,var>(const Eigen::Matrix<double, -1,1>& p,
const Eigen::Matrix<var, -1, 1>& x,
std::ostream* pstream__) {
vector<var> x_std(x.data(), x.data() + x.rows() * x.cols());
Eigen::VectorXd a = value_of(x);
Eigen::MatrixXd J = linear_sde::spec_fun_jacobian(p, a);
Eigen::Matrix<var, -1, 1> result(a.rows());
std::vector<double> grads(x.rows());
for (int i = 0; i < J.rows(); ++i) {
for (int j = 0; j < a.rows(); ++j)
grads[j] = J(i, j);
result(i) = precomputed_gradients(a(i), x_std, grads);
}
return result;
}
}
}
```

These wrappers were written using the guidance here under the heading of ‘Functions with more than one output’.

Note that in the guidance an `Eigen::VectorXd`

is passed to the function `precomputed_gradients(...)`

. When I tried this I got a compiler error, so I now copy the data into an `std::vector`

before passing it to `precomputed_gradients(...)`

.

My usual way of debugging issues like this is to use the Visual Studio debugger, and this typically gives me a call stack with line numbers. Is it worth trying to do this with a Stan program?

Do you have any other suggestions for how to debug Stan programs? Or what the specific issue with my program might be?