diff --git a/project/housQR.jl b/project/housQR.jl new file mode 100644 index 0000000..8a2e45e --- /dev/null +++ b/project/housQR.jl @@ -0,0 +1,138 @@ +module housQR + +import Base: size, show, getproperty, getfield, propertynames, \, * +using LinearAlgebra: norm, I, triu, diagm, dot + +export QRhous, qrfact + +mutable struct QRhous{T <: Real} + A::AbstractVecOrMat{T} + d::AbstractArray{T} + AQ::Union{Nothing, AbstractVecOrMat{T}} + AR::Union{Nothing, AbstractVecOrMat{T}} + + QRhous(A, d, AQ=nothing, AR=nothing) = new{eltype(A)}(A, d, AQ, AR) +end + + + +function householder_vector(x::Vector{T})::Tuple{Vector{T}, T} where T + # returns the normalized vector u such that H*x is a multiple of e_1 + s = norm(x) + if x[1] ≥ 0 + s = -s + end + u = copy(x) + u[1] -= s + u ./= norm(u) + return u, s +end + +function householder_vector(x::Matrix{T})::Tuple{Matrix{T}, T} where T + # returns the normalized vector u such that H*x is a multiple of e_1 + s = norm(x) + if x[1] ≥ 0 + s = -s + end + u = copy(x) + u[1] -= s + u ./= norm(u) + return u, s +end + +function qrfact(A::Matrix{T})::QRhous where T + (m, n) = size(A) + R = deepcopy(A) + d = zeros(n) + + for k ∈ 1:n + (u, s) = householder_vector(R[k:end, k]) + # construct R + d[k] = s + R[k:end, k+1:end] -= 2 * u * (u' * R[k:end, k+1:end]) + R[k:end, k] .= u + end + return QRhous(R, d) +end + + +function qyhous(A::QRhous{T}, y::AbstractArray{T}) where T + m, n = size(A) + z = deepcopy(y) + for j ∈ n:-1:1 + # z[j:m] = z[j:m] - 2 * A.A[j:m, j] * (A.A[j:m, j]' * z[j:m]) + z[j:m] -= 2 * A.A[j:m, j] .* dot(A.A[j:m, j], z[j:m]) + end + return z +end + +function qyhoust(A::QRhous{T}, y::AbstractArray{T}) where T + m, n = size(A.A) + z = deepcopy(y) + for j ∈ 1:n + # z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m]) + z[j:m] -= 2 * A.A[j:m, j] .* dot(A.A[j:m, j], z[j:m]) + end + return z +end + +function calculateQ(A::QRhous{T}) where T + if A.AQ != nothing + return A.AQ + end + m, n = size(A.A) + + A.AQ = zeros(m, 0) + id = Matrix{eltype(A.A)}(I, m, n) + for i ∈ eachcol(id) + A.AQ = [A.AQ qyhous(A, i)] + end + return A.AQ +end + +function calculateR(A::QRhous{T}) where T + if A.AR != nothing + return A.AR + end + m, n = size(A.A) + A.AR = triu(A.A[1:n, :], 1) + diagm(A.d) + return A.AR +end + +function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, A::QRhous{T}) where T + summary(io, A); println(io) + print(io, "Q factor: ") + show(io, mime, A.Q) + print(io, "\nR factor: ") + show(io, mime, A.R) +end + +function Base.getproperty(A::QRhous{T}, d::Symbol) where T + if d === :R + return calculateR(A) + elseif d === :Q + return calculateQ(A) + else + getfield(A, d) + end +end + +Base.propertynames(A::QRhous, private::Bool=false) = (:R, :Q, (private ? fieldnames(typeof(A)) : ())...) + +Base.size(A::QRhous) = size(getfield(A, :A)) + +function (\)(A::QRhous{T}, b::AbstractVector{T}) where T + n, m = size(A) + v = qyhoust(A, b) + x = zeros(m) + for j ∈ m:-1:1 + x[j] = (v[j] - dot(x[j+1:m], A.A[j, j+1:m])) * A.d[j]^-1 + end + return x +end + +function (*)(A::QRhous{T}, x::AbstractVecOrMat{T}) where T + return qyhous(A, (A.R * x)) +end + +end \ No newline at end of file diff --git a/project/testing.ipynb b/project/testing.ipynb index c6e6b41..a7f668f 100644 --- a/project/testing.ipynb +++ b/project/testing.ipynb @@ -11,6 +11,78 @@ "using .thinQR" ] }, + { + "cell_type": "code", + "execution_count": 2, + "id": "8b3d71f8-2f4f-44f0-8c97-dd9d78706163", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x = [-1.5012955407352218e-16, 0.5000000000000001]\n" + ] + }, + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 0.9999999999999991\n", + " 1.9999999999999996\n", + " 2.9999999999999996" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = [1. 2; 3. 4.; 5. 6.]\n", + "A = qrhous!(A)\n", + "b = [1., 2, 3]\n", + "x = A \\ b\n", + "@show x\n", + "A * x" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c7f143dc-a80d-4326-83c1-31dac3a6e3e0", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x = [-1.0000000000000033, 1.0000000000000002, 4.0792198665315527e-16]\n" + ] + }, + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 1.0000000000000009\n", + " 2.000000000000001\n", + " 3.000000000000001" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = [1. 2. 3.; 1. 3. 6.; 1. 4. 5.]\n", + "A = qrhous!(A)\n", + "b = [1., 2, 3]\n", + "x = A \\ b\n", + "@show x\n", + "A * x" + ] + }, { "cell_type": "code", "execution_count": 2, @@ -49,58 +121,1560 @@ }, { "cell_type": "code", - "execution_count": 3, - "id": "8b3d71f8-2f4f-44f0-8c97-dd9d78706163", + "execution_count": 7, + "id": "d7d24c53-5485-4f8a-bea3-3fe4ec1f456b", + "metadata": {}, + "outputs": [], + "source": [ + "using BenchmarkTools\n", + "A=randn(10000, 100)\n", + "a=deepcopy(A)\n", + "y=randn(10000);" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "652fe3c4-6b44-45e1-9bdc-5d7d41b55c6a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "BenchmarkTools.Trial: 11 samples with 1 evaluation.\n", + " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m458.557 ms\u001b[22m\u001b[39m … \u001b[35m496.613 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m18.11% … 20.40%\n", + " Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m475.247 ms \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m20.13%\n", + " Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m475.772 ms\u001b[22m\u001b[39m ± \u001b[32m 8.688 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m19.88% ± 0.73%\n", + "\n", + " \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m█\u001b[39m \u001b[34m \u001b[39m\u001b[39m▃\u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n", + " \u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m▇\u001b[34m▇\u001b[39m\u001b[39m█\u001b[32m▁\u001b[39m\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m \u001b[39m▁\n", + " 459 ms\u001b[90m Histogram: frequency by time\u001b[39m 497 ms \u001b[0m\u001b[1m<\u001b[22m\n", + "\n", + " Memory estimate\u001b[90m: \u001b[39m\u001b[33m1.54 GiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m4804\u001b[39m." + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = deepcopy(a)\n", + "@benchmark begin x = qrhous!($A) \\ $y; end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c8e93f5c-b32c-4334-bd97-9779360d5d81", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "BenchmarkTools.Trial: 77 samples with 1 evaluation.\n", + " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m56.414 ms\u001b[22m\u001b[39m … \u001b[35m98.485 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 0.00%\n", + " Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m62.403 ms \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m0.00%\n", + " Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m65.010 ms\u001b[22m\u001b[39m ± \u001b[32m 8.780 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m0.75% ± 1.84%\n", + "\n", + " \u001b[39m \u001b[39m█\u001b[39m▃\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[34m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m▄\u001b[39m▂\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \n", + " \u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m▃\u001b[39m▃\u001b[39m▅\u001b[39m▅\u001b[39m▃\u001b[39m▃\u001b[39m▁\u001b[39m▃\u001b[34m▃\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m▃\u001b[39m▁\u001b[39m█\u001b[39m█\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▄\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▅\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▄\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▃\u001b[39m▃\u001b[39m \u001b[39m▁\n", + " 56.4 ms\u001b[90m Histogram: frequency by time\u001b[39m 88 ms \u001b[0m\u001b[1m<\u001b[22m\n", + "\n", + " Memory estimate\u001b[90m: \u001b[39m\u001b[33m7.92 MiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m624\u001b[39m." + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using LinearAlgebra\n", + "\n", + "@benchmark begin x = $a \\ $y; end" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "43d521a7-6bf4-45c6-abb6-0a07a666570d", + "metadata": {}, + "outputs": [], + "source": [ + "using Profile, ProfileSVG" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "9443e1c4-03dd-46cb-a8ea-652c6862dbcf", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n" + ], + "text/html": [ + "\n", + "\n", + "
\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "ProfileSVG.FGConfig(Node(FlameGraphs.NodeData(ip:0x0, 0x01, 1:382)), Dict{Symbol, Any}(), FlameGraphs.FlameColors(ColorTypes.RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(0.882,0.698,1.0), RGB{N0f8}(0.435,0.863,0.569), RGB{N0f8}(0.0,0.71,0.545), RGB{N0f8}(0.173,0.639,1.0)], RGB{N0f8}(1.0,1.0,1.0), RGB{N0f8}(0.0,0.0,0.0), ColorTypes.RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(0.953,0.0,0.302), RGB{N0f8}(0.894,0.0,0.255), RGB{N0f8}(0.831,0.129,0.216), RGB{N0f8}(0.773,0.192,0.184)], ColorTypes.RGB{FixedPointNumbers.N0f8}[RGB{N0f8}(1.0,0.627,0.0), RGB{N0f8}(1.0,0.643,0.0), RGB{N0f8}(0.965,0.651,0.039), RGB{N0f8}(0.894,0.655,0.11)]), :fcolor, :fcolor, 1.0, false, 50, 2000, 960.0, 0.0, 2.0, \"inherit\", 12.0, false, :none, 0.001)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = a\n", + "@profview qrhous!(A) \\ y" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "fb2484bd-fd33-4b42-bc14-1006c4d9811d", + "metadata": {}, + "outputs": [], + "source": [ + "include(\"housQR.jl\")\n", + "using .housQR" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "eb827167-7fa3-4106-9bac-71a3c237ad5b", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", - " -1.5012955407352218e-16\n", - " 0.5000000000000001" + " 5.810770900369532e-16\n", + " 0.4999999999999994" ] }, - "execution_count": 3, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "A = [1. 2; 3. 4.; 5. 6.]\n", - "A = qrhous!(A)\n", - "b = [1., 2, 3]\n", - "x = A \\ b" + "A = [1. 2; 3 4; 5 6]\n", + "qr = qrfact(A)\n", + "\n", + "A \\ [1., 2, 3]" ] }, { "cell_type": "code", - "execution_count": 4, - "id": "1883485c-b026-4a8a-8c34-b7c74ad79868", + "execution_count": 9, + "id": "41d9dca0-6880-4133-8dc3-c35d5c0c9a49", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "3-element Vector{Float64}:\n", - " 0.9999999999999991\n", - " 1.9999999999999996\n", - " 2.9999999999999996" + "BenchmarkTools.Trial: 21 samples with 1 evaluation.\n", + " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m): \u001b[39m\u001b[36m\u001b[1m222.384 ms\u001b[22m\u001b[39m … \u001b[35m274.339 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m14.47% … 24.11%\n", + " Time \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m): \u001b[39m\u001b[34m\u001b[1m250.843 ms \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m): \u001b[39m18.85%\n", + " Time \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m): \u001b[39m\u001b[32m\u001b[1m245.919 ms\u001b[22m\u001b[39m ± \u001b[32m 16.588 ms\u001b[39m \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m): \u001b[39m19.02% ± 2.58%\n", + "\n", + " \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▃\u001b[39m▃\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[32m \u001b[39m\u001b[39m \u001b[39m \u001b[39m \u001b[34m \u001b[39m\u001b[39m█\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▃\u001b[39m \u001b[39m \n", + " \u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[32m▁\u001b[39m\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[34m▇\u001b[39m\u001b[39m█\u001b[39m▇\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▇\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m█\u001b[39m \u001b[39m▁\n", + " 222 ms\u001b[90m Histogram: frequency by time\u001b[39m 274 ms \u001b[0m\u001b[1m<\u001b[22m\n", + "\n", + " Memory estimate\u001b[90m: \u001b[39m\u001b[33m1.54 GiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m4617\u001b[39m." ] }, - "execution_count": 4, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "A * x" + "@benchmark qrfact($A) \\ $y" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f505daf2-0a23-415a-bfb7-b282acee8294", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Julia 1.9.3", + "display_name": "Julia 1.9.4", "language": "julia", "name": "julia-1.9" }, @@ -108,7 +1682,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.9.3" + "version": "1.9.4" } }, "nbformat": 4, diff --git a/project/thinQR.jl b/project/thinQR.jl index cf4d561..114b9d6 100644 --- a/project/thinQR.jl +++ b/project/thinQR.jl @@ -16,7 +16,6 @@ end function qrhous!(A::AbstractMatrix{T})::QRthin{T} where T m, n = size(A) - d = zeros(n) @inbounds begin for j ∈ 1:n @@ -25,7 +24,7 @@ function qrhous!(A::AbstractMatrix{T})::QRthin{T} where T d[j]= copysign(s, -A[j,j]) fak = sqrt(s * (s + abs(A[j,j]))) A[j,j] -= d[j] - + A[j:m, j] ./= fak if j < n @@ -44,7 +43,8 @@ function qyhous(A::QRthin{T}, y::AbstractArray{T}) where T m, n = size(A.A) z = deepcopy(y) for j ∈ n:-1:1 - z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m]) + # z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m]) + z[j:m] -= A.A[j:m, j] .* dot(A.A[j:m, j], z[j:m]) end return z end