added project folder
This commit is contained in:
94
project/thinQR.jl
Normal file
94
project/thinQR.jl
Normal file
@ -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
|
||||
Reference in New Issue
Block a user