added views decorator

This commit is contained in:
elvis
2023-11-25 20:17:20 +01:00
parent f36dec52dd
commit 5161464087
2 changed files with 1928 additions and 602 deletions

View File

@ -16,19 +16,7 @@ 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
function householder_vector(x::AbstractVecOrMat{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
@ -45,12 +33,12 @@ function qrfact(A::Matrix{T})::QRhous where T
R = deepcopy(A)
d = zeros(n)
for k 1:n
(u, s) = householder_vector(R[k:end, k])
for k 1:min(m, n)
@views (u, s) = householder_vector(R[k:m, 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
@views R[k:m, k+1:n] -= 2 * u * (u' * R[k:m, k+1:n])
R[k:m, k] .= u
end
return QRhous(R, d)
end
@ -59,19 +47,19 @@ end
function qyhous(A::QRhous{T}, y::AbstractArray{T}) where T
m, n = size(A)
z = deepcopy(y)
for j n:-1:1
for k 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])
@views z[k:m] -= 2 * A.A[k:m, k] .* dot(A.A[k:m, k], z[k:m])
end
return z
end
function qyhoust(A::QRhous{T}, y::AbstractArray{T}) where T
m, n = size(A.A)
m, n = size(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])
for k 1:n
# z[j:m] = z[j:m] - 2 * A.A[j:m, j] * (A.A[j:m, j]' * z[j:m])
@views z[k:m] -= 2 * A.A[k:m, k] .* dot(A.A[k:m, k], z[k:m])
end
return z
end
@ -125,9 +113,18 @@ 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
@show size(v)
@show size(x)
@show size(A.A)
@show size(A.d)
@show m
#@inbounds begin
for j m:-1:1
@views x[j] = (v[j] - dot(x[j+1:m], A.A[j, j+1:m])) * A.d[j]^-1
end
#end
return x
end