/-Operator for LU

[judober]

PR to add the /-Operator for LU-Types refering to this.
I'm fairly new to julia, so it is likely that there is a better way to calculate pivot (line 146/147 of lu.jl).
The PR also includes a testset for the division with LU.
Feel free to suggest improvements.

Fixes #582

[coveralls]

Coverage increased (+0.03%) to 81.614% when pulling c1325ba on judober:LuForRdiv into 1c53f4a on JuliaArrays:master.

[andyferris]

Will this allocation be elided?

[judober]

By elided you mean omitted? (Sorry, I'm not familiar with the technical terms).
If that's the case, than I think no. It is needed for the next step to be accessed at the indices F.p.
The problem is that the pivot indices cannot be used the same was as in \. I wrote down the whole operation with matrices and the difference is that in case of / the P-matrix is multiplied from the right, while it is normally from the left (as is case of \). So the p-indices (F.p) have to be ordered differently. There might exist a smarter way to do this but I don't know it.

[c42f]

Yes, in this context, "elided" means "optimized away".

This operation with the indices looks like Base.invperm(). I'd suggest adding the following utility function:

Base.@propagate_inbounds function Base.invperm(p::StaticVector)
     ip = similar(p)
     ip[p] = 1:length(p)
     similar_type(p)(ip)
end

Checking with @code_native indicates that the MVector which is created here by similar is optimized away :-D And so is the range object, which is very cool and kind of surprising.

Then inside your function, use pivot = @inbounds invperm(F.p) (we can use @inbounds because we know that F.p is a permutation).

[judober]

Thanks, I was not aware of this function.
I added your suggestion to util. The performance increases a bit by using @inbounds.
However, in Base invperm also checks if p is a valid pivot vector (by checking if every value is unique). This might be worth implementing though it probably adds computation time.

[andyferris]

I think we just have to assume the data in the LU is correct (that the pivot vector is a valid permutation).

[andyferris]

@c42f I'm not sure when to "trust" the optimizer - my understanding is that once the function is inlined into a more complex one, these optimizations may or may not occur (I'd love to know what is guaranteed...).

[andyferris]

I think you need to make this just an internal-facing static arrays utility function (not an overload of Base.invperm because I feel that one should check the permutation is valid).

[judober]

I found the time to remove the overloading of Base.invperm. Not sure why the checks are failing though. Probably due to incompatibility to code changes from the last months?

[c42f]

The test errors look like they might be real. At least, they're related to the LU tests and fail on all julia versions.

[c42f]

I'd love to merge this but we need to understand / fix the test failures. Do you have any time to investigate?

[judober]

I found some time and I think I know the problem and possibly a solution:

• The error occurs only when the size of A is 2x2 or 3x3. This is related to the special case treatment in solve.jl for these kinds of matrices. This treatment is only applied, when lu() is not used und has slightly different results.
• When I change the test from int to float, it passes. I think that this is due to the different value range of rand() for int and float.

So for now I changed the type of the test to float. Maybe a different kind of test would be advisable?

[judober]

Just curious: is there something left to do?
I think the difference in the values of int vs float is a problem for \approx. It does only check for the absolute error which increases for big values, while the (in this case more relevant) relative error stays as small as expected.
So I think this ready.

[c42f]

Thanks for bumping this, I must have missed your previous update.

It does only check for the absolute error which increases for big values, while the (in this case more relevant) relative error stays as small as expected.

Actually isapprox uses relative tolerance rtol=√eps and absolute tolerance atol=0 by default. It could be something to do with the distribution of random integer matrices being poorly behaved in some way though I find this a bit surprising.

It would be nice to understand why that was failing but it's hard to imagine it's related to this PR.