Have matrix square root account for Rii=Rjj=0 to prevent NaN values (#…

…35758)
JuliaLang · May 8, 2020 · 880c731 · 880c731
1 parent ba91a98
commit 880c731
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 1 deletion.
diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl
@@ -2607,7 +2607,9 @@ function sqrt(A::UpperTriangular{T},::Val{realmatrix}) where {T,realmatrix}
             @simd for k = i+1:j-1
                 r -= R[i,k]*R[k,j]
             end
-            iszero(r) || (R[i,j] = sylvester(R[i,i],R[j,j],-r))
+            if !(iszero(r) || (iszero(R[i,i]) && iszero(R[j,j])))
+                R[i,j] = sylvester(R[i,i],R[j,j],-r)
+            end
         end
     end
     return UpperTriangular(R)

diff --git a/stdlib/LinearAlgebra/test/triangular.jl b/stdlib/LinearAlgebra/test/triangular.jl
@@ -651,4 +651,9 @@ end
     @test_throws ArgumentError LinearAlgebra.powm(A, 2.2)
 end
 
+# Issue 35058
+let 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]
+    @test sqrt(A*B*A')^2 ≈ A*B*A'
+end
 end # module TestTriangular