It is easier to help you if you provide an example that can actually run.
What is A when you get the error.
A guess could be that it is something like a Matrix{Int} and therefore R =
copy(A) is also Matrix{Int}. Hence, the line R[k:m, :] = P * R[k:m, :]
fails as the right hand side consists of floating point values that are not
integer valued. It works when you use eye() because it is integer valued
even though it is a Matrix{Float64}.
Med venlig hilsen
Andreas Noack
2014-09-20 20:24 GMT-04:00 Clemens Heitzinger <[email protected]>
:
> [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
>
