spurious NaN values from matrix square root

[stevengj]

As reported on discourse:

A = [
0.9999999999999998	4.649058915617843e-16	-1.3149405273715513e-16	9.9959579317056e-17
-8.326672684688674e-16	1.0000000000000004	2.9280733590254494e-16	-2.9993900031619594e-16
9.43689570931383e-16	-1.339206523454095e-15	1.0000000000000007	-8.550505126287743e-16
-6.245004513516506e-16	-2.0122792321330962e-16	1.183061278035052e-16	1.0000000000000002
]

B = [
0.09648289218436859	0.023497875751503007	0.0	0.0
0.023497875751503007	0.045787575150300804	0.0	0.0
0.0	0.0	0.0	0.0
0.0	0.0	0.0	0.0
]
sqrt(A*B*A')

gives

4×4 Array{Complex{Float64},2}:
 NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
 NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
 NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
 NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im

This is incorrect. In fact, up to roundoff error, the correct answer is:

4×4 Array{Float64,2}:
 0.307265   0.0455076  0.0  0.0
 0.0455076  0.209085   0.0  0.0
 0.0        0.0        0.0  0.0
 0.0        0.0        0.0  0.0

Roundoff errors may give a complex result here (see also JuliaLang/julia#35057), but should not lead to NaN values.

[andreasnoack]

I think we'd have to handle exact zeros in the diagonal in https://github.com/JuliaLang/julia/blob/c7f1ffbaf0a67983dadca1855a287f58425e3c55/stdlib/LinearAlgebra/src/triangular.jl#L2580-L2597. We end up with a call like sylvester(0, 0, ε) which would have to return NaN. Maybe we should just set the strict triangular elements above zero diagonal elements to zero.