Merge pull request #1558 from isaacsas/return_pivot_cols

optionally return col ordering in nullspace

Author: YingboMa

src/structural_transformation/bareiss.jl

@@ -207,11 +207,22 @@ function bareiss!(M::AbstractMatrix{T}, swap_strategy=bareiss_colswap;
     return (n, pivot, column_permuted)
 end
 
-function nullspace(A)
+function nullspace(A; col_order=nothing)
     column_pivots = collect(1:size(A, 2))
     B = copy(A)
     (rank, d, column_permuted) = bareiss!(B; column_pivots)
     reduce_echelon!(B, rank, d)
+
+    # The first rank entries in col_order are columns that give a basis
+    # for the column space. The remainder give the free variables.
+    if col_order !== nothing
+        resize!(col_order, size(A,2))
+        col_order .= 1:size(A,2)
+        for (i,cp) in enumerate(column_pivots)
+            @swap(col_order[i],col_order[cp])
+        end
+    end
+    
     N = ModelingToolkit.reduced_echelon_nullspace(rank, B)
     apply_inv_pivot_rows!(N, column_pivots)
 end

test/linalg.jl

@@ -18,3 +18,15 @@ N = ModelingToolkit.nullspace(A)
 @test size(N, 2) == 3
 @test rank(N) == 3
 @test iszero(A * N)
+
+A = [0 1 2 0 1 0; 
+     0 0 0 0 0 1;
+     0 0 0 0 0 1;
+     1 0 1 2 0 1;
+     0 0 0 2 1 0]
+col_order = Int[]
+N = ModelingToolkit.nullspace(A; col_order)
+colspan = A[:,col_order[1:4]]   # rank is 4
+@test iszero(ModelingToolkit.nullspace(colspan))
+@test !iszero(ModelingToolkit.nullspace(A[:,col_order[1:5]]))
+@test !iszero(ModelingToolkit.nullspace(A[:,[col_order[1:4]...,col_order[6]]]))