|
| 1 | +module MatrixAlgebraKitGenericLinearAlgebraExt |
| 2 | + |
| 3 | +using MatrixAlgebraKit |
| 4 | +using MatrixAlgebraKit: sign_safe, check_input, diagview |
| 5 | +using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! |
| 6 | +using LinearAlgebra: I, Diagonal, lmul! |
| 7 | + |
| 8 | +function MatrixAlgebraKit.default_svd_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} |
| 9 | + return GLA_QRIteration() |
| 10 | +end |
| 11 | + |
| 12 | +for f! in (:svd_compact!, :svd_full!, :svd_vals!) |
| 13 | + @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GLA_QRIteration) = nothing |
| 14 | +end |
| 15 | + |
| 16 | +function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) |
| 17 | + F = svd!(A) |
| 18 | + U, S, Vᴴ = F.U, Diagonal(F.S), F.Vt |
| 19 | + return MatrixAlgebraKit.gaugefix!(svd_compact!, U, S, Vᴴ, size(A)...) |
| 20 | +end |
| 21 | + |
| 22 | +function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, ::GLA_QRIteration) |
| 23 | + F = svd!(A; full = true) |
| 24 | + U, Vᴴ = F.U, F.Vt |
| 25 | + S = MatrixAlgebraKit.zero!(similar(F.S, (size(U, 2), size(Vᴴ, 1)))) |
| 26 | + diagview(S) .= F.S |
| 27 | + return MatrixAlgebraKit.gaugefix!(svd_full!, U, S, Vᴴ, size(A)...) |
| 28 | +end |
| 29 | + |
| 30 | +function MatrixAlgebraKit.svd_vals!(A::AbstractMatrix, S, ::GLA_QRIteration) |
| 31 | + return svdvals!(A) |
| 32 | +end |
| 33 | + |
| 34 | +function MatrixAlgebraKit.default_eigh_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} |
| 35 | + return GLA_QRIteration(; kwargs...) |
| 36 | +end |
| 37 | + |
| 38 | +for f! in (:eigh_full!, :eigh_vals!) |
| 39 | + @eval MatrixAlgebraKit.initialize_output(::typeof($f!), A::AbstractMatrix, ::GLA_QRIteration) = nothing |
| 40 | +end |
| 41 | + |
| 42 | +function MatrixAlgebraKit.eigh_full!(A::AbstractMatrix, DV, ::GLA_QRIteration) |
| 43 | + eigval, eigvec = eigen!(Hermitian(A); sortby = real) |
| 44 | + return Diagonal(eigval::AbstractVector{real(eltype(A))}), eigvec::AbstractMatrix{eltype(A)} |
| 45 | +end |
| 46 | + |
| 47 | +function MatrixAlgebraKit.eigh_vals!(A::AbstractMatrix, D, ::GLA_QRIteration) |
| 48 | + return eigvals!(Hermitian(A); sortby = real) |
| 49 | +end |
| 50 | + |
| 51 | +function MatrixAlgebraKit.default_qr_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} |
| 52 | + return GLA_HouseholderQR(; kwargs...) |
| 53 | +end |
| 54 | + |
| 55 | +function MatrixAlgebraKit.qr_full!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) |
| 56 | + check_input(qr_full!, A, QR, alg) |
| 57 | + Q, R = QR |
| 58 | + return _gla_householder_qr!(A, Q, R; alg.kwargs...) |
| 59 | +end |
| 60 | + |
| 61 | +function MatrixAlgebraKit.qr_compact!(A::AbstractMatrix, QR, alg::GLA_HouseholderQR) |
| 62 | + check_input(qr_compact!, A, QR, alg) |
| 63 | + Q, R = QR |
| 64 | + return _gla_householder_qr!(A, Q, R; alg.kwargs...) |
| 65 | +end |
| 66 | + |
| 67 | +function _gla_householder_qr!(A::AbstractMatrix, Q, R; positive = false, blocksize = 1, pivoted = false) |
| 68 | + pivoted && throw(ArgumentError("Only pivoted = false implemented for GLA_HouseholderQR.")) |
| 69 | + (blocksize == 1) || throw(ArgumentError("Only blocksize = 1 implemented for GLA_HouseholderQR.")) |
| 70 | + |
| 71 | + m, n = size(A) |
| 72 | + k = min(m, n) |
| 73 | + Q̃, R̃ = qr!(A) |
| 74 | + lmul!(Q̃, MatrixAlgebraKit.one!(Q)) |
| 75 | + |
| 76 | + if positive |
| 77 | + @inbounds for j in 1:k |
| 78 | + s = sign_safe(R̃[j, j]) |
| 79 | + @simd for i in 1:m |
| 80 | + Q[i, j] *= s |
| 81 | + end |
| 82 | + end |
| 83 | + end |
| 84 | + |
| 85 | + computeR = length(R) > 0 |
| 86 | + if computeR |
| 87 | + if positive |
| 88 | + @inbounds for j in n:-1:1 |
| 89 | + @simd for i in 1:min(k, j) |
| 90 | + R[i, j] = R̃[i, j] * conj(sign_safe(R̃[i, i])) |
| 91 | + end |
| 92 | + @simd for i in (min(k, j) + 1):size(R, 1) |
| 93 | + R[i, j] = zero(eltype(R)) |
| 94 | + end |
| 95 | + end |
| 96 | + else |
| 97 | + R[1:k, :] .= R̃ |
| 98 | + MatrixAlgebraKit.zero!(@view(R[(k + 1):end, :])) |
| 99 | + end |
| 100 | + end |
| 101 | + return Q, R |
| 102 | +end |
| 103 | + |
| 104 | +function MatrixAlgebraKit.default_lq_algorithm(::Type{T}; kwargs...) where {T <: StridedMatrix{<:Union{BigFloat, Complex{BigFloat}}}} |
| 105 | + return MatrixAlgebraKit.LQViaTransposedQR(GLA_HouseholderQR(; kwargs...)) |
| 106 | +end |
| 107 | + |
| 108 | +end |
0 commit comments