Walking through the debugger for this code:
%{
aaa1=[4.546541072119, -19.52465262];
aaa2=[4.548657475153, -17.361228432903];
Mew = [3.6604, -5.0249];
Sigma =[3.1531, -4.7301; -4.7301, 7.1776];
mvncdf(aaa1,aaa2,Mew,Sigma)
The line 70 in internal.stats.bvncdf is actually a continuation line. The full line is
p2 = exp(-shk/2 + ...
log(t1.*exp(-bsq./(2*asq)) - t2.*sqrt(2*pi).*Phi(-b./a)));
%}
The argument to the log function is small but negative, presumably rounding error (copied from my command window):
arglog = -2.964393875047479e-323;
The log of a negative has an imaginary component equal to pi
log(arglog)
and the argument to the exp() is (again, copied from my command window)
argexp = -7.438000043200474e+02 + 3.141592653589793e+00i
Normally, I would expect the exp(argexp) to have a small imaginary component, because mathematically, exp(argexp) can be expressed as exp(real(argexp))*exp(1i*imag(argexp))
format long e
exp(real(argexp)), exp(1i*imag(argexp))
But, presumably the product of 1e-324 and 1e-16 rounds down to zero, so we get
exp(real(argexp))*exp(1i*imag(argexp))
whicih is the same as what exp() calculates on its own
exp(argexp)
So, mvncdf doesn't suffer with complex numbers for these double inputs because the imaginary part magically disappears (which isn't to say that bvncdf shouldn't be more careful about the argument to that log, maybe there are some inputs where the imaginary component will still stick around, unless an imaginary piece is dealt with later in the code)
From the error message, it looks like that on GPU code needs to be more careful to account for the possibility of a negative input to log.
When this code is executed
%{
aaa1=z_gridvals(1976,:)-z_gridspacing_down(1976,:);
aaa2=z_gridvals(1976,:)-z_gridspacing_up(1976,:);
%}
are aaa1 and aaa2 ordinary Matlab doubles (i.e., not gpuarray)?
