diff --git a/project/testing.ipynb b/project/testing.ipynb new file mode 100644 index 0000000..ac57874 --- /dev/null +++ b/project/testing.ipynb @@ -0,0 +1,83 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "57a2b335-a7f3-4ace-a2d8-52219c4febc5", + "metadata": {}, + "outputs": [], + "source": [ + "include(\"thinQR.jl\")\n", + "using .thinQR" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "426735f0-859a-4f7d-83cf-f52f464a2737", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 2.8498844747360086\n", + " -2.404529848146403\n", + " -0.3104749426792832" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A = [1. 2; 3 4; 5 6]\n", + "A = qrhous!(A)\n", + "A * [1., 2, 3]" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "afd34800-7b9a-4006-9849-071fa197d962", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 2.8498844747360064\n", + " -2.404529848146404\n", + " -0.31047494267928366" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using LinearAlgebra: qr!\n", + "A = [1. 2.; 3. 4.; 5. 6.]\n", + "(Q, R) = qr!(A)\n", + "Q * [1., 2, 3]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.3", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/project/thinQR.jl b/project/thinQR.jl new file mode 100644 index 0000000..3df71ec --- /dev/null +++ b/project/thinQR.jl @@ -0,0 +1,94 @@ +module thinQR + +import Base: size, show, getproperty, getfield, propertynames, * +using LinearAlgebra: norm, I, triu, diagm + +export QRthin, qrhous!, qyhous + +mutable struct QRthin{T <: Real} + A::AbstractVecOrMat{T} + d::AbstractArray{T} + AQ::Union{Nothing, AbstractVecOrMat{T}} + AR::Union{Nothing, AbstractVecOrMat{T}} + + QRthin(A, d, AQ=nothing, AR=nothing) = new{eltype(A)}(A, d, AQ, AR) +end + +function qrhous!(A::AbstractMatrix{T})::QRthin{T} where T + m, n = size(A) + + d = zeros(n) + for j ∈ 1:n + s = norm(A[j:m, j]) + # iszero(s) && throw(ArgumentError("The matrix A is singular")) + 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 + A[j:m, j+1:n] -= A[j:m, j] * (A[j:m, j]' * A[j:m, j+1:n]) + end + end + + return QRthin(A, d) +end + +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]) + end + return z +end + +function calculateQ(A::QRthin{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::QRthin{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::QRthin{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::QRthin{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::QRthin, private::Bool=false) = (:R, :Q, (private ? fieldnames(typeof(A)) : ())...) + +Base.size(A::QRthin) = size(getfield(A, :A)) + +function *(A::QRthin{T}, b::AbstractVector{T}) where {T} + return qyhous(A, b) +end + +end # module \ No newline at end of file