Inconsistency in treatment of qz_criterium in dyn_first_order_solver.m / mjdgges.c

I was under the impression that qz_criterium was intended to be a threshold on the absolute value of dr.eigval. See e.g. line 85 of AIM_first_order_solver.m:

nba = nd-sum( abs(dr.eigval) < qz_criterium );

However, it appears that in dyn_first_order_solver.m it is instead a threshold on the squared absolute value of dr.eigval. See line 190 of dyn_first_order_solver.m:

[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium, DynareOptions.qz_zero_threshold);

line 31 of mjdgges.c:

return ((*alphar **alphar + *alphai **alphai) < criterium **beta **beta);

and line 107 of mjdgges.c:

criterium = *mxGetPr(prhs[2]);

The eigenvalues are ( alphar[j] + alphai[j] * 1i ) / beta[j] according to the MKL manual. Thus, it seems that this code is comparing the squared absolute value of the eigenvalue to qz_criterium.

The easiest fix would be to replace line 190 of dyn_first_order_solver.m with:

[err, ss, tt, w, sdim, dr.eigval, info1] = mjdgges(E, D, DynareOptions.qz_criterium .^ 2, DynareOptions.qz_zero_threshold);

so you're comparing squared to squared.