Merge pull request #6 from PBrdng/master

RowEchelon_with_pivots also returns the pivot vector

src/RowEchelon.jl

@@ -79,4 +79,7 @@ rref(A::Matrix{T}) where {T <: Complex} = rrefconv(Complex128, A)
 rref(A::Matrix{Complex128}) = rref!(copy(A))
 rref(A::Matrix{T}) where {T <: Union{Integer, Float16, Float32}} = rrefconv(Float64, A)
 
+
+include("RowEchelon_with_pivots.jl")
+
 end # module

src/RowEchelon_with_pivots.jl

@@ -0,0 +1,93 @@
+
+# The function rref(), but with the modification that it also returns the pivot vector
+
+export rref_with_pivots, rref_with_pivots!
+
+function rref_with_pivots!(A::Matrix{T}, ɛ=T <: Union{Rational,Integer} ? 0 : eps(norm(A,Inf))) where T
+    nr, nc = size(A)
+    pivots = Vector{Int64}()
+    i = j = 1
+    while i <= nr && j <= nc
+        (m, mi) = findmax(abs.(A[i:nr,j]))
+        mi = mi+i - 1
+        if m <= ɛ
+            if ɛ > 0
+                A[i:nr,j] = 0
+            end
+            j += 1
+        else
+            for k=j:nc
+                A[i, k], A[mi, k] = A[mi, k], A[i, k]
+            end
+            d = A[i,j]
+            for k = j:nc
+                A[i,k] /= d
+            end
+            for k = 1:nr
+                if k != i
+                    d = A[k,j]
+                    for l = j:nc
+                        A[k,l] -= d*A[i,l]
+                    end
+                end
+            end
+            append!(pivots,j)
+            i += 1
+            j += 1
+        end
+    end
+    return A, pivots
+end
+
+rref_with_pivots_conv{T}(::Type{T}, A::Matrix) = rref_with_pivots!(copy!(similar(A, T), A))
+
+"""
+    rref_with_pivots(A)
+Compute the reduced row echelon form of the matrix A together with the
+position of the pivots.
+Since this algorithm is sensitive to numerical imprecision,
+* Complex numbers are converted to Complex128
+* Integer, Float16 and Float32 numbers are converted to Float64
+* Rational are kept unchanged
+
+```jldoctest
+julia> rref_with_pivots([ 1  2 -1  -4;
+              2  3 -1 -11;
+             -2  0 -3  22])
+3×4 Array{Float64,2}:
+ 1.0  0.0  0.0  -8.0
+ 0.0  1.0  0.0   1.0
+ 0.0  0.0  1.0  -2.0
+ Int64[3]:
+ 1
+ 2
+ 3
+
+julia> rref_with_pivots([16  2  3  13;
+              5 11 10   8;
+              9  7  6  12;
+              4 14 15   1])
+4×4 Array{Float64,2}:
+ 1.0  0.0  0.0   1.0
+ 0.0  1.0  0.0   3.0
+ 0.0  0.0  1.0  -3.0
+ 0.0  0.0  0.0   0.0
+ Int64[3]:
+ 1
+ 2
+ 3
+
+julia> rref_with_pivots([ 1  2  0   3;
+              2  4  0   7])
+2×4 Array{Float64,2}:
+ 1.0  2.0  0.0  0.0
+ 0.0  0.0  0.0  1.0
+ Int64[3]:
+ 1
+ 4
+```
+"""
+rref_with_pivots(A::Matrix{T}) where {T} = rref_with_pivots!(copy(A))
+rref_with_pivots(A::Matrix{T}) where {T <: Complex} = rref_with_pivots_conv(Complex128, A)
+rref_with_pivots(A::Matrix{Complex128}) = rref_with_pivots!(copy(A))
+rref_with_pivots(A::Matrix{T}) where {T <: Union{Integer, Float16, Float32}} = rref_with_pivots_conv(Float64, A)

test/runtests.jl

@@ -3,30 +3,35 @@ using Base.Test
 
 As = Vector{Matrix{Int}}()
 Rs = Vector{Matrix{Int}}()
+Ps = Vector{Vector{Int}}()
 # Issue #518 of Julia Base
 push!(As, [ 1  2  0   3;
             2  4  0   7])
 push!(Rs, [ 1  2  0   0;
             0  0  0   1])
+push!(Ps, [ 1; 4])
 # Example from Wikipedia
 push!(As, [ 1  3  -1;
             0  1   7])
 push!(Rs, [ 1  0 -22;
             0  1   7])
+push!(Ps, [ 1; 2])
 # Example from Wikibooks
 push!(As, [ 0  3 -6   6  4  -5;
             3 -7  8  -5  8   9;
             3 -9 12  -9  6  15])
 push!(Rs, [ 1  0 -2   3  0 -24;
             0  1 -2   2  0  -7;
             0  0  0   0  1   4])
+push!(Ps, [ 1; 2; 5])
 # Example from Rosetta Code
 push!(As, [ 1  2 -1  -4;
             2  3 -1 -11;
            -2  0 -3  22])
 push!(Rs, [ 1  0  0  -8;
             0  1  0   1;
             0  0  1  -2])
+push!(Ps, [ 1; 2; 3])
 # Magic Square
 push!(As, [16  2  3  13;
             5 11 10   8;
@@ -36,7 +41,7 @@ push!(Rs, [ 1  0  0   1;
             0  1  0   3;
             0  0  1  -3;
             0  0  0   0])
-
+push!(Ps, [ 1; 2; 3])
 
 @testset "Matrices of integers (treated as Float64)" begin
     for (A,R) in zip(As,Rs)
@@ -69,3 +74,42 @@ end
         @test C ≈ R
     end
 end
+
+@testset "Matrices of integers (treated as Float64) with Pivots" begin
+    for (A,X) in zip(As, zip(Rs,Ps))
+        R=X[1]
+        P=X[2]
+        C = rref_with_pivots(A)
+        @test C[1] isa Matrix{Float64}
+        @test C[1] ≈ R
+        @test C[2] == P
+    end
+end
+
+@testset "Matrices of rationals with Pivots" begin
+    for (A,X) in zip(As, zip(Rs,Ps))
+        R=X[1]
+        P=X[2]
+        B = Matrix{Rational{Int}}(A)
+        C = rref_with_pivots(B)
+        @test C[1] isa Matrix{Rational{Int}}
+        @test C[1] == R
+        @test C[2] == P
+    end
+end
+
+@testset "Matrix of Complex numbers with Pivots" begin
+    A = [1+ im 1-im  4;
+         3+2im 2-im  2]
+    R = [1     0    -2- im;
+         0     1     1+4im]
+    for conv in [false, true]
+        if conv
+            A = Matrix{Complex128}(A)
+        end
+        C = rref_with_pivots(A)
+        @test C[1] isa Matrix{Complex128}
+        @test C[1] ≈ R
+        @test C[2] == [1; 2]
+    end
+end