Is this expected?

```
julia> sqrt([0 1; 0 0])
2Γ2 Matrix{Float64}:
0.0 0.0
0.0 0.0
```

I thought it would give an error.

Is this expected?

```
julia> sqrt([0 1; 0 0])
2Γ2 Matrix{Float64}:
0.0 0.0
0.0 0.0
```

I thought it would give an error.

Throwing an error here would be hard. The `sqrt`

is defined as long as the eigenvalues are unique, and numerically determining this would be really unreliable.

WolframAlpha produces the same result happily and quietly:

```
MatrixPower[{{0, 1}, {0, 0}}, 1/2]
#result:
{{0, 0}, {0, 0}}
```

Itβs fine to have repeated eigenvalues, even for a defective matrix:

```
julia> A = [4 1; 0 4]
2Γ2 Matrix{Int64}:
4 1
0 4
julia> S = sqrt(A)
2Γ2 Matrix{Float64}:
2.0 0.25
0.0 2.0
julia> S^2
2Γ2 Matrix{Float64}:
4.0 1.0
0.0 4.0
```

The problem only arises for repeated eigenvalues that are exactly zero and defective, I think?

It seems like we could check for this case in `_sqrt_quasitriu_diag_block!`

β perhaps you could just set a flag on this line when `iszero(A[i,i])`

, and throw an error if you hit that condition twice in a row (assuming the diagonal entries are sorted, as usual for Schur).

Of course, roundoff errors may defeat this too if the matrix is not in Schur-like form to start with:

```
julia> A = [0 1; 0 0]
2Γ2 Matrix{Int64}:
0 1
0 0
julia> sqrt(A) # fail
2Γ2 Matrix{Float64}:
0.0 0.0
0.0 0.0
julia> B = rand(2,2); # random change of basis
julia> A2 = B * A / B # similar to A, should also not have a square root
2Γ2 Matrix{Float64}:
0.68983 -0.154736
3.07534 -0.68983
julia> S = B \ sqrt(A2) * B # change basis from βA2 back to original basis
2Γ2 Matrix{ComplexF64}:
7.63631e-5+1.68025e-20im 6547.66+5.55286e-13im
-1.34717e-12-3.78013e-21im 7.63631e-5-4.38691e-13im
julia> S^2 # should == A if square root existed
2Γ2 Matrix{ComplexF64}:
-2.9895e-9-2.4751e-17im 1.0-2.8724e-9im
-2.05748e-16+1.36679e-26im -2.9895e-9-9.17507e-17im
julia> norm(S^2 - A) / norm(A) # loses half of the significant digits due to ill-conditioning
3.147765149816803e-8
```

3 Likes