generic lufact requires more than it should

[arghhhh]

I've made a minimal Mod{N} with only the operations required to be a field. I needed to make a few modifications to the linalg code to get A \ b to work for a square matrix A with Mod{N} as the element type.

Specifically: In a few places, a literal "0" is used instead of zero( .. ). I'd rather not add a conversion / promotion. It should work without this. In one place "real( .. )" is used to initialize a variable that will hold the result of abs( .. )

If I turn off pivoting to reduce the requirements on my type further (so I don't need < or abs) then it doesn't work due to a SingularException error. I added a search for the first non-zero index and pivot on that - so its now pivoting, but without the requirement to have "abs" and "<". I use this to line to turn off pivoting for matrices of my type Mod{N}: import Base: lufact lufact{N}(A::AbstractMatrix{Mod{N}}, pivot::Union{Type{Val{false}}, Type{Val{true}}} = Val{false}) = lufact!(copy(A), pivot)

Below is a complete diff of my changes relative to v0.6.0-rc2.

Its pretty amazing that you can do this in Julia. It'd be great if we could fix these minor issues. I'm not a mathematician, so I'm not completely sure about the pivoting on non-zero element, but it worked for me on my one(!) test case.

diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl
index 6f44208..21e032e 100644
--- a/base/linalg/generic.jl
+++ b/base/linalg/generic.jl
@@ -955,7 +955,7 @@ true
 function istriu(A::AbstractMatrix)
     m, n = size(A)
     for j = 1:min(n,m-1), i = j+1:m
-        if A[i,j] != 0
+        if A[i,j] != zero( eltype(A) )
             return false
         end
     end
@@ -990,7 +990,7 @@ true
 function istril(A::AbstractMatrix)
     m, n = size(A)
     for j = 2:n, i = 1:min(j-1,m)
-        if A[i,j] != 0
+        if A[i,j] != zero( eltype(A) )
             return false
         end
     end
diff --git a/base/linalg/lu.jl b/base/linalg/lu.jl
index 7105bf3..ac2cf7f 100644
--- a/base/linalg/lu.jl
+++ b/base/linalg/lu.jl
@@ -39,7 +39,7 @@ function generic_lufact!(A::StridedMatrix{T}, ::Type{Val{Pivot}} = Val{true}) wh
             # find index max
             kp = k
             if Pivot
-                amax = real(zero(T))
+                amax = abs(zero(T))
                 for i = k:m
                     absi = abs(A[i,k])
                     if absi > amax
@@ -47,9 +47,16 @@ function generic_lufact!(A::StridedMatrix{T}, ::Type{Val{Pivot}} = Val{true}) wh
                         amax = absi
                     end
                 end
+        else # find first index with non-zero
+                for i = k:m
+                    if A[i,k] != zero(T)
+                        kp = i
+                        break
+                    end
+                end
             end
             ipiv[k] = kp
-            if A[kp,k] != 0
+            if A[kp,k] != zero(eltype(A))
                 if k != kp
                     # Interchange
                     for i = 1:n

[dpsanders]

Could you please submit a Pull Request with these changes so that people can easily test them out? It would also be great to include some tests, for example with your ModN type.

[arghhhh]

OK. I am new to this (never contributed to any open source project). There are really two issues here - one is relatively trivial uses of zero() and one() instead of 0 and 1 literals and the use of abs(zero(T)) instead of real(zero(T)) to avoid having to define real() which might not make sense for some matrix eltypes (eg for modulo{polynomial} ).

I put an example derived from the ModInt example in the Julia source here: https://gist.github.com/arghhhh/abdb8b6d28039f683f5287e124903670 and will submit a pull request shortly.

The second issue is how to control pivoting in lufact. It seems to me that there could be three options here:

1. • No pivoting
2. • Non-zero pivoting - ensure that the pivot is non-zero, but don't attempt to choose largest
3. • Pivoting - requires abs and < which might not make sense for some eltypes

The current implementation uses Val{false} and Val{true} to choose between 1. and 3. Adding 2. would require a bigger discussion.

[andreasnoack]

It would be great with the other pivoting option as well. The reason for the current design is that we try to use the same API for LU and QR. The pivoted QR returns a different type which is the reason why we use Vals here. For LU, we use the same type so Val is not really needed but we should probably be consistent across LU and QR. We could consider adding pivoting types instead or use symbols instead of Booleans in the Val, e.g. Pivot{:anyrow}, Val{:maxrow}, or NormCol. I'm not super happy with any of these. Thoughts?

[arghhhh]

I had wondered why the Val{} argument had been used. I was thinking that using a symbol instead of the Bool would work - but I only realized this needed to be a separate option after my original issue submission. My original proposal hijacked Pivot=Val{false} to mean pivot on !iszero.

I'm not sure how to proceed. I'm new to contributing to Julia, so it was exciting(!) to see my PR make progress - but I haven't actually solved my original problem. My motivation is that I'm learning about finite fields and abstract algebra to enable understanding and implementation of error control coding - similar to https://github.com/andrewcooke/IntModN.jl