diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl index 00f9cbbd..8d224503 100644 --- a/src/implementations/polar.jl +++ b/src/implementations/polar.jl @@ -123,6 +123,7 @@ function _left_polarnewton!(A::AbstractMatrix, W, P = similar(A, (0, 0)); tol = Rc .= R Rᴴinv = ldiv!(UpperTriangular(Rc)', one!(Rᴴinv)) else # m == n + Q = nothing R = A Rc = view(W, 1:n, 1:n) Rc .= R @@ -163,6 +164,7 @@ function _right_polarnewton!(A::AbstractMatrix, Wᴴ, P = similar(A, (0, 0)); to copy!(Lc, L) Lᴴinv = ldiv!(LowerTriangular(Lc)', one!(Lᴴinv)) else # m == n + Q = nothing L = A Lc = view(Wᴴ, 1:m, 1:m) Lc .= L diff --git a/src/implementations/qr.jl b/src/implementations/qr.jl index f78d8c44..6b1cc58b 100644 --- a/src/implementations/qr.jl +++ b/src/implementations/qr.jl @@ -141,6 +141,8 @@ function householder_qr!( (inplaceQ && (computeR || positive || blocksize > 1 || m < n)) && throw(ArgumentError("inplace Q only supported if matrix is tall (`m >= n`), R is not required, and using the unblocked algorithm (`blocksize = 1`) with `positive = false`")) + jpvt = Vector{Int}(undef, 0) + τ = Vector{eltype(A)}(undef, 0) # Compute QR in packed form if blocksize > 1 nb = min(minmn, blocksize) diff --git a/src/pullbacks/eig.jl b/src/pullbacks/eig.jl index a03eb3c4..0cb90cb8 100644 --- a/src/pullbacks/eig.jl +++ b/src/pullbacks/eig.jl @@ -135,15 +135,12 @@ function eig_trunc_pullback!( (n, n) == size(ΔA) || throw(DimensionMismatch()) G = V' * V + VᴴΔV = !iszerotangent(ΔV) ? V' * ΔV : zero(G) + ΔVperp = ΔV - V * inv(G) * VᴴΔV if !iszerotangent(ΔV) (n, p) == size(ΔV) || throw(DimensionMismatch()) - VᴴΔV = V' * ΔV check_eig_cotangents(D, VᴴΔV; degeneracy_atol, gauge_atol) - - ΔVperp = ΔV - V * inv(G) * VᴴΔV VᴴΔV .*= conj.(inv_safe.(transpose(D) .- D, degeneracy_atol)) - else - VᴴΔV = zero(G) end if !iszerotangent(ΔDmat) diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl index 01fdc4f7..fe20608b 100644 --- a/src/pullbacks/svd.jl +++ b/src/pullbacks/svd.jl @@ -47,6 +47,8 @@ function svd_pullback!( Ur = view(U, :, 1:r) Vᴴr = view(Vᴴ, 1:r, :) Sr = view(S, 1:r) + indU = axes(U, 2) + indV = axes(Vᴴ, 1) # Extract and check the cotangents ΔU, ΔSmat, ΔVᴴ = ΔUSVᴴ diff --git a/src/yalapack.jl b/src/yalapack.jl index 576fe3c5..78dd81a8 100644 --- a/src/yalapack.jl +++ b/src/yalapack.jl @@ -2033,6 +2033,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in for i in 1:2 # first call returns lwork as work[1] #! format: off if eltype(A) <: Complex + rwork_ = isnothing(rwork) ? Vector{$relty}(undef, 0) : rwork ccall((@blasfunc($gesvd), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -2040,7 +2041,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in Ptr{BlasInt}, Clong, Clong), jobu, jobvt, m, n, A, lda, S, U, ldu, Vᴴ, ldv, - work, lwork, rwork, + work, lwork, rwork_, info, 1, 1) else ccall((@blasfunc($gesvd), libblastrampoline), Cvoid, @@ -2135,6 +2136,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in for i in 1:2 # first call returns lwork as work[1] #! format: off if eltype(A) <: Complex + rwork_ = isnothing(rwork) ? Vector{$relty}(undef, 0) : rwork ccall((@blasfunc($gesdd), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -2142,7 +2144,7 @@ for (gesvd, gesdd, gesvdx, gejsv, gesvj, elty, relty) in Ptr{BlasInt}, Clong), job, m, n, A, lda, S, U, ldu, Vᴴ, ldv, - work, lwork, rwork, iwork, + work, lwork, rwork_, iwork, info, 1) else ccall((@blasfunc($gesdd), libblastrampoline), Cvoid,