[I am reposting this from julia-dev, where it was suggested to try here.]
Dear all:
I am seeing a curious problem when assigning submatrices, e.g., in
R[k:m, :] = P * R[k:m, :]
I was implementing QR factorization and here is the code.
==========
function bug(A)
(m, n) = size(A)
Q = eye(m)
R = copy(A)
for k in 1:n
z = R[k:m, k]
v = [-sign(z[1]) * norm(z) - z[1]; -z[2:end]]
P = eye(m-k+1) - 2 * (v * v') / dot(v, v) # does not work
##P = eye(m-k+1) - (v * v') # does not work
##P = eye(m-k+1) # works
##P = rand(m-k+1, m-k+1) # does not work
R[k:m, :] = P * R[k:m, :]
Q[k:m, :] = P * Q[k:m, :]
end
(Q', triu(R))
end
==========
This is the error I get:
==========
ERROR: InexactError()
in setindex! at multidimensional.jl:76
in bug at [...].jl:75
in include at ./boot.jl:245
in include_from_node1 at ./loading.jl:128
in eval at [...]/src/julia-0.3/base/sysimg.jl:7
in anonymous at multi.jl:1310
in run_work_thunk at multi.jl:621
in run_work_thunk at multi.jl:630
in anonymous at task.jl:6
while loading [...].jl, in expression starting on line 103
==========
The offending line seems to be:
R[k:m, :] = P * R[k:m, :]
Curiously enough, the InexactError does not occur when P is the identity
matrix. Here are some data points:
P = eye(m-k+1) - 2 * (v * v') / dot(v, v) # does not work
P = eye(m-k+1) - (v * v') # does not work
P = eye(m-k+1) # works
P = rand(m-k+1, m-k+1) # does not work
I'd appreciate any insights... It may be a problem with setindex! in
multidimensional.jl.
Cheers,
Clemens
