Skip to content

Commit

Permalink
Fix to_boundary for complex numbers (#928)
Browse files Browse the repository at this point in the history
  • Loading branch information
amontoison authored Nov 5, 2024
1 parent 8cdcf6c commit c5174fb
Showing 1 changed file with 6 additions and 6 deletions.
12 changes: 6 additions & 6 deletions src/krylov_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -336,12 +336,12 @@ kcopy!(n :: Integer, y :: AbstractVector{T}, x :: AbstractVector{T}) where T <:

kfill!(x :: AbstractArray{T}, val :: T) where T <: FloatOrComplex = fill!(x, val)

kref!(n, x, y, c, s) = reflect!(x, y, c, s)

kgeqrf!(A :: AbstractMatrix{T}, tau :: AbstractVector{T}) where T <: BLAS.BlasFloat = LAPACK.geqrf!(A, tau)
korgqr!(A :: AbstractMatrix{T}, tau :: AbstractVector{T}) where T <: BLAS.BlasFloat = LAPACK.orgqr!(A, tau)
kormqr!(side :: Char, trans :: Char, A :: AbstractMatrix{T}, tau :: AbstractVector{T}, C :: AbstractMatrix{T}) where T <: BLAS.BlasFloat = LAPACK.ormqr!(side, trans, A, tau, C)

kref!(n, x, y, c, s) = reflect!(x, y, c, s)

macro kswap!(x, y)
quote
local tmp = $(esc(x))
Expand Down Expand Up @@ -377,10 +377,10 @@ function to_boundary(n :: Int, x :: AbstractVector{FC}, d :: AbstractVector{FC},
else
# (dᴴMd) σ² + (xᴴMd + dᴴMx) σ + (xᴴMx - Δ²).
mulorldiv!(z, M, x, ldiv)
rxd = dot(z, d)
xNorm2 = dot(z, x)
rxd = kdot(n, z, d)
xNorm2 = kdotr(n, z, x)
mulorldiv!(z, M, d, ldiv)
dNorm2 = dot(z, d)
dNorm2 = kdotr(n, z, d)
end
dNorm2 == zero(T) && error("zero direction")
flip && (rxd = -rxd)
Expand All @@ -390,7 +390,7 @@ function to_boundary(n :: Int, x :: AbstractVector{FC}, d :: AbstractVector{FC},

# q₂ = ‖d‖², q₁ = xᴴd + dᴴx, q₀ = ‖x‖² - Δ²
# ‖x‖² ≤ Δ² ⟹ (q₁)² - 4 * q₂ * q₀ ≥ 0
roots = roots_quadratic(dNorm2, 2 * rxd, xNorm2 - radius2)
roots = roots_quadratic(dNorm2, 2 * real(rxd), xNorm2 - radius2)

return roots # `σ1` and `σ2`
end
Expand Down

0 comments on commit c5174fb

Please sign in to comment.